【BZOJ4555】【TJOI2016】【HEOI2016】求和 (第二类斯特林数+NTT卷积)

Description

在2016年,佳媛姐姐刚刚学习了第二类斯特林数,非常开心。

现在他想计算这样一个函数的值:

$$f(n)=\sum_{i=0}^n\sum_{j=0}^i S(i,j)\times 2^j\times(j!)$$

$S(i,j)$表示第二类斯特林数,递推公式为:
$S(i,j)=j\times S(i-1,j)+S(i-1,j-1),1\leq j\leq i-1$。
边界条件为:$S(i,i)=1(0\leq i),S(i,0)=0(1\leq i)$
你能帮帮他吗?

Input

输入只有一个正整数$n$

Output

输出$f(n)$。

由于结果会很大,输出$f(n)$对$998244353(7×17×223+1)$取模的结果即可。

题解:

递推式给你就是玩你的,一点关系都没有!

第二类斯特林数通项公式:

$$S(i,j)=\frac{1}{j!}\sum_{k=0}^j(-1)^k C_j^k(j-k)^i$$

代入原式得

$$\sum_{i=0}^n\sum_{j=0}^i \frac{1}{j!}\sum_{k=0}^j(-1)^k C_j^k(j-k)^i\times 2^j\times(j!)$$

组合数展开得

$$\sum_{i=0}^n\sum_{j=0}^i 2^j\sum_{k=0}^j(-1)^k\frac{j!}{k!(j-k)!}(j-k)^i$$

交换枚举顺序得

$$\sum_{j=0}^n(j!)2^j\sum_{k=0}^n{(-1)^k\over k!}{1\over (j-k)!}\sum_{i=0}^n(j-k)^i$$

看到又有$k$又有$j-k$,这不是卷积吗?

令 $a(x)=\sum_{x=0}^n{(-1)^x\over x!}$

    $b(x)={1\over (x)!}\sum_{i=0}^n(x)^i$

则 $f(x)=(x!)2^x\sum_{i=0}^n a(x)b(x-i)$

不过循环边界到n不会越界吗?这对答案没有影响,因为斯特林数和组合数都为0。

那就直接NTT了,刚好给了一个NTT模数。(这莫非是提示?)

CODE:

 1 #include<iostream>
 2 #include<cstdio>
 3 using namespace std;
 4 
 5 #define mod 998244353
 6 int n,bit=1,rev[400005],ans=0;
 7 long long fac[100005],inv[100005],ivf[100005];
 8 long long a[400005],b[400005];
 9 
10 int qpow(int x,int y){
11     int ans=1;
12     while(y){
13         if(y&1)ans=1LL*ans*x%mod;
14         y>>=1,x=1LL*x*x%mod;
15     }
16     return ans;
17 }
18 
19 void init(){
20     fac[0]=ivf[0]=inv[0]=inv[1]=1;
21     for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
22     for(int i=2;i<=n;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
23     for(int i=1;i<=n;i++)ivf[i]=ivf[i-1]*inv[i]%mod;
24     a[0]=b[0]=1;
25     for(int i=1,f=-1;i<=n;i++,f=-f){
26         if (i==1)a[i]=n+1;
27         else a[i]=ivf[i]*(qpow(i,n+1)-1)%mod*inv[i-1]%mod;
28         b[i]=f*ivf[i]%mod;
29         if(a[i]<0)a[i]+=mod;
30         if(b[i]<0)b[i]+=mod;
31     }
32 }
33 
34 void get_rev(){
35     while(bit<=n+n)bit<<=1;
36     for(int i=0;i<bit;i++)
37     rev[i]=(rev[i>>1]>>1)|(i&1)*(bit>>1);
38 }
39 
40 void NTT(long long a[],int dft){
41     for(int i=0;i<bit;i++)
42         if(i<rev[i])swap(a[i],a[rev[i]]);
43     for(int i=1;i<bit;i<<=1){
44         long long W=qpow(3,(mod-1)/i/2);
45         if(dft==-1)W=qpow(W,mod-2);
46         for(int j=0;j<bit;j+=i<<1){
47             long long w=1;
48             for(int k=j;k<i+j;k++,w=w*W%mod){
49                 long long x=a[k];
50                 long long y=w*a[k+i]%mod;
51                 a[k]=(x+y)%mod,a[k+i]=(x+mod-y)%mod;
52             }
53         }
54     }
55     int inv=qpow(bit,mod-2);
56     if(dft==-1)for(int i=0;i<bit;i++)a[i]=a[i]*inv%mod;
57 }
58 
59 int main(){
60     scanf("%d",&n);
61     init();
62     get_rev();
63     NTT(a,1),NTT(b,1);
64     for(int i=0;i<bit;i++)(a[i]*=b[i])%=mod;
65     NTT(a,-1);
66     for(int i=0;i<=n;i++){
67         ans+=qpow(2,i)*fac[i]%mod*a[i]%mod;
68         if(ans>=mod)ans-=mod;
69     }
70     printf("%d",ans);
71 }
作者:ezoiLZH
本文版权归作者和博客园所有,欢迎转载,只要写明原文链接即可(^_^)。
posted @ 2018-08-05 16:24  ezoiLZH  阅读(377)  评论(0编辑  收藏  举报