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

【题意】给定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;
}