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

题目大意

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

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

给定正整数(n,(nleq 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}^mfrac{(-1)^k}{k!}frac{(m-k)^n}{(m-k)!} ]

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

[ans = sum_{n=0}^xsum_{m=0}^n2^m(m!)sum_{k=0}^mfrac{(-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}^mfrac{(-1)^k}{k!}g(m-k) ]

于是我们在(sum_{k=0}^mfrac{(-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
",ans);
	getchar();getchar();
	return 0;
}
原文地址:https://www.cnblogs.com/Skyminer/p/6402254.html