[Luogu5162]WD与积木(多项式求逆)

不要以为用上Stirling数就一定离正解更近，FFT都是从DP式本身出发的。

设f[i]为i个积木的所有方案的层数总和，g[i]为i个积木的方案数，则答案为$frac{f[i]}{g[i]}$

转移枚举第一层是哪些积木：$$f_n=g_n+sum_{i=1}^{n}inom{n}{i}f_{n-i},f_0=0$$$$g_n=sum_{i=1}^{n}inom{n}{i}g_{n-i},g_0=1$$

转化成卷积形式：$$frac{f_n}{n!}=frac{g_n}{n!}+sum_{i=1}^{n}frac{1}{i!} imes frac{f_{n-i}}{i!}$$$$frac{g_n}{n!}=sum_{i=1}^{n}frac{1}{i!} imes frac{g_{n-i}}{(n-i)!}$$

构造生成函数：$F(x)=sumfrac{f_i}{i!}x^i$,$G(x)=sumfrac{g_i}{i!}x^i$,$H(x)=sumfrac{x^i}{i!}$。

则根据上式有：$F(x)=G(x)+F(x)(H(x)-1)-1$,$G(x)=G(x)(H(x)-1)+1$。

移项得：$F(x)=frac{G(x)-1}{2-H(x)}=G(x)(G(x)-1)$,$G(x)=frac{1}{2-H(x)}$。

多项式求逆即可。

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=500010,mod=998244353;
int n,m,T,fac[N],h[N],g[N],f[N],inv[N],rev[N],tmp[N];

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

void DFT(int a[],int n,bool f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
        for (int p=i<<1,j=0; j<n; j+=p){
            int w=1;
            for (int k=0; k<i; k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
                a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    if (f) return;
    int inv=ksm(n,mod-2);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
}

void Inv(int a[],int b[],int l){
    if (l==1){ b[0]=ksm(a[0],mod-2); return; }
    Inv(a,b,l>>1); int n=1,L=0;
    for (; n<=l; n<<=1) L++;
    for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0; i<l; i++) tmp[i]=a[i],tmp[i+l]=0;
    DFT(tmp,n,1); DFT(b,n,1);
    for (int i=0; i<n; i++) tmp[i]=1ll*b[i]*(2-1ll*tmp[i]*b[i]%mod+mod)%mod;
    DFT(tmp,n,0);
    for (int i=0; i<l; i++) b[i]=tmp[i],b[i+l]=0;
}

int main(){
    n=100000;
    fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod;
    inv[n]=ksm(fac[n],mod-2);
    for (int i=n-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
    rep(i,1,n) h[i]=mod-inv[i]; h[0]=1;
    for (m=1; m<=n; m<<=1); Inv(h,g,m);
    for (int i=1; i<m; i++) f[i]=g[i];
    f[0]=(g[0]-1+mod)%mod; m<<=1;
    DFT(f,m,1); DFT(g,m,1);
    for (int i=0; i<m; i++) f[i]=1ll*f[i]*g[i]%mod;
    DFT(f,m,0); DFT(g,m,0);
    for (scanf("%d",&T); T--; ) scanf("%d",&n),printf("%lld
",1ll*f[n]*ksm(g[n],mod-2)%mod);
    return 0;
}