• 洛谷P4705 玩游戏(生成函数+多项式运算)


    题面

    传送门

    题解

    妈呀这辣鸡题目调了我整整三天……最后发现竟然是因为分治(NTT)之后的多项式长度不是(2)的幂导致把多项式的值存下来的时候发生了一些玄学错误……玄学到了我(WA)的点全都是(WA)(2)的幂次行里……

    看到这种题目二话不说先推倒

    [egin{aligned} [x^k]Ans &={1over nm}sum_{i=1}^nsum_{j=1}^mleft(a_i+b_j ight)^k\ &={1over nm}sum_{i=1}^nsum_{j=1}^msum_{p=0}^k{kchoose p}{a_i}^p{b_j}^{k-p}\ &={k!over nm}sum_{p=0}^k{sum_{i=1}^n{a_i}^pover p!}{sum_{j=1}^m{b_j}^{k-p}over (k-p)!}\ end{aligned} ]

    然后这就被画成了一个卷积的形式

    定义两个多项式(A(x)=sum_{i=0}^infty x^isum_{j=1}^n{a_j}^i),和(B(x)=sum_{i=0}^infty x^isum_{j=1}^m{b_j}^i),只要我们能求出这两个多项式的系数,然后一通乱搞之后就能求出(Ans)

    然后继续推倒

    [egin{aligned} A(x) &=sum_{i=0}^infty x^isum_{j=1}^n{a_j}^i\ &=sum_{j=1}^nsum_{i=0}^infty {a_j}^ix^i\ &=sum_{i=1}^n{1over 1-a_ix}\ &=sum_{i=1}^n {a_i}^0+{a_i}^1x^1+{a_i}^2x^2+... end{aligned} ]

    所以……这玩意儿该咋算啊……

    我们设

    [egin{aligned} G(x) &=sum_{i=1}^n{-a_iover 1-a_ix}\ &=sum_{i=1}^n-{a_i}^1-{a_i}^2x-{a_i}^3x^2-...\ end{aligned} ]

    那么就有(A(x)=-xG(x)+n)

    然而我还是不会算(G)啊……

    那就继续推倒

    [egin{aligned} G(x) &=sum_{i=1}^n{-a_iover 1-a_ix}\ &=sum_{i=1}^nln'left(1-a_ix ight)\ &=ln'left(prod_{i=1}^n (1-a_ix) ight) end{aligned} ]

    分治(NTT)就行啦

    然后没有然后了

    我错了多项式比计算几何难调多了

    //minamoto
    #include<bits/stdc++.h>
    #define R register
    #define fp(i,a,b) for(R int i=(a),I=(b)+1;i<I;++i)
    #define fd(i,a,b) for(R int i=(a),I=(b)-1;i>I;--i)
    #define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
    using namespace std;
    char buf[1<<21],*p1=buf,*p2=buf;
    inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
    int read(){
        R int res,f=1;R char ch;
        while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
        for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
        return res*f;
    }
    char sr[1<<21],z[20];int C=-1,Z=0;
    inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
    void print(R int x){
        if(C>1<<20)Ot();if(x<0)sr[++C]='-',x=-x;
        while(z[++Z]=x%10+48,x/=10);
        while(sr[++C]=z[Z],--Z);sr[++C]='
    ';
    }
    const int N=(1<<18)+5,P=998244353,Gi=332748118;
    inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
    inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
    inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
    int ksm(R int x,R int y){
        R int res=1;
        for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
        return res;
    }
    vector<int>r[21];int rt[2][N<<1],inv[N],fac[N],ifac[N],lim,d;
    inline void init(R int len){lim=1,d=0;while(lim<len)lim<<=1,++d;}
    void Pre(){
        fp(d,1,18){
        	r[d].resize(1<<d);
        	fp(i,1,(1<<d)-1)r[d][i]=(r[d][i>>1]>>1)|((i&1)<<(d-1));
    	}
        inv[0]=inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
        fp(i,2,262144){
            fac[i]=mul(fac[i-1],i),
            inv[i]=mul(P-P/i,inv[P%i]),
            ifac[i]=mul(ifac[i-1],inv[i]);
        }
        for(R int t=(P-1)>>1,i=1,x,y;i<=262144;i<<=1,t>>=1){
            x=ksm(3,t),y=ksm(Gi,t);
            rt[1][i]=rt[0][i]=1;
            fp(k,1,i-1){
                rt[1][k+i]=mul(rt[1][k+i-1],x),
                rt[0][k+i]=mul(rt[0][k+i-1],y);
            }
        }
    }
    int rev[N];
    void NTT(int *A,int ty){
        fp(i,0,lim-1)if(i<r[d][i])swap(A[i],A[r[d][i]]);
        for(R int mid=1;mid<lim;mid<<=1)
            for(R int j=0;j<lim;j+=(mid<<1))
                for(R int k=0,t;k<mid;++k)
                    A[j+k+mid]=dec(A[j+k],t=mul(rt[ty][mid+k],A[j+k+mid])),
                    A[j+k]=add(A[j+k],t);
        if(!ty)for(R int i=0,inv=ksm(lim,P-2);i<lim;++i)A[i]=mul(A[i],inv);
    }
    void Inv(int *a,int *b,int len){
        if(len==1)return b[0]=ksm(a[0],P-2),void();
        Inv(a,b,len>>1),init(len<<1);
        static int A[N],B[N];
        fp(i,0,len-1)A[i]=a[i],B[i]=b[i];
        fp(i,len,lim-1)A[i]=B[i]=0;
        NTT(A,1),NTT(B,1);
        fp(i,0,lim-1)A[i]=mul(A[i],mul(B[i],B[i]));
        NTT(A,0);
        fp(i,0,len-1)b[i]=dec(add(b[i],b[i]),A[i]);
        fp(i,len,lim-1)b[i]=0;
    }
    void Ln(int *a,int *b,int len){
        static int A[N],B[N];
        fp(i,1,len-1)A[i-1]=mul(a[i],i);A[len-1]=0;
        Inv(a,B,len),init(len<<1);fp(i,len,lim-1)A[i]=B[i]=0;
        NTT(A,1),NTT(B,1);
        fp(i,0,lim-1)A[i]=mul(A[i],B[i]);
        NTT(A,0);
        fp(i,1,len-1)b[i]=mul(A[i-1],inv[i]);b[0]=0;
        fp(i,len,lim-1)b[i]=0;
    }
    int D[25][N];
    void solve(int *a,int d,int l,int r){
        if(l==r)return D[d][0]=1,D[d][1]=P-a[l],void();
        int mid=(l+r)>>1;
        solve(a,d,l,mid),solve(a,d+1,mid+1,r),init(r-l+1+1);
        static int A[N],B[N];
        fp(i,0,mid-l+1)A[i]=D[d][i];fp(i,mid-l+2,lim-1)A[i]=0;
        fp(i,0,r-mid)B[i]=D[d+1][i];fp(i,r-mid+1,lim-1)B[i]=0;
        NTT(A,1),NTT(B,1);
        fp(i,0,lim-1)A[i]=mul(A[i],B[i]);
        NTT(A,0);
        fp(i,0,r-l+1)D[d][i]=A[i];
        fp(i,r-l+2,lim-1)D[d][i]=0;
    }
    int a[N],b[N],ak[N],bk[N];
    int n,m,t;
    void calc(int *a,int *b,int n,int t){
        static int A[N],B[N];
        solve(a,1,1,n);
        init(t+1);int len=lim;
        fp(i,0,n)A[i]=D[1][i];
        fp(i,n+1,len-1)A[i]=0;
        Ln(A,B,len);
        fp(i,1,len-1)B[i-1]=mul(B[i],i);B[len-1]=0;
        b[0]=n;
        fp(i,1,t)b[i]=mul(P-B[i-1],ifac[i]);
    }
    void Mul(int *a,int *b){
        init(t<<1);
        NTT(a,1),NTT(b,1);
        fp(i,0,lim-1)a[i]=mul(a[i],b[i]);
        NTT(a,0);
        int invm=ksm(mul(n,m),P-2);
        fp(i,1,t)print(1ll*a[i]*fac[i]%P*invm%P);
    }
    int main(){
    //	freopen("testdata.in","r",stdin);
        Pre();
        n=read(),m=read();
        fp(i,1,n)a[i]=read();
        fp(i,1,m)b[i]=read();
        t=read();
        calc(a,ak,n,t),calc(b,bk,m,t);
        Mul(ak,bk);
        return Ot(),0;
    }
    
  • 相关阅读:
    12 Source Code Profilers for C & C++
    HttpWebRequest的使用方法
    MSDN Windows 下载
    Qt 4.7 在VS2010环境下的编译
    [转].NET Logging Tools and Libraries
    硬盘崩溃之后
    .net core 下使用 logdashboard 日志面板
    工具收藏 年终工作总结必备工具之ppt利器
    Dapper 的应用和Dapper.Contrib 的方法封装(一)
    Dapper 的应用和Dapper.Contrib 的方法封装(二)
  • 原文地址:https://www.cnblogs.com/bztMinamoto/p/10548277.html
Copyright © 2020-2023  润新知