洛谷P4705 玩游戏（生成函数+多项式运算）

题面

传送门

题解

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

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

\[\begin{aligned} [x^k]Ans &={1\over nm}\sum_{i=1}^n\sum_{j=1}^m\left(a_i+b_j\right)^k\\ &={1\over nm}\sum_{i=1}^n\sum_{j=1}^m\sum_{p=0}^k{k\choose p}{a_i}^p{b_j}^{k-p}\\ &={k!\over nm}\sum_{p=0}^k{\sum_{i=1}^n{a_i}^p\over p!}{\sum_{j=1}^m{b_j}^{k-p}\over (k-p)!}\\ \end{aligned} \]

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

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

然后继续推倒

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

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

我们设

\[\begin{aligned} G(x) &=\sum_{i=1}^n{-a_i\over 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\)啊……

那就继续推倒

\[\begin{aligned} G(x) &=\sum_{i=1}^n{-a_i\over 1-a_ix}\\ &=\sum_{i=1}^n\ln'\left(1-a_ix\right)\\ &=\ln'\left(\prod_{i=1}^n (1-a_ix)\right) \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]='\n';
}
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;
}