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

View Code