玩游戏

XXV.玩游戏

我们考虑令\(f(p)\)表示游戏的“\(p\)次价值”的期望。

则按照期望定义,我们有

\[f(p)=\dfrac{\sum\limits_{i=0}^n\sum\limits_{j=0}^m(a_i+b_j)^p}{nm} \]

考虑二项式暴力展开,得到

\[f(p)=\dfrac{\sum\limits_{i=0}^n\sum\limits_{j=0}^m\sum\limits_{k=0}^pa_i^kb_j^{p-k}\dbinom{p}{k}}{nm} \]

于是我们拆开组合数并调换枚举顺序,得到

\[f(p)=\dfrac{p!}{nm}\sum\limits_{k=0}^p\Big(\sum\limits_{i=0}^n\dfrac{a_i^k}{k!}\Big)\Big(\sum\limits_{j=0}^m\dfrac{b_j^{p-k}}{(p-k)!}\Big) \]

则我们考虑预处理

\[A_k=\sum\limits_{i=0}^na_i^k,B_k=\sum\limits_{j=0}^mb_j^k \]

则就有

\[f(p)=\dfrac{p!}{nm}\sum\limits_{k=0}^p\dfrac{A_k}{k!}\dfrac{B_{p-k}}{(p-k)!} \]

卷积的形式,可以直接NTT求出。

我们现在考虑\(A\)的OGF

\[A(x)=\sum\limits_{i=0}^{\infty}x^i\sum\limits_{j=1}^na_j^i \]

调换枚举顺序就是

\[A(x)=\sum\limits_{j=1}^n\sum\limits_{i=0}^{\infty}a_j^ix^i=\sum\limits_{i=1}^n\dfrac{1}{1-a_ix} \]

直接求右边那玩意会挂掉;我们考虑通过求\(\ln\)的方法来计算。

我们有

\[\Big(\ln(1-ax)\Big)'=\dfrac{-a}{1-ax} \]

所以有

\[\dfrac{1}{1-ax}=1+\dfrac{ax}{1-ax}=1-x\Big(\ln(1-ax)\Big)' \]

故而

\[A(x)=\sum\limits_{i=1}^n1-x\Big(\ln(1-a_ix)\Big)'=n-x\sum\limits_{i=1}^n\Big(\ln(1-a_ix)\Big)' \]

因为求导具有线性性(即和的导等于导的和),故我们有

\[A(x)=n-x\Big(\sum\limits_{i=1}^n\ln(1-a_ix)\Big)' \]

我们用\(\ln\)化加为乘,得到

\[A(x)=n-x\Bigg(\ln\Big(\prod\limits_{i=1}^n1-a_ix\Big)\Bigg)' \]

完成!则现在唯一的问题在于求出

\[\prod\limits_{i=1}^n1-a_ix \]

它可以使用分治求出(类比分治NTT)

故复杂度为\(O(n\log^2n)\)。

代码(要开O2):

#include<bits/stdc++.h>
using namespace std;
const int N=1<<20;
const int mod=998244353;
int ksm(int x,int y){
	int rt=1;
	for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)rt=(1ll*rt*x)%mod;
	return rt;
}
namespace Poly{
	const int G=3;
	int rev[N];
	void NTT(int *a,int tp,int LG){
		int lim=(1<<LG),invlim=ksm(lim,mod-2);
		for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
		for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int md=1;md<lim;md<<=1){
			int rt=ksm(G,(mod-1)/(md<<1));
			if(tp==-1)rt=ksm(rt,mod-2);
			for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
				int w=1;
				for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
					int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
					a[pos+i]=(x+y)%mod;
					a[pos+md+i]=(x-y+mod)%mod;
				}
			}
		}
		if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
	}
	int A[N],B[N],C[N];
	void mul(int *a,int *b,int *c,int LG){//using: Array A and B
		int lim=(1<<LG);
		for(int i=0;i<lim;i++)A[i]=B[i]=0;
		for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
		NTT(A,1,LG),NTT(B,1,LG);
		for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
		NTT(A,-1,LG);
		for(int i=0;i<lim;i++)c[i]=A[i];
	}
	void inv(int *a,int *b,int LG){//using: Array C
		b[0]=ksm(a[0],mod-2);
		for(int k=1;k<=LG+1;k++){
			mul(b,a,C,k);
			for(int i=0;i<(1<<k);i++)C[i]=(mod-C[i])%mod;
			(C[0]+=2)%=mod;
			mul(C,b,b,k);
		}
	}
	void diff(int *a,int *b,int lim){
		for(int i=0;i<lim;i++)b[i]=1ll*a[i+1]*(i+1)%mod;
		b[lim-1]=0;
	}
	void ln(int *a,int *b,int LG){//using: Array C
		inv(a,b,LG);
		diff(a,C,1<<LG);
		mul(b,C,b,LG+1);
	}
}
void solve(int *a,int *A,int LG){
	if(!LG){A[0]=1,A[1]=(mod-a[0])%mod;return;}
	solve(a,A,LG-1);
	solve(a+(1<<(LG-1)),A+(2<<(LG-1)),LG-1);
	Poly::mul(A,A+(2<<(LG-1)),A,LG+1);
}
int tmp[N],p;
int fac[N],inv[N];
void func(int *a,int *A,int n){
	int all=0;
	while((1<<all)<n)all++;
	solve(a,A,all);
//	for(int i=0;i<(4<<all);i++)printf("%d ",A[i]);puts("");
	all=0;
	while((1<<all)<=p)all++;
	Poly::ln(A,tmp,all);
	A[0]=n;
	for(int i=1;i<=p;i++)A[i]=(mod-tmp[i-1])%mod;
	for(int i=0;i<=p;i++)A[i]=1ll*A[i]*inv[i]%mod;
}
int n,m;
int a[N],b[N],A[N],B[N];
int main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<n;i++)scanf("%d",&a[i]);
	for(int i=0;i<m;i++)scanf("%d",&b[i]);
	scanf("%d",&p);
	fac[0]=1;
	for(int i=1;i<=p;i++)fac[i]=1ll*fac[i-1]*i%mod;
	inv[p]=ksm(fac[p],mod-2);
	for(int i=p-1;i>=0;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod;
	func(a,A,n),func(b,B,m);
	int all=0;
	while((1<<all)<=p)all++;
	Poly::mul(A,B,tmp,all+1);
	int INV=ksm(1ll*n*m%mod,mod-2);
	for(int i=1;i<=p;i++)printf("%d\n",1ll*tmp[i]*fac[i]%mod*INV%mod);
	return 0;
} 

上一篇:LJJ爱数数


下一篇:【高等数学】第 1 讲 元素和极限