Luogu5162 WD与积木（生成函数+多项式求逆）

Posted 2021-01-31 gloid

　　显然的做法是求出斯特林数，但没有什么优化空间。

　　考虑一种暴力dp，即设f[i]为i块积木的所有方案层数之和，g[i]为i块积木的方案数。转移时枚举第一层是哪些积木，于是有f[i]=g[i]+ΣC(i,j)·f[i-j]，g[i]=ΣC(i,j)·g[i-j] (j=1~i)。

　　考虑优化 。我们发现这个转移非常像卷积。写成卷积形式，有f[i]=g[i]+Σi!·Σf[i-j]/j!/(i-j)!，g[i]=i!·Σg[i-j]/j!/(i-j)!。直接分治NTT即可。

　　诶是不是强行多了个log？考虑构造出生成函数。g比较简单，将该式写成g[i]/i!=Σg[i-j]/j!/(i-j)!，设g[i]/i!生成函数为G(x)，inv(i!)生成函数为H(x)，则有G(x)=G(x)·H(x)-G(x)+1，也即G(x)=1/(2-H(x))。f同样写成f[i]/i!=g[i]/i!+Σf[i-j]/j!/(i-j)!，则有F(x)=G(x)+F(x)·H(x)-F(x)-1，也即F(x)=G(x)·(G(x)-1)。多项式求逆即可。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define P 998244353
#define N 100010
char getc(){char c=getchar();while ((c<‘A‘||c>‘Z‘)&&(c<‘a‘||c>‘z‘)&&(c<‘0‘||c>‘9‘)) c=getchar();return c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<‘0‘||c>‘9‘) {if (c==‘-‘) f=-1;c=getchar();}
    while (c>=‘0‘&&c<=‘9‘) x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int T,n,t,r[1<<19],f[1<<19],g[1<<19],h[1<<19],tmp[1<<19];
int ksm(int a,int k)
{
    int s=1;
    for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P;
    return s;
}
int inv(int a){return ksm(a,P-2);}
void DFT(int *a,int n,int g)
{
    for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (int i=2;i<=n;i<<=1)
    {
        int wn=ksm(g,(P-1)/i);
        for (int j=0;j<n;j+=i)
        {
            int w=1;
            for (int k=j;k<j+(i>>1);k++,w=1ll*w*wn%P)
            {
                int x=a[k],y=1ll*w*a[k+(i>>1)]%P;
                a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P;
            }
        }
    }
}
void mul(int *a,int *b,int n,int op)
{
    for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1);
    DFT(a,n,3),DFT(b,n,3);
    if (op==0) for (int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%P;
    else for (int i=0;i<n;i++) a[i]=1ll*a[i]*(P+2-1ll*a[i]*b[i]%P)%P;
    DFT(a,n,inv(3));
    int t=inv(n);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*t%P;
}
void getinv(int n)
{
    g[0]=1;
    for (int i=2;i<=n;i<<=1)
    {
        for (int j=i;j<(i<<1);j++) tmp[j]=0;
        for (int j=0;j<i;j++) tmp[j]=h[j];
        mul(g,tmp,i<<1,1);
        for (int j=i;j<(i<<1);j++) g[j]=0;
    }
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("a.in","r",stdin);
    freopen("a.out","w",stdout);
    const char LL[]="%I64d
";
#else
    const char LL[]="%lld
";
#endif
    T=read();int t=1;while (t<=(N<<1)) t<<=1;
    h[0]=h[1]=1;for (int i=2;i<=N;i++) h[i]=P-1ll*(P/i)*h[P%i]%P;
    for (int i=2;i<=N;i++) h[i]=1ll*h[i]*h[i-1]%P;
    for (int i=0;i<=N;i++) h[i]=(P-h[i])%P;
    h[0]=(h[0]+2)%P;
    getinv(t);
    for (int i=N+1;i<t;i++) g[i]=0;
    memcpy(f,g,sizeof(f));f[0]--;if (f[0]<0) f[0]+=P;
    mul(f,g,t,0);DFT(g,t,inv(3));int u=inv(t);for (int i=0;i<t;i++) g[i]=1ll*g[i]*u%P;
    while (T--)
    {
        n=read();
        printf("%d
",1ll*f[n]*inv(g[n])%P);
    }
    return 0;
}