【生成函数+多项式求逆】LGP5162 WD与积木

【题目】
原题地址
有$n$块积木，给每块积木随机一个大小并标号，然后将相同大小的积木放在一层，再从大到小堆起来。我们只关心积木的相对大小，因此所有堆法等概率出现，求期望层数。$T,n\le 10^5$

【解题思路】
从期望的定义入手，我们先考虑一个朴素的$\text{DP}$：设$f_i$表示有$i$块积木时产生层数和，$g_i$表示$i$块积木不同堆法的数量，答案就是$\frac{f_n}{g_n}$，其中$f_0=0,g_0=1$。我们枚举第一层放几块积木可以得到转移：
$$g_n=\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)g_{n-i}$$
$$f_n=\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)(f_{n-i}+g_{n-i})=g_n+\sum _{i=1}^n\left(\genfrac{}{}{0ex}{}{n}{i}\right)f_{n-i}$$
这样做是$O(n^2)$的，这里好像还有一个标号顺序的问题我们好像没有考虑，但由于方案和层数是一一对应的，所以相当于同时约掉了。

我们考虑一个更为正确的写法：不妨设
$$F(x)=\sum _i\frac{f_i}{i!}{x}^i,G(x)=\sum _i\frac{g_i}{i!}{x}^i,H(x)=\sum _i\frac{x^i}{i!}$$
由于没有组合数的计算，我们这里写出的生成函数要考虑顺序（实际上和前几天的生成函数题类似，答案的形式都是$\frac{n!}{k_1!k_2!\dots }$）

我们可以得到
$$2F(x)=G(x)+F(x)H(x)-1,2G(x)=G(x)H(x)+1$$
$$F(x)=G(x)(G(x)-1),G(x)=\frac{1}{2-H(x)}$$

于是多项式求逆可以做到$O(n\log n)$

【参考代码】

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N=1e5+5,M=264233;
const int mod=998244353,G=3;
int fac[M],ifac[M];

int read()
{
	int ret=0;char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) ret=ret*10+(c^48),c=getchar();
	return ret;
}
int qpow(int x,int y)
{
	int res=1;
	for(;y;y>>=1,x=(ll)x*x%mod) if(y&1) res=(ll)res*x%mod;
	return res;
}
int upm(int x){return x>=mod?x-mod:(x<0?x+mod:x);}
void up(int &x,int y){x=upm(x+y);}

namespace NTT
{
	int m,L;
	int f[M],g[M],h[M],t[M],t2[M];
	int rev[M];
	void reget(int n)
	{
		for(L=0,m=1;m<2*n;++L,m<<=1);
		for(int i=0;i<m;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	}
	void ntt(int *a,int n,int op)
	{
		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=qpow(G,(mod-1)/(i<<1));
			if(!~op) wn=qpow(wn,mod-2);
			for(int j=0;j<n;j+=(i<<1))
			{
				int w=1;
				for(int k=0;k<i;++k,w=(ll)w*wn%mod)
				{
					int x=a[j+k],y=(ll)w*a[i+j+k]%mod;
					a[j+k]=upm(x+y);a[i+j+k]=upm(x-y);
				}
			}
		}
		int inv=qpow(n,mod-2);
		if(!~op) for(int i=0;i<n;++i) a[i]=(ll)a[i]*inv%mod;
	}
	void polyinv(int *a,int *b,int deg)
	{
		if(deg==1){b[0]=qpow(a[0],mod-2);return;}
		polyinv(a,b,(deg+1)>>1);
		reget(deg);
		for(int i=0;i<deg;++i) t[i]=a[i];
		for(int i=deg;i<m;++i) t[i]=0;
		ntt(b,m,1);ntt(t,m,1);
		for(int i=0;i<m;++i) b[i]=(ll)(upm(2-(ll)b[i]*t[i]%mod))*b[i]%mod;
		ntt(b,m,-1);
		for(int i=deg;i<m;++i) b[i]=0;
	}
	void mul(int *a,int *b,int *c,int deg)
	{
		reget(deg);
		for(int i=0;i<deg;++i) t[i]=a[i],c[i]=b[i];
		for(int i=deg;i<m;++i) t[i]=c[i]=0;
		ntt(t,m,1);ntt(c,m,1);
		for(int i=0;i<m;++i) c[i]=(ll)t[i]*c[i]%mod;
		ntt(c,m,-1);
		for(int i=deg;i<m;++i) c[i]=0;
	}
}
using namespace NTT;

void init()
{
	fac[0]=1;for(int i=1;i<N;++i)fac[i]=(ll)fac[i-1]*i%mod; 
	for(int i=0;i<N;++i)ifac[i]=qpow(fac[i],mod-2);
	for(int i=0;i<N;++i) h[i]=mod-ifac[i]; h[0]=upm(h[0]+2);
	polyinv(h,g,N);
	for(int i=0;i<N;++i) f[i]=g[i]; f[0]=upm(g[0]-1);
	mul(f,g,f,N);
}
void solve()
{
	int T=read();
	while(T--)
	{
		int x=read();
		printf("%d\n",(ll)f[x]*qpow(g[x],mod-2)%mod);
	}
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("LGP5162.in","r",stdin);
	freopen("LGP5162.out","w",stdout);
#endif
	init();solve();
	return 0;
}

【总结】
事实证明我的$\text{DP}$学得并不好，开始的式子想了很久。
然而后面的生成函数一下子就搞出来了。。。