「TJOI/HEOI2016」求和 - 指数型生成函数 + 多项式求逆 / 分治FFT

题目大意

    在2016年，佳媛姐姐刚刚学习了第二类斯特林数，非常开心。

    现在他想计算这样一个函数的值:

$$f(n)=\sum_{i=0}^n \sum_{j=0}^i S(i, j) \cdot 2^j \cdot j!$$

    $S(i,j)$表示第二类斯特林数，递推公式为:$S(i,j)=j\cdot S(i-1,j)+S(i-1,j-1),1\le j\le i-1$。

    边界条件为：$S(i, i) = 1(0 \leq i), \ S(i, 0) = 0(1 \leq i)$

    你能帮帮她吗?

---

题目分析

关于分治FFT直接看Candy?dalao的题解吧。
下面将多项式求逆的方法。

令$g(n)=\sum\limits_{i=1}^n {n \brace i}\cdot2^i\cdot i!$，那么答案就是其前缀和。
考虑$g(n)$的组合意义，代表将$n$个$1−n$的数分成若干集合，给每个集合黑白染色的方案数。

$$g(n)=\sum_{i=1}^n 2\cdot\binom ni\cdot g(n-i)$$

拆分组合数得到：

$$\frac{g(n)}{n!}=\sum_{i=1}^n \frac{2}{i!}\cdot\frac{g(n-i)}{(n-i)!}$$

又出现了，指数型生成函数！
令：

$$G(x)=\sum_{i=0}^{\infty}g(i)\frac{x^i}{i!}$$ $$F(x)=\sum_{i=1}^{\infty}2\frac{x^i}{i!}$$

那么有：

$$G(x)=F(x)G(x)+1$$

这里的$+1$是因为$G(x)$与$F(x)$卷积后第$0$次项（F从$1$开始卷积）没了，然而$G(x)$是有的，因此需要补$g(0)$回去，而根据组合意义（并非表达式）可以知道$g(0)=1$。

因此：

$$G(x)=\frac{1}{1-F(x)}$$

直接上多项式求逆即可。

---

代码

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

inline int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(!isdigit(x)) {if(x=='-')bj=-1;x=getchar();}
	while(isdigit(x)) {num=num*10+x-'0';x=getchar();}
	return num*bj;
}

const int maxn=262144+5;
const LL mod=998244353,g=3;

LL Quick_Pow(LL a,LL b) {
	LL sum=1;
	for(; b; b>>=1,a=a*a%mod)if(b&1)sum=sum*a%mod;
	return sum;
}

LL inv(LL x) {return Quick_Pow(x,mod-2);}

struct NumberTheoreticTransform {
	int n,rev[maxn];
	LL omega[maxn],iomega[maxn];
	void init(int n) {
		this->n=n;
		int x=Quick_Pow(g,(mod-1)/n);
		omega[0]=iomega[0]=1;
		for(int i=1; i<n; i++)omega[i]=omega[i-1]*x%mod,iomega[i]=inv(omega[i]);
		int k=log2(n);
		for(int i=0; i<n; i++) {
			int t=0;
			for(int j=0; j<k; j++)if(i&(1<<j))t|=1<<(k-j-1);
			rev[i]=t;
		}
	}
	void transform(LL *a,LL *omega) {
		for(int i=0; i<n; i++)if(i<rev[i])swap(a[i],a[rev[i]]);
		for(int len=2; len<=n; len*=2) {
			int mid=len>>1;
			for(LL *p=a; p!=a+n; p+=len)
				for(int i=0; i<mid; i++) {
					LL t=omega[n/len*i]*p[mid+i]%mod;
					p[mid+i]=(p[i]-t+mod)%mod;
					p[i]=(p[i]+t)%mod;
				}
		}
	}
	void dft(LL *a) {transform(a,omega);}
	void idft(LL *a) {
		transform(a,iomega);
		LL x=inv(n);
		for(int i=0; i<n; i++)a[i]=a[i]*x%mod;
	}
} ntt;

void polynomial_inverse(const LL* a,const int n,LL* b) {
	if(n==1) {b[0]=inv(a[0]);return;}
	polynomial_inverse(a,n>>1,b);
	int p=n<<1;
	static LL x[maxn];
	copy(a,a+n,x),fill(x+n,x+p,0);
	ntt.init(p),ntt.dft(x),ntt.dft(b);
	for(int i=0; i<p; i++)b[i]=b[i]*((2-x[i]*b[i]%mod+mod)%mod)%mod;
	ntt.idft(b),fill(b+n,b+p,0);
}

int n;
LL fac[maxn],invf[maxn],F[maxn],G[maxn];

int main() {
	n=Get_Int();
	int T=1;
	while(T<=n)T<<=1;
	fac[0]=1;
	for(int i=1; i<=T; i++)fac[i]=fac[i-1]*i%mod;
	invf[T]=inv(fac[T]);
	for(int i=T; i>=1; i--)invf[i-1]=invf[i]*i%mod;
	for(int i=1; i<=n; i++)F[i]=(mod-invf[i])*2%mod;
	F[0]=1;
	polynomial_inverse(F,T,G);
	LL ans=0;
	for(int i=n; i>=0; i--)ans=(ans+G[i]*fac[i]%mod)%mod;
	printf("%lld\n",ans);
	return 0;
}