CF621641A Majsoul

https://codeforces.com/group/MIyYz3rj9b/contest/621641/problem/A

如果 $a+b=n$ 我们知道

\[\sum_{x+y=m}\frac{\binom ax\binom by}{\binom nm}xy=\frac{abm^{\underline2} }{n^{\underline2} }\]

证明：

\[\begin{aligned} \sum_{x+y=m}\frac{\binom ax\binom by}{\binom nm}xy&=\sum_{x+y=m}\frac{\binom{a-1}{x-1}\binom{b-1}{y-1} }{\binom nm}ab\\ &=\frac{\binom{a+b-2}{x+y-2} }{\binom nm}ab\\ &=\frac{\binom{n-2}{m-2} }{\binom nm}ab\\ &=\frac{abm^{\underline2} }{n^{\underline2} } \end{aligned}\]

假设四个人分别摸到了 $a,b,c,d$ 张牌，分别有 $x,y,z,w$ 张是宝牌。
如果 $a+b+c+d=n$，上述式子的四元形式就是

\[\mathbb E[xyzw]=\sum_{x+y+z+w=m}\frac{\binom ax\binom by\binom cz\binom dw }{\binom nm}xyzw=\frac{abcdm^{\underline4} }{n^{\underline4 } }\]

证明类似，略。

因此 $m$ 的限制实际上没用，我们只需要知道 $\mathbb E[abcd]$ 的期望，乘以 $m^{\underline4}/n^{\underline4}$ 即得 $\mathbb E[xyzw]$。

在期望中，只有独立变量具有可乘性，但这里 $x,y,z,w$ 显然不独立。
但注意到 $x,y,z,w$ 每次操作至多只有 $1$ 个会加一，所以能把最后的乘法变成每次操作进行加法。

\[(x+1)yzw=xyzw+yzw\\ (x+1)yz=xyz+yz\\ (x+1)y=xy+y\]

因此只需设 $f_{i,j,k}$ 表示第 $i$ 次操作（摸牌）是由 $k$ 来完成的，集合 $j$ 内的乘积的期望与到达该种情况的概率的乘积即可。

#include<bits/stdc++.h>
#define siz(x) int((x).size())
#define all(x) std::begin(x),std::end(x)
#define fi first
#define se second
using namespace std;
using unt=unsigned;
using loli=long long;
using lolu=unsigned long long;
using pii=pair<int,int>;
mt19937_64 rng(random_device{}());
constexpr int P=998244353;
struct mint{
	int d;
	mint()=default;
	mint(int x):d(x){}
	friend std::istream&operator>>(std::istream&x,mint&y){return x>>y.d;}
	friend std::ostream&operator<<(std::ostream&x,mint y){return x<<y.d;}
	friend mint operator+(mint x,mint y){return (x.d+=y.d)<P?x.d:x.d-P;}
	mint&operator+=(mint z){return (d+=z.d)<P?d:d-=P,*this;}
	friend mint operator-(mint x,mint y){return (x.d-=y.d)<0?x.d+P:x.d;}
	mint&operator-=(mint z){return (d-=z.d)<0?d+=P:d,*this;}
	friend mint operator*(mint x,mint y){return int(1ll*x.d*y.d%P);}
	mint&operator*=(mint z){return d=int(1ll*d*z.d%P),*this;}
	static mint qpow(int x,int y=P-2){int z=1;for(;y;y>>=1,x=int(1ll*x*x%P))if(y&1)z=int(1ll*x*z%P);return z;}
	friend mint operator/(mint x,mint y){return x*=qpow(y.d);}
	mint&operator/=(mint z){return (*this)*=qpow(z.d);}
	friend mint operator^(mint x,mint y){return qpow(x.d,y.d);}
	mint&operator^=(mint z){return *this=qpow(d,z.d);}
	mint operator()(mint z)const{return qpow(d,z.d);}
	mint&operator[](mint z){return *this=qpow(d,z.d);}
	mint inv()const{return qpow(d);}
	mint pow(mint z)const{return qpow(d,z.d);}
	int operator+()const{return d;}
	mint operator-()const{return d?P-d:0;}
	int operator~()const{return ~d;}
};
mint operator""_m(lolu x){return mint(int(x%P));}
constexpr int N=1007;
int n,m;
mint f[N][1<<4][4],p[4][5];
signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	ios::sync_with_stdio(false);cin.tie(nullptr);
	cin>>n>>m;
	for(int i=0;i<4;i++){
		p[i][4]=1;
		for(int j=0;j<4;j++)
			cin>>p[i][(j+1)%4],p[i][4]-=(p[i][(j+1)%4]/=100);
	}
	f[1][1][0]=f[1][0][0]=1;
	for(int i=1;i<n;i++)for(int k=0;k<(1<<4);k++){
		for(int j1=0;j1<4;j1++){
			for(int j2=0;j2<4;j2++)
				if(((1<<j2)|k)==k)
					f[i+1][k][j2]+=(f[i][k][j1]+f[i][k^(1<<j2)][j1])*p[j1][j2];
				else
					f[i+1][k][j2]+=f[i][k][j1]*p[j1][j2];
			int j2=(j1+1)%4;
			if(((1<<j2)|k)==k)
				f[i+1][k][j2]+=(f[i][k][j1]+f[i][k^(1<<j2)][j1])*p[j1][4];
			else
				f[i+1][k][j2]+=f[i][k][j1]*p[j1][4];
		}
	}
	mint ans=0;
	for(int i=0;i<4;i++)
		ans+=f[n][15][i];
	cout<<ans*m*(m-1)*(m-2)*(m-3)/(n*(n-1)*(n-2)*(n-3));
	return 0;
}

要跑 $10^{18}$ 只需矩阵快速幂优化该 dp 转移。

#include<bits/stdc++.h>
#define siz(x) int((x).size())
#define all(x) std::begin(x),std::end(x)
#define fi first
#define se second
using namespace std;
using unt=unsigned;
using loli=long long;
using lolu=unsigned long long;
using pii=pair<int,int>;
mt19937_64 rng(random_device{}());
constexpr int P=998244353;
struct mint{
	int d;
	mint()=default;
	mint(int x):d(x){}
	friend istream&operator>>(istream&x,mint&y){return x>>y.d;}
	friend ostream&operator<<(ostream&x,mint y){return x<<y.d;}
	friend mint operator+(mint x,mint y){return (x.d+=y.d)<P?x.d:x.d-P;}
	mint&operator+=(mint z){return (d+=z.d)<P?d:d-=P,*this;}
	friend mint operator-(mint x,mint y){return (x.d-=y.d)<0?x.d+P:x.d;}
	mint&operator-=(mint z){return (d-=z.d)<0?d+=P:d,*this;}
	friend mint operator*(mint x,mint y){return int(1ll*x.d*y.d%P);}
	mint&operator*=(mint z){return d=int(1ll*d*z.d%P),*this;}
	static mint qpow(int x,int y=P-2){int z=1;for(;y;y>>=1,x=int(1ll*x*x%P))if(y&1)z=int(1ll*x*z%P);return z;}
	friend mint operator/(mint x,mint y){return x*=qpow(y.d);}
	mint&operator/=(mint z){return (*this)*=qpow(z.d);}
	friend mint operator^(mint x,mint y){return qpow(x.d,y.d);}
	mint&operator^=(mint z){return *this=qpow(d,z.d);}
	mint operator()(mint z)const{return qpow(d,z.d);}
	mint&operator[](mint z){return *this=qpow(d,z.d);}
	mint inv()const{return qpow(d);}
	mint pow(mint z)const{return qpow(d,z.d);}
	int operator+()const{return d;}
	mint operator-()const{return d?P-d:0;}
	int operator~()const{return ~d;}
};
mint operator""_m(lolu x){return mint(int(x%P));}
constexpr int N=107;
struct matrix:array<array<mint,N>,N>{
	matrix()=default;
	matrix(int x){
		memset(this,0,sizeof(*this));
		if(x==1)for(int i=0;i<N;i++)(*this)[i][i]=1;
	}
	friend matrix operator*(const matrix &x,const matrix &y){
		matrix z=0;
		for(int i=0;i<N;i++)for(int k=0;k<N;k++)for(int j=0;j<N;j++)
			z[i][j]+=x[i][k]*y[k][j];
		return z;
	}
	matrix&operator*=(const matrix&t){return *this=*this*t;}
	matrix&operator^=(loli y){
#define x (*this)
		matrix z=1;
		for(;y;y>>=1,x*=x)
			if(y&1)z*=x;
		return x=z;
#undef x
	}
}a,b;
loli n,m;
mint f[1<<4][4],p[4][4];
signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	ios::sync_with_stdio(false);cin.tie(nullptr);
	cin>>n>>m;
	for(int i=0;i<4;i++){
		for(int j=0;j<4;j++){
			cin>>p[i][(j+1)%4];
			p[i][(j+1)%4]/=100;
		}
		p[i][(i+1)%4]=1;
		for(int j=0;j<4;j++){
			if(i!=j)p[i][(i+1)%4]-=p[i][(j+1)%4];
		}
	}
	a[0][0]=a[0][1]=1;
	for(int k=0;k<16;k++)for(int j1=0;j1<4;j1++)for(int j2=0;j2<4;j2++){
		if(k>>j2&1)b[j1*16+k-(1<<j2)][j2*16+k]+=p[j1][j2];
		b[j1*16+k][j2*16+k]+=p[j1][j2];
	}
	a*=(b^=(n-1));
	mint ans=0;
	for(int i=0;i<4;i++)
		ans+=a[0][i*16+15];
	n%=P;m%=P;
	cout<<ans*m*(m-1_m)*(m-2_m)*(m-3_m)/(n*(n-1_m)*(n-2_m)*(n-3_m));
	return 0;
}