[HEOI2016/TJOI2016][bzoj4555] 求和 [斯特林数+NTT]

Posted 2020-11-12 Orion545

题面

传送门

思路

首先，我们发现这个式子中大部分的项都和\(j\)有关（尤其是后面的\(2^j\ast j!\)），所以我们更换一下枚举方式，把这道题的枚举方式变成先\(j\)再\(i\)

\(f(n)=\sum_{j=0}^n2^j\ast j!\sum_{i=0}^nS_i^j\)

第二类斯特林数有一个基于组合意义的公式：

\(S_i^j=\frac1{j!}\sum_{k=0}^j(-1)^kC_j^k(j-k)^i=\sum_{k=0}^j\frac{(-1)^k(j-k)^i}{k!(j-k)!}\)

把这个公式代回原式中，得到：

\(f(n)=\sum_{j=0}^n2^j\ast j!\sum_{i=0}^n\sum_{k=0}^j\frac{(-1)^k(j-k)^i}{k!(j-k)!}\)

再次更换一下枚举方式，变成：

\(f(n)=\sum_{j=0}^n2^j\ast j!\sum_{k=0}^j\frac{(-1)^k}{k!}\sum_{i=0}^n\frac{(j-k)^i}{(j-k)!}\)

\(f(n)=\sum_{j=0}^n2^j\ast j!\sum_{k=0}^j\frac{(-1)^k}{k!}\ast\frac{\sum_{i=0}^n(j-k)^i}{(j-k)!}\)

此时，设两个函数\(a\)和\(b\)，令：

\(a(i)=\frac{(-1)^i}{i!}\)

\(b(i)=\frac{\sum_{j=0}^ni^j}{i!}=\frac{i^{n+1}-1}{(i-1)i!}\)

那么，

\(f(n)=\sum_{j=0}^n 2^j\ast j!\ast(a\ast b)(j)\)

其中(a\ast b)(j)表示\(a\)和\(b\)的\(0-j\)项的卷积

模数为\(998244353\)，用\(NTT\)做一遍卷积即可，时间效率为\(O(nlog_2n)\)

注意事项

\(b(0)=1,b(1)=n+1\)

这两个要提前保存一下，因为用公式推的话会div 0

还有一个奇怪的问题我没有解决，具体看代码最后面吧

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
inline int read(){
    int re=0,flag=1;char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') flag=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9') re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
    return re*flag;
}
#define ll long long
ll MOD=998244353,g=3,inv[400010],f[400010],finv[400010];
int qpow(ll a,ll b){//快速幂
    ll re=1;
    while(b){
        if(b&1) re=re*a%MOD;
        a=a*a%MOD;b>>=1;
    }
    return re;
}
ll n,A[400010],B[400010],C[400010],r[400010],limit,cnt;
void ntt(ll *a,ll type){
    int i,j,k,mid;ll y,w,wn;
    for(i=0;i<limit;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    for(mid=1;mid<limit;mid<<=1){
        wn=qpow((type==1)?g:inv[g],(MOD-1)/(mid<<1));
        for(j=0;j<limit;j+=(mid<<1)){
            w=1;
            for(k=0;k<mid;k++,w=w*wn%MOD){
                y=a[j+k+mid]*w%MOD;
                a[j+k+mid]=(a[j+k]-y+MOD)%MOD;
                a[j+k]=(a[j+k]+y)%MOD;
            }
        }
    }
    if(type==-1) for(i=0;i<limit;i++) a[i]=a[i]*inv[limit]%MOD;
}
void init(){
    limit=1;cnt=0;int i;
    while(limit<=(n<<1)) limit<<=1,cnt++;
    for(i=0;i<limit;i++) r[i]=((r[i>>1]>>1)|((i&1)<<(cnt-1)));
    inv[1]=A[0]=B[0]=f[1]=finv[1]=1;A[1]=MOD-1;B[1]=n+1;
    for(i=2;i<=limit;i++) inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
    for(i=2;i<=limit;i++){
        f[i]=f[i-1]*i%MOD;
        finv[i]=finv[i-1]*inv[i]%MOD;
    }
}
int main(){
    n=read();
    init();int i;
    for(i=2;i<=n;i++) A[i]=(((i%2)?-1:1)*finv[i]+MOD)%MOD;
    for(i=2;i<=n;i++) B[i]=((qpow(i,n+1)-1)*inv[i-1]%MOD*finv[i])%MOD;
    ntt(A,1);ntt(B,1);
    for(i=0;i<limit;i++) C[i]=A[i]*B[i]%MOD;
    ntt(C,-1);
    ll ans=0;
    for(i=0;i<=n;i++) ans=(ans+qpow(2,i)*f[i]%MOD*C[i]%MOD)%MOD;
    printf("%lld\n",(ans+1)%MOD);//这里不知道为什么，一定要加个1，我也没有搞明白
}