【題意】給定n，求Σi=0~nΣj=1~i s(i,j)*2^j*j!，n<=10^5。

【算法】生成函數+排列組合+多項式求逆

【題解】參考： [BZOJ4555][Tjoi2016&Heoi2016]求和-NTT-多項式求逆

$ans=\sum_{i=0}^{n}\sum_{j=0}^{i}s(i,j)*2^j*j!$

令$g(n)=\sum_{j=0}^{n}s(n,j)*2^j*j!$

則ans是Σg(i)，只要計算出g(i)的生成函數就可以統計答案。

g(n)可以理解為將n個數劃分成若幹集合，每個集合有2個屬性的排列數。基於此實際意義，通過枚舉第一個集合來遞推g(n)。

$g(n)=\sum_{i=1}^{n}2*C(n,i)*g(n-i)$

特別的，g(0)=1。

兩邊乘n!（令人窒息的操作），得

$\frac{g(n)}{n!}=\sum_{i=1}^{n}\frac{2}{i!}*\frac{g(n-i)}{(n-i)!}$

這已經是卷積的形式了：

$F(n)=\sum_{i=1}^{n}\frac{2}{n!}*x^i$

$G(n)=\sum_{i=0}^{n}\frac{g(n)}{n!}*x^i$

註意此時F*G卷積後，G(0)的位置是0，所以

$G(n)=F(n)G(n)+1$

移項得，

$G(n)=\frac{1}{1-F(n)}$

最後，多項式求逆即可，復雜度O(n log n)。

#include<cstdio>
#include
 
<cstring>
#include<algorithm>
using namespace std;
const int maxn=300010,MOD=998244353;
void gcd(int a,int b,int &x,int &y){
    if(!b){x=1;y=0;}else{gcd(b,a%b,y,x);y-=x*(a/b);}
}
int inv(int a){int x,y;gcd(a,MOD,x,y);return (x%MOD+MOD)%MOD;}
int power(int x,int k){
    int ans=1;
    
 
while(k){
        if(k&1)ans=1ll*ans*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1;
    }
    return ans;
}
namespace ntt{
    int o[maxn],oi[maxn],f[maxn];
    void init(int n){
        int x=1,p=power(3,(MOD-1)/n);
        for(int k=0;k<n;k++){
            o[k]=x;oi[k]=inv(o[k]);
            x=1ll*x*p%MOD;
        }
    }
    void transform(int *a,int n,int *o){
        int k=0;
        while((1<<k)<n)k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++)if((1<<j)&i)t|=(1<<(k-j-1));
            if(i<t)swap(a[i],a[t]);
        }
        for(int l=2;l<=n;l*=2){
            int m=l>>1;
            for(int *p=a;p!=a+n;p+=l){
                for(int i=0;i<m;i++){
                    int t=1ll*p[i+m]*o[n/l*i]%MOD;
                    p[i+m]=(p[i]-t+MOD)%MOD;
                    p[i]=(p[i]+t)%MOD;
                }
            }
        }
    }
    void dft(int *a,int n){transform(a,n,o);}
    void idft(int *a,int n){
        transform(a,n,oi);int v=inv(n);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*v%MOD;
    }
    void pinv(int *F,int *g,int n){
        if(n==1){g[0]=inv(F[0]);return;}
        pinv(F,g,n>>1);n<<=1;
        init(n);// 
        for(int i=0;i<n/2;i++)f[i]=F[i],f[i+n/2]=0;
        dft(f,n);dft(g,n);
        for(int i=0;i<n;i++)g[i]=1ll*g[i]*(2-1ll*f[i]*g[i]%MOD+MOD)%MOD;//1ll
        idft(g,n);for(int i=n/2;i<n;i++)g[i]=0;//
        
    }
}
int F[maxn],G[maxn],n,fac[maxn];
int main(){
    int n,N=1;
    scanf("%d",&n);n++;
    while(N<n)N*=2;
    fac[0]=1;
    for(int i=1;i<n;i++)fac[i]=1ll*fac[i-1]*i%MOD;
    for(int i=1;i<n;i++)F[i]=((-2*inv(fac[i]))%MOD+MOD)%MOD;
    F[0]++;
    ntt::pinv(F,G,N);
    int ans=0;
    for(int i=0;i<n;i++)ans=(ans+1ll*G[i]*fac[i]%MOD)%MOD;
    printf("%d",ans);
    return 0;
}

註意：多項式求逆過程中每次都要對不同的n進行一次預處理omega[]。

另一種做法

【算法】斯特林數+NTT

【題解】首先有第二類斯特林數的通項公式。

$s(n,m)=\frac{1}{m!}\sum_{k=0}^{m}(-1)^k*C(m,k)*(m-k)^n$

當斯特林數s(n,m)滿足m>n時，上述公式計算結果為0，所以第二個Σ可以擴展到n。

$ans=\sum_{i=0}^{n}\sum_{j=0}^{n}s(i,j)*s^j*j!$

代入第二類斯特林數公式。

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

通過組合數的分解，向卷積靠攏。

註意到Σi只對最後一個括號有影響，所以移動到最後。

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

這已經是標準的卷積形式（第二個函數分子是等比數列可以直接計算）。

使用NTT計算。

註意：

1.n=0，只有使0^0=1，斯特林數通項公式才能處理s(0,0)的情況。

2.n=1，等比數列求和公式不能處理Σ1^i即q=1的情況。

所以特殊地，g[0]=1，g[1]=n+1。先計算完再n++。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=300010,MOD=998244353;

int power(int x,int k){
    int ans=1;
    while(k){
        if(k&1)ans=1ll*ans*x%MOD;
        x=1ll*x*x%MOD;
        k>>=1;
    }
    return ans;
}
int inv(int x){return power(x,MOD-2);}
namespace ntt{
    int o[maxn],oi[maxn];
    void init(int n){
        int x=1,g=power(3,(MOD-1)/n);
        for(int k=0;k<n;k++){
            o[k]=x;oi[k]=inv(o[k]);
            x=1ll*x*g%MOD;
        }
    }
    void transform(int *a,int n,int *o){
        int k=0;
        while((1<<k)<n)k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++)if((1<<j)&i)t|=(1<<(k-j-1));
            if(i<t)swap(a[i],a[t]);
        }
        for(int l=2;l<=n;l*=2){
            int m=l>>1;
            for(int *p=a;p!=a+n;p+=l){
                for(int i=0;i<m;i++){
                    int t=1ll*p[i+m]*o[n/l*i]%MOD;
                    p[i+m]=(p[i]-t+MOD)%MOD;
                    p[i]=(p[i]+t)%MOD;
                }
            }
        }
    }
    void dft(int *a,int n){transform(a,n,o);}
    void idft(int *a,int n){
        transform(a,n,oi);
        int x=inv(n);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*x%MOD;
    }
}
int n,fac[maxn],f[maxn],g[maxn];
int main(){
    scanf("%d",&n);
    fac[0]=1;
    for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%MOD;
    for(int i=0;i<=n;i++)f[i]=1ll*((i&1)?MOD-1:1)*inv(fac[i])%MOD;
    for(int i=2;i<=n;i++)g[i]=1ll*(1-power(i,n+1)+MOD)*inv((1-i+MOD)%MOD)%MOD*inv(fac[i])%MOD;
    g[0]=1;g[1]=n+1;//
    n++;int N=1;//
    while(N<n+n)N*=2;
    ntt::init(N);
    ntt::dft(f,N);ntt::dft(g,N);
    for(int i=0;i<N;i++)f[i]=1ll*f[i]*g[i]%MOD;
    ntt::idft(f,N);
    int ans=0;
    for(int i=0;i<n;i++)ans=(ans+1ll*power(2,i)*fac[i]%MOD*f[i]%MOD)%MOD;
    printf("%d",ans);
    return 0;
}

【BZOJ】4555: [Tjoi2016&Heoi2016]求和 排列組合+多項式求逆 或 斯特林數+NTT