1. 程式人生 > >NTT+分治FFT--P4091 [HEOI2016/TJOI2016]求和

NTT+分治FFT--P4091 [HEOI2016/TJOI2016]求和

傳送門

這道題很妙啊

首先看題目中的式子,令新的 f ( n ) = i =

0 n S ( n , i ) × 2
i
× ( i ! ) f(n)=\sum_{i=0}^nS(n,i)\times 2^i\times (i!) ,如果能快速求出這個式子的值,那麼 a
n s = i = 0 n f ( i ) ans=\sum_{i=0}^n f(i)

首先設 g ( n ) = i = 0 n S ( n , i ) × i ! g(n)=\sum_{i=0}^nS(n,i)\times i! ,考慮組合意義,第二類 s t r i l i n g striling 數是 n n 個不同小球放入 m m 個相同盒子,那麼 g ( n ) g(n) 就代表把 n n 個小球放入任意多 ( n ) (\le n) 的不同盒子的方案數。
也就可以轉化成 g ( n ) = i = 0 n C ( n , i ) × g ( i ) g(n)=\sum_{i=0}^nC(n,i)\times g(i)

考慮原式: f ( n ) = i = 0 n C ( n , i ) × g ( i ) × 2 i f(n)=\sum_{i=0}^n C(n,i)\times g(i)\times 2^i
其實就相當於給每個盒子黑白染色

也可以轉化成: f ( n ) = i = 0 n 1 C ( n , i ) × f ( i ) f(n)=\sum_{i=0}^{n-1} C(n,i)\times f(i)

把組合數展開,發現式子可以化簡為: f ( n ) n ! = i = 0 n 1 2 ( n i ) ! × f ( i ) i ! \frac{f(n)}{n!}=\sum_{i=0}^{n-1}\frac{2}{(n-i)!}\times \frac{f(i)}{i!}

到這步就可以看到右邊的兩部分下標和為 i i ,就可以分治 F F T FFT 來求出 f f ,最後求答案的時候只需要 × i ! \times i!
因為題目要取模,模數剛好是 998244353 998244353 ,原根為 3 3 ,直接 N T T NTT 即可

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define maxn 400005
#define LL long long
using namespace std;
int n,limit=1,lim,rev[maxn];
LL a[maxn],b[maxn],f[maxn],fac[maxn],inv[maxn],ans;
const int mod=998244353,g=3;

inline LL qpow(LL x,int k){
	LL ret=1;
	while(k){
		if(k&1) (ret*=x)%=mod;
		(x*=x)%=mod; k>>=1;
	} return ret%mod;
}

inline void NTT(LL *F,int type){
	for(int i=0;i<limit;i++)
		if(i<rev[i]) swap(F[i],F[rev[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		LL Wn=qpow(g,type==1?(mod-1)/(mid<<1):(mod-1-(mod-1)/(mid<<1)));
		for(int r=mid<<1,j=0;j<limit;j+=r){
			LL w=1;
			for(int k=0;k<mid;k++,(w*=Wn)%=mod){
				LL x=F[j+k],y=F[j+mid+k]*w%mod;
				F[j+k]=(x+y)%mod,F[j+mid+k]=(x-y+mod)%mod;
			}
		}
	}
	if(type==-1){
		LL INV=qpow(limit,mod-2);
		for(int i=0;i<limit;i++) F[i]=(F[i]*INV)%mod;
	}
}

inline void work(){
	NTT(a,1); NTT(b,1);
	for(int i=0;i<limit;i++) a[i]=a[i]*b[i]%mod;
	NTT(a,-1);
}

inline void cdqFFT(int l,int r){
	if(l==r) return;
	int mid=(l+r)>>1,len=r-l-1;
	cdqFFT(l,mid);
	limit=1,lim=0;
	while(limit<=len) limit<<=1,++lim;
	for(int i=0;i<limit;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1)),a[i]=b[i]=0;//清零 
	for(int i=l;i<=mid;i++) a[i-l]=f[i];
	for(int i=1;i<=r-l;i++) b[i-1]=2LL*inv[i]%mod;
	work();
	for(int i=mid+1;i<=r;i++) f[i]=(f[i]+a[i-l-1]%mod)%mod;
	cdqFFT(mid+1,r);
}

int main(){
	scanf("%d",&n); fac[1]=fac[0]=1; f[0]=1;
	for(int i=2;i<=n;i++) fac[i]=1LL*fac[i-1]*i%mod;
	inv[n]=qpow(fac[n],mod-2);
	for(int i=n;i>1;i--) inv[i-1]=1LL*inv[i]*i%mod;
	cdqFFT(0,n);
	for(int i=0;i<=n;i++) (ans+=f[i]*fac[i]%mod)%=mod;
	printf("%lld\n",ans);
	return 0;
}