不要以為用上Stirling數就一定離正解更近，FFT都是從DP式本身出發的。

設f[i]為i個積木的所有方案的層數總和，g[i]為i個積木的方案數，則答案為$\frac{f[i]}{g[i]}$

轉移列舉第一層是哪些積木：$$f_n=g_n+\sum_{i=1}^{n}\binom{n}{i}f_{n-i},f_0=0$$$$g_n=\sum_{i=1}^{n}\binom{n}{i}g_{n-i},g_0=1$$

轉化成卷積形式：$$\frac{f_n}{n!}=\frac{g_n}{n!}+\sum_{i=1}^{n}\frac{1}{i!}\times \frac{f_{n-i}}{i!}$$$$\frac{g_n}{n!}=\sum_{i=1}^{n}\frac{1}{i!}\times \frac{g_{n-i}}{(n-i)!}$$

構造生成函式：$F(x)=\sum\frac{f_i}{i!}x^i$,$G(x)=\sum\frac{g_i}{i!}x^i$,$H(x)=\sum\frac{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)}$。

多項式求逆即可。

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define
 
 rep(i,l,r) for (int i=(l); i<=(r); i++)
 4 using namespace std;
 5 
 6 const int N=500010,mod=998244353;
 7 int n,m,T,fac[N],h[N],g[N],f[N],inv[N],rev[N],tmp[N];
 8 
 9 int ksm(int a,int b){
10     int res=1;
11     for (; b; a=1ll*a*a%mod,b>>=1)
12         if (b & 1) res=1ll*res*a%mod;
 
13     return res;
14 }
15 
16 void DFT(int a[],int n,bool f){
17     for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
18     for (int i=1; i<n; i<<=1){
19         int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
20         for (int p=i<<1,j=0; j<n; j+=p){
21             int w=1;
22             for (int k=0; k<i; k++,w=1ll*w*wn%mod){
23                 int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
24                 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
25             }
26         }
27     }
28     if (f) return;
29     int inv=ksm(n,mod-2);
30     for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
31 }
32 
33 void Inv(int a[],int b[],int l){
34     if (l==1){ b[0]=ksm(a[0],mod-2); return; }
35     Inv(a,b,l>>1); int n=1,L=0;
36     for (; n<=l; n<<=1) L++;
37     for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
38     for (int i=0; i<l; i++) tmp[i]=a[i],tmp[i+l]=0;
39     DFT(tmp,n,1); DFT(b,n,1);
40     for (int i=0; i<n; i++) tmp[i]=1ll*b[i]*(2-1ll*tmp[i]*b[i]%mod+mod)%mod;
41     DFT(tmp,n,0);
42     for (int i=0; i<l; i++) b[i]=tmp[i],b[i+l]=0;
43 }
44 
45 int main(){
46     n=100000;
47     fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod;
48     inv[n]=ksm(fac[n],mod-2);
49     for (int i=n-1; ~i; i--) inv[i]=1ll*inv[i+1]*(i+1)%mod;
50     rep(i,1,n) h[i]=mod-inv[i]; h[0]=1;
51     for (m=1; m<=n; m<<=1); Inv(h,g,m);
52     for (int i=1; i<m; i++) f[i]=g[i];
53     f[0]=(g[0]-1+mod)%mod; m<<=1;
54     DFT(f,m,1); DFT(g,m,1);
55     for (int i=0; i<m; i++) f[i]=1ll*f[i]*g[i]%mod;
56     DFT(f,m,0); DFT(g,m,0);
57     for (scanf("%d",&T); T--; ) scanf("%d",&n),printf("%lld\n",1ll*f[n]*ksm(g[n],mod-2)%mod);
58     return 0;
59 }