P5162 WD與積木

題目背景

WD整日沉浸在積木中，無法自拔……

題目描述

WD想買\(n\)塊積木，商場中每塊積木的高度都是\(1\)，俯檢視為正方形（邊長不一定相同）。由於一些特殊原因，商家會給每個積木隨機一個大小並標號，發給WD。

接下來WD會把相同大小的積木放在一層，並把所有層從大到小堆起來。WD希望知道所有不同的堆法中層數的期望。兩種堆法不同當且僅當某個積木在兩種堆法中處於不同的層中，由於WD只關心積木的相對大小，因此所有堆法等概率出現，而不是隨機的大小等概率（可以看樣例理解）。輸出結果\(\bmod 998244353\)即可。

（如果還是不能夠理解題意，請看樣例）

輸入輸出格式

輸入格式：

第一行一個數\(T\)，表示詢問個數。

接下來\(T\)行每行一個數\(n\)，表示WD希望使用\(n\)塊積木。

輸出格式：

共\(T\)行，每行一個數表示答案\(\bmod 998244353\)。

說明

subtask1(21pts): \(1≤T≤1,000, 1≤n≤1,000\)

subtask2(37pts):\(~1\le T\le 10,~1\le n\le 100,000\)

subtask3(42pts):\(~1\le T\le 100,000,~1\le n\le 100,000\)

當時拿斯特林數推了一會兒沒推出來就走了，原來只是部分分啊。

姿勢水平屑了，今天稍微給自己普及了一些指數型生成函式。

樸素DP

\(g_n\)代表\(n\)個東西的堆法，就是有標號球和盒子不準空的方案數，就是有序貝爾數，遞推的時候列舉第一堆大小。
\[ g_n=\sum_{i=1}^n \binom{n}{i}g_{n-i} \]
\(f_i\)代表\(n\)個東西所有堆法的總層數。
\[ f_n=g_n+\sum_{i=1}^n\binom{n}{i}f_{n-i} \]
從意義上看，\(f_0=0,g_0=1\)，不想從意義上看就稍微推一下。

實質還是列舉第一堆的大小，\(g_n\)算的是列舉的這一堆的貢獻和，剩下的是其他堆的貢獻。

然後構造指數生成函式
\[ F(x)=\sum_{i=0}^\infty\frac{f_i}{i!}x^i,G(x)=\sum_{i=0}^\infty\frac{g_i}{i!}x^i,H(x)=\sum_{i=0}^\infty\frac{1}{i!}x^i \]

\(1\)實際上是一個多項式單位元，像數論卷積定義的\(\epsilon(n)=[n=1]\)那樣。

然後＋1還是為了\(g_0=1\)考慮的，因為本身這個東西已經封閉了，你得給它開啟一個開口啊（天吶我在扯什麼

後面-1減的實際上還是\(g_0=1\)

化簡一下式子
\[ F=G(G-1) \]
求逆求出\(G\)就行了

Code:

#include <cstdio>
#include <algorithm>
const int mod=998244353,Gi=332748118;
const int M=(1<<18)+10;
const int N=1e5;
#define mul(a,b) (1ll*(a)*(b)%mod)
#define add(a,b) ((a+b)%mod)
int qp(int d,int k){int f=1;while(k){if(k&1)f=mul(f,d);d=mul(d,d),k>>=1;}return f;}
int fac[M],inv[M],ans[M],H[M],A[M],B[M],D[N],turn[M],len=1;
void NTT(int *a,int typ,int len)
{
    int L=-1;for(int i=1;i<len;i<<=1) ++L;
    for(int i=1;i<len;i++)
    {
        turn[i]=turn[i>>1]>>1|(i&1)<<L;
        if(i<turn[i]) std::swap(a[i],a[turn[i]]);
    }
    for(int le=1;le<len;le<<=1)
    {
        int wn=qp(typ?3:Gi,(mod-1)/(le<<1));
        for(int p=0;p<len;p+=le<<1)
        {
            int w=1;
            for(int i=p;i<p+le;i++,w=mul(w,wn))
            {
                int tx=a[i],ty=mul(w,a[i+le]);
                a[i]=add(tx,ty);
                a[i+le]=add(tx,mod-ty);
            }
        }
    }
    if(!typ)
    {
        int inv=qp(len,mod-2);
        for(int i=0;i<len;i++) a[i]=mul(a[i],inv);
    }
}
void polyinv(int *a,int *b,int len)
{
    if(len==1){b[0]=qp(a[0],mod-2);return;}
    polyinv(a,b,len>>1);
    for(int i=0;i<len;i++) A[i]=b[i],B[i]=a[i],A[i+len]=B[i+len]=0;
    NTT(A,1,len<<1),NTT(B,1,len<<1);
    for(int i=0;i<len<<1;i++) A[i]=mul(A[i],add(2,mod-mul(A[i],B[i])));
    NTT(A,0,len<<1);
    for(int i=0;i<len;i++) b[i]=A[i];
}
int main()
{
    fac[0]=1;for(int i=1;i<=N;i++) fac[i]=mul(fac[i-1],i);
    inv[N]=qp(fac[N],mod-2);
    for(int i=N-1;~i;i--) inv[i]=mul(inv[i+1],i+1);
    for(int i=0;i<=N;i++) H[i]=mod-inv[i];H[0]=1;
    while(len<=N) len<<=1;
    polyinv(H,ans,len);
    for(int i=0;i<len;i++) D[i]=mul(ans[i],fac[i]),H[i]=ans[i];
    --ans[0];
    NTT(ans,1,len<<1),NTT(H,1,len<<1);
    for(int i=0;i<len<<1;i++) ans[i]=mul(ans[i],H[i]);
    NTT(ans,0,len<<1);
    for(int i=0;i<=N;i++) ans[i]=mul(qp(D[i],mod-2),mul(ans[i],fac[i]));
    int T,n;scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&n);
        printf("%d\n",ans[n]);
    }
    return 0;
}

2018.12.31