【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 }