[BZOJ4555][TJOI2016&HEOI2016]求和(分治FFT)
4555: [Tjoi2016&Heoi2016]求和
Time Limit: 40 Sec Memory Limit: 128 MB
Submit: 525 Solved: 418
[Submit][Status][Discuss]Description
在2016年,佳媛姐姐刚刚学习了第二类斯特林数,非常开心。
现在他想计算这样一个函数的值:S(i, j)表示第二类斯特林数,递推公式为:S(i, j) = j ∗ S(i − 1, j) + S(i − 1, j − 1), 1 <= j <= i − 1。边界条件为:S(i, i) = 1(0 <= i), S(i, 0) = 0(1 <= i)你能帮帮他吗?Input
输入只有一个正整数
Output
输出f(n)。由于结果会很大,输出f(n)对998244353(7 × 17 × 223 + 1)取模的结果即可。1 ≤ n ≤ 100000
Sample Input
3Sample Output
87HINT
Source
容易得到递推式,可以用CDQ分治+FFT
[l,mid]和[mid+1,r]卷起来怎么处理呢?平移数组变成[0,mid-l]和[mid-l+1,r-l+1]卷,次数界设为r-l+1即可。
代码用时:1h 比较顺利,没有低级错误。
实现比较简单,11348ms
1 #include<cstdio> 2 #include<algorithm> 3 #define rep(i,l,r) for (int i=l; i<=r; i++) 4 typedef long long ll; 5 using namespace std; 6 7 const int N=(1<<18)+100,P=998244353,g=3; 8 int n,rev[N]; 9 ll inv[N],fac[N],facinv[N],f[N],a[N],b[N]; 10 11 ll ksm(ll a,ll b){ 12 ll ans=1; 13 for (; b; b>>=1,a=a*a%P) 14 if (b & 1) ans=ans*a%P; 15 return ans; 16 } 17 18 void DFT(ll a[],int n,int f){ 19 rep(i,0,n-1) if (i<rev[i]) swap(a[i],a[rev[i]]); 20 for (int i=1; i<n; i<<=1){ 21 int wn=ksm(g,(f==1) ? (P-1)/(i<<1) : (P-1)-(P-1)/(i<<1)); 22 for (int p=i<<1,j=0; j<n; j+=p){ 23 int w=1; 24 for (int k=0; k<i; k++,w=1ll*w*wn%P){ 25 int x=a[j+k],y=1ll*w*a[i+j+k]%P; 26 a[j+k]=(x+y)%P; a[i+j+k]=(x-y+P)%P; 27 } 28 } 29 } 30 if (f==-1){ 31 int inv=ksm(n,P-2); 32 rep(i,0,n-1) a[i]=1ll*a[i]*inv%P; 33 } 34 } 35 36 void cdq(int l,int r){ 37 if (l==r) return; 38 int mid=(l+r)>>1,lim=r-l+1,n=1,L=0; 39 cdq(l,mid); 40 while (n<lim) n<<=1,L++; 41 rep(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1)); 42 rep(i,0,n-1) a[i]=b[i]=0; 43 rep(i,l,mid) a[i-l]=f[i]; 44 rep(i,0,r-l) b[i]=facinv[i]; 45 DFT(a,n,1); DFT(b,n,1); 46 rep(i,0,n-1) a[i]=a[i]*b[i]%P; 47 DFT(a,n,-1); 48 rep(i,mid+1,r) f[i]=(f[i]+2*a[i-l])%P; 49 cdq(mid+1,r); 50 } 51 52 int main(){ 53 freopen("bzoj4555.in","r",stdin); 54 freopen("bzoj4555.out","w",stdout); 55 scanf("%d",&n); inv[1]=1; fac[0]=facinv[0]=1; 56 rep(i,1,n){ 57 if (i!=1) inv[i]=(P-P/i)*inv[P%i]%P; 58 fac[i]=fac[i-1]*i%P; 59 facinv[i]=facinv[i-1]*inv[i]%P; 60 } 61 f[0]=1; cdq(0,n); ll ans=0; 62 rep(i,0,n) ans=(ans+f[i]*fac[i]%P)%P; 63 if (ans<0) ans+=P; 64 printf("%lld\n",ans); 65 return 0; 66 }