bzoj4555: 求和sum 快速傅立叶变换

题目大意

给定\(S(n,m)\)表示第二类斯特林数,定义函数\(f(n)\)

\[f(n) = \sum_{i=0}^n\sum_{j=0}^iS(i,j)*2^j*(j!) \]

给定正整数\(n,(n\leq 10^5)\),求\(f(n)\)

题解

我们都知道第二类斯特林数的递推公式为

\[S(i,j) = S(i-1,j-1) + j*S(i-1,j),(1 \leq j \leq i-1) \]

且有边界\(S(i,i) = 1(0 \leq i),S(i,0) = 0(1 \leq i)\)
第二类斯特林数\(S(i,j)\)的含义是把\(i\)个元素划分成\(j\)个无序的集合的方案
假设允许空集合的存在的话,方案即为\(m^n\)
我们应用容斥原理,枚举至少有多少空集合空集合,那么有

\[S(n,m) = \frac{1}{m!}\sum_{k=1}^{m}C_m^k(m-k)^n(-1)^k \]

设\(g(n) = \sum_{i=0}^nS(n,i)2^i(i!)\)
那么我们将\(S(n,m)\)代入\(g(n)\)化简得

\[g(n) = \sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!} \]

那么将\(g(n)\)带入答案表达式中,有

\[ans = \sum_{n=0}^x\sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!} \]

这时我们发现每次最外层的\(n -> (n+1)\)时,都相当于在内部的\(\frac{(m-k)^n}{(m-k)!}\)一项上又加上了一个\(\frac{(m-k)^{n+1}}{(m-k)!}\)
所以我们把这一项做等比数列求和.
设\(g(x) = \frac{x^{n+1} - x}{(x-1)(x!)}\)
那么上式变成了

\[ans = \sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k) \]

于是我们在\(\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k)\)进行FFT计算卷积
这样就只剩下了一个sigma式,for循环一边即可.
复杂度\(O(nlogn)\)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
template<typename T>inline void read(T &x){
	x=0;char ch;bool flag = false;
	while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
	while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
const int maxn = 600010;
const int mod = 998244353;
const int pri_rt = 3;
int w[maxn];
inline int qpow(int x,int p){
	int ret = 1;
	for(;p;p>>=1,x=1LL*x*x%mod) if(p&1) ret=1LL*ret*x % mod;
	return ret;
}
inline void FNT(int *x,int n,int p){
	for(int i=0,t=0;i<n;++i){
		if(i > t) swap(x[i],x[t]);
		for(int j=n>>1;(t^=j)<j;j>>=1);
	}
	for(int m=2;m<=n;m<<=1){
		int k = m>>1;
		int wn = qpow(pri_rt,p == 1 ? (mod-1)/m : (mod-1) - (mod-1)/m);
		for(int i=1;i<k;++i) w[i] = 1LL*w[i-1]*wn % mod;
		w[0] = 1;
		for(int i=0;i<n;i+=m){
			for(int j=0;j<k;++j){
				int u = 1LL*x[i+j+k]*w[j] % mod;
				x[i+j+k] = x[i+j] - u;
				if(x[i+j+k] < 0) x[i+j+k] += mod;
				x[i+j] += u;
				if(x[i+j] >= mod) x[i+j] -= mod;
			}
		}
	}
	if(p == -1){
		int inv = qpow(n,mod-2);
		for(int i=0;i<n;++i) x[i] = 1LL*x[i]*inv % mod;
	}
}
int fac[maxn],inv[maxn];
inline void init(int n){
	fac[0] = 1;
	for(int i=1;i<=n;++i) fac[i] = 1LL*fac[i-1]*i % mod;
	inv[n] = qpow(fac[n],mod-2);
	for(int i = n-1;i>=0;--i) inv[i] = 1LL*inv[i+1]*(i+1) % mod;
}
int A[maxn],B[maxn];
int main(){
	int n;read(n);
	int len;for(len=1;len <= (n+1);len<<=1);len<<=1;
	init(n);
	for(int i=0;i<=n;++i){
		if(i&1) A[i] = -inv[i] + mod;
		else A[i] = inv[i];
	}
	for(int i=2;i<=n;++i){
		B[i] = qpow(i,n+1) - i + mod;
		if(B[i] < 0) B[i] += mod;
		B[i] = (1LL*B[i]*qpow(i-1,mod-2)%mod*inv[i]) % mod;
	}B[1] = n;
	FNT(A,len,1);FNT(B,len,1);
	for(int i=0;i<len;++i) A[i] = 1LL*A[i]*B[i] % mod;
	FNT(A,len,-1);
	int ans = 1;
	for(int i=1,f2=2;i<=n;++i){
		ans = (ans + 1LL*A[i]*f2%mod*fac[i]) % mod;
		f2 = (f2<<1) % mod;
	}printf("%d\n",ans);
	getchar();getchar();
	return 0;
}