题解 P4705 【玩游戏】

题意

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

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

题解

看到还能发就来水一篇

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

我们记(x)的答案为(ans_x),那么就有:

[ans_x imes nm=sum_{i=1}^nsum_{j=1}^m(a_i+b_j)^x ]

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

[ans_x imes nm=sum_{i=1}^nsum_{j=1}^msum_{k=0}^x binom{x}{k}a_i^{k}b_j^{x-k} ]

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

[ans_x imes nm=sum_{k=0}^xsum_{i=1}^nsum_{j=1}^m binom{x}{k}a_i^{k}b_j^{x-k}=sum_{k=0}^xsum_{i=1}^nsum_{j=1}^mfrac{x!}{k!(x-k)!}a_i^{k}b_j^{x-k} ]

[frac{ans_x imes nm}{x!}=sum_{k=0}^xfrac{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 imes nm}{x!}=sum_{k=0}^xfrac{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}^nln(1+a_ix) ]

两边求导:

[(ln F(x))^prime=sum_{i=1}^nfrac{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_isum_{j=0}^{+infty}(-a_ix)^j ]

再积回去:

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

[=sum_{i=1}^nsum_{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;
}
原文地址:https://www.cnblogs.com/juruo-cjl/p/14220920.html