zoukankan      html  css  js  c++  java
  • 题解 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;
    }
    
  • 相关阅读:
    latex How do I know what symbols/characters are available in a font package
    fun字形
    inspection tool
    msys2 安装 基本配置
    R语言 测试 训练步骤
    R 语言描述性 数据分析 步骤
    如何构建分类模型
    史上最酷的数学动态图
    极大似然估计四个步骤
    react 组件间参数传递
  • 原文地址:https://www.cnblogs.com/juruo-cjl/p/14220920.html
Copyright © 2011-2022 走看看