[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

3

Sample Output

87

HINT

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
",ans);
65    return 0;
66 }