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

Description

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

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

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

$S(i,j)$表示第二类斯特林数，递推公式为:
$S(i,j)=j imes S(i-1,j)+S(i-1,j-1),1leq jleq i-1$。
边界条件为：$S(i,i)=1(0leq i),S(i,0)=0(1leq 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}^nsum_{j=0}^i frac{1}{j!}sum_{k=0}^j(-1)^k C_j^k(j-k)^i imes 2^j imes(j!)$$

组合数展开得

$$sum_{i=0}^nsum_{j=0}^i 2^jsum_{k=0}^j(-1)^kfrac{j!}{k!(j-k)!}(j-k)^i$$

交换枚举顺序得

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

看到又有$k$又有$j-k$，这不是卷积吗？

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

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

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

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

那就直接NTT了，刚好给了一个NTT模数。（这莫非是提示？）

CODE：

#include<iostream>
#include<cstdio>
using namespace std;

#define mod 998244353
int n,bit=1,rev[400005],ans=0;
long long fac[100005],inv[100005],ivf[100005];
long long a[400005],b[400005];

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

void init(){
    fac[0]=ivf[0]=inv[0]=inv[1]=1;
    for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
    for(int i=2;i<=n;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;
    for(int i=1;i<=n;i++)ivf[i]=ivf[i-1]*inv[i]%mod;
    a[0]=b[0]=1;
    for(int i=1,f=-1;i<=n;i++,f=-f){
        if (i==1)a[i]=n+1;
        else a[i]=ivf[i]*(qpow(i,n+1)-1)%mod*inv[i-1]%mod;
        b[i]=f*ivf[i]%mod;
        if(a[i]<0)a[i]+=mod;
        if(b[i]<0)b[i]+=mod;
    }
}

void get_rev(){
    while(bit<=n+n)bit<<=1;
    for(int i=0;i<bit;i++)
    rev[i]=(rev[i>>1]>>1)|(i&1)*(bit>>1);
}

void NTT(long long a[],int dft){
    for(int i=0;i<bit;i++)
        if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=1;i<bit;i<<=1){
        long long W=qpow(3,(mod-1)/i/2);
        if(dft==-1)W=qpow(W,mod-2);
        for(int j=0;j<bit;j+=i<<1){
            long long w=1;
            for(int k=j;k<i+j;k++,w=w*W%mod){
                long long x=a[k];
                long long y=w*a[k+i]%mod;
                a[k]=(x+y)%mod,a[k+i]=(x+mod-y)%mod;
            }
        }
    }
    int inv=qpow(bit,mod-2);
    if(dft==-1)for(int i=0;i<bit;i++)a[i]=a[i]*inv%mod;
}

int main(){
    scanf("%d",&n);
    init();
    get_rev();
    NTT(a,1),NTT(b,1);
    for(int i=0;i<bit;i++)(a[i]*=b[i])%=mod;
    NTT(a,-1);
    for(int i=0;i<=n;i++){
        ans+=qpow(2,i)*fac[i]%mod*a[i]%mod;
        if(ans>=mod)ans-=mod;
    }
    printf("%d",ans);
}