题解 P4705 【玩游戏】

题意

给出两个序列\(a,b\)，长度分别为\(n,m\)，请求出\(x=1,2,\ldots,t\)，下式的值

\[\dfrac{\sum_{i=1}^n\sum_{j=1}^m(a_i+b_j)^x}{nm} \]

题解

看到还能发就来水一篇

本人数学比较菜有错误请联系这个屑

我们记\(x\)的答案为\(ans_x\)，那么就有：

\[ans_x\times nm=\sum_{i=1}^n\sum_{j=1}^m(a_i+b_j)^x \]

好像也做不了什么，那就把\((a_i+b_j)^x\)用二项式定理展开再说：

\[ans_x\times nm=\sum_{i=1}^n\sum_{j=1}^m\sum_{k=0}^x\tbinom{x}{k}a_i^{k}b_j^{x-k} \]

交换顺序，再把万恶的组合数展开。

\[ans_x\times nm=\sum_{k=0}^x\sum_{i=1}^n\sum_{j=1}^m\tbinom{x}{k}a_i^{k}b_j^{x-k}=\sum_{k=0}^x\sum_{i=1}^n\sum_{j=1}^m\frac{x!}{k!(x-k)!}a_i^{k}b_j^{x-k} \]

\[\frac{ans_x\times nm}{x!}=\sum_{k=0}^x\frac{\sum_{i=1}^na_i^k}{k!}\frac{\sum_{j=1}^mb_j^{x-k}}{(x-k)!} \]

如果我们记\(A(k)=\sum_{i=1}^na_i^k\),\(B(k)=\sum_{i=1}^mb_i^k\)，那么就有：

\[\frac{ans_x\times nm}{x!}=\sum_{k=0}^x\frac{A(k)}{k!}\frac{B(x-k)}{(x-k)!} \]

这样卷积就很明显了。

于是剩下的问题就变成求\(A(0),A(1),A(2),\ldots,A(t)\)，考虑构造下面一个多项式：

\[F(x)=\prod_{i=1}^n(1+a_ix) \]

这个可以分治得到，复杂度显然正确

两边同取\(\ln\)得到：

\[\ln F(x)=\sum_{i=1}^n\ln(1+a_ix) \]

两边求导：

\[(\ln F(x))^\prime=\sum_{i=1}^n\frac{a_i}{1+a_ix} \]

又知，\(\frac{1}{1+a_ix}=\frac{1}{1-(-a_ix)}\)，即为一个首项为\(1\)，公比为\(-a_ix\)的数列的和，那么\(\frac{1}{1+a_ix}=\sum_{j=0}^{+\infty}(-a_ix)^j\)

带带带回去：

\[(\ln F(x))^\prime=\sum_{i=1}^na_i\sum_{j=0}^{+\infty}(-a_ix)^j \]

再积回去：

\[\ln F(x)=\int \sum_{i=1}^na_i\sum_{j=0}^{+\infty}(-a_ix)^jdx=\sum_{i=1}^na_i\sum_{j=1}^{+\infty}\frac{(-a_i)^{j-1}x^j}{j} \]

\[=\sum_{i=1}^n\sum_{j=1}^{+\infty}\frac{(-1)^{j-1}a_i^{j}x^j}{j}=\sum_{j=1}^{+\infty}\frac{(-1)^{j-1}x_j}{j}\sum_{i=1}^na_i^j=\sum_{j=1}^{+\infty}\frac{(-1)^{j-1}A(j)}{j}x^j \]

于是此时求出\(A\)也不难了，随便乘一下就好了。

\(B\)也同理。把\(A,B\)带回去，不难得到答案。

代码

#include<bits/stdc++.h>
namespace in{
    #ifdef slow
    inline int getc(){return getchar();}
    #else
    char buf[1<<21],*p1=buf,*p2=buf;
    inline int getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
    #endif
    template <typename T>inline void read(T& t){
		t=0;int f=0;char ch=getc();while (!isdigit(ch)){if(ch=='-')f = 1;ch=getc();}
	    while(isdigit(ch)){t=t*10+(ch-48);ch = getc();}if(f)t=-t;
	}
    template <typename T,typename... Args> inline void read(T& t, Args&... args){read(t);read(args...);}
}
namespace out{
	char buffer[1<<21];int p1=-1;const int p2 = (1<<21)-1;
	inline void flush(){fwrite(buffer,1,p1+1,stdout),p1=-1;}
	inline void putc(const char &x) {if(p1==p2)flush();buffer[++p1]=x;}
	template <typename T>void write(T x) {
		static char buf[15];static int len=-1;if(x>=0){do{buf[++len]=x%10+48,x/=10;}while (x);}else{putc('-');do {buf[++len]=-(x%10)+48,x/=10;}while(x);}
		while (len>=0)putc(buf[len]),--len;
	}
}
using namespace std;
template<const int mod>
struct modint{
    int x;
    modint<mod>(int o=0){x=o;}
    modint<mod> &operator = (int o){return x=o,*this;}
    modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
    modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
    modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
    modint<mod> &operator ^=(int b){
        modint<mod> a=*this,c=1;
        for(;b;b>>=1,a*=a)if(b&1)c*=a;
        return x=c.x,*this;
    }
    modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
    modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
    modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
    modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
    modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
	template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
    template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
    template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
    template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
    friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
    friend bool operator ==(modint<mod> a,int b){return a.x==b;}
    friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
    bool operator ! () {return !x;}
    modint<mod> operator - () {return x?mod-x:0;}
	modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;

const int mod=998244353;
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
	vector<modint<mod>>a;
	modint<mod>&operator[](int i){return a[i];}
	int size(){return a.size();}
	void resize(int n){a.resize(n);}
};
int rev[N];
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
	int n=1<<k;
	for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int mid=1;mid<n;mid<<=1){
		modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
		for(int r=mid<<1,j=0;j<n;j+=r){
			modint<mod> w=1;
			for(int k=0;k<mid;k++,w=w*wn){
				modint<mod> x=a[j+k],y=w*a[j+k+mid];
				a[j+k]=x+y,a[j+k+mid]=x-y;
			}
		}
	}
	if(typ<0){
		modint<mod> inv=modint<mod>(1)/n;
		for(int i=0;i<n;i++)a[i]*=inv;
	}
}
inline poly one(){poly a;a.a.push_back(1);return a;}
poly operator +(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
	int n=a.size()+b.size()-1,k=ext(n);
	a.resize(1<<k),b.resize(1<<k),init(k);
	ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
	ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }
poly inv(poly F,int k){
	int n=1<<k;F.resize(n);
	if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
	poly G,H=inv(F,k-1);
	G.resize(n),H.resize(n<<1),F.resize(n<<1);
	for(int i=0;i<n/2;i++)G[i]=H[i]*2;
	init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
	for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
	ntt(H,k+1,-1),H.resize(n);
	for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
	int n=a.size();
	a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly deriv(poly a){//求导 
	int n=a.size()-1;
	for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
	a.resize(n);return a;
}
inline poly inter(poly a){//求原 
	int n=a.size()+1;a.resize(n);
	for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
	a[0]=0;return a;
}
inline poly ln(poly a){
	int n=a.size();
	a=inter(deriv(a)*inv(a));
	a.resize(n);return a;
}
#define mid (l+r>>1)
inline poly F(int l,int r,int *a){
	poly G;
	if(l==r){G.resize(2);G[0]=1;G[1]=a[l];}
	else{
		G=F(l,mid,a)*F(mid+1,r,a);
	}
	//printf("[%d,%d] ",l,r);for(int i=0;i<G.size();i++)printf("%d ",G[i].x);puts("");
	return G;
}
int n,m,t,a[N],b[N];
poly A,B;
modint<mod>fac[N];
signed main(){
	fac[0]=1;for(int i=1;i<N;i++)fac[i]=fac[i-1]*i;
	in::read(n,m);
	for(int i=1;i<=n;i++)in::read(a[i]);
	for(int i=1;i<=m;i++)in::read(b[i]);
	in::read(t);
	A=F(1,n,a);A.resize(t+1);A=ln(A);
	A[0]=n;for(int i=1;i<=t;i++)if(i&1)A[i]*=i;else A[i]*=i,A[i]=-A[i];
	B=F(1,m,b);B.resize(t+1);B=ln(B);
	B[0]=m;for(int i=1;i<=t;i++)if(i&1)B[i]*=i;else B[i]*=i,B[i]=-B[i];
	//for(int i=0;i<=t;i++)cout<<A[i].x<<" ";cout<<endl;
	//for(int i=0;i<=t;i++)cout<<B[i].x<<" ";cout<<endl;
	for(int i=0;i<=t;i++)A[i]/=fac[i];
	for(int i=0;i<=t;i++)B[i]/=fac[i];
	poly ans=A*B;for(int i=0;i<=t;i++)ans[i]=(ans[i]*fac[i]/n)/m;
	for(int i=1;i<=t;i++)cout<<ans[i].x<<endl;
}