[BZOJ 4555][Tjoi2016&Heoi2016]求和

题意

给定 $n$ , 求下式的值:

$$ f(n)= sum_{i=0}^nsum_{j=0}^iegin{Bmatrix}i\ jend{Bmatrix} imes 2^j imes j!$$

题解

这题比较神仙...

那么我们可以思考如何来求一个比较简单的转移式.

首先我们发现, $f(n)$ 表达式中的第一重和式包含了 $f(n-1)$, 那么我们对 $f$ 的值做差分, 于是我们有 $f(n)-f(n-1)=sumlimits_{i=0}^negin{Bmatrix}i\ jend{Bmatrix} imes 2^j imes j!$ .设 $g(n)=f(n)-f(n-1)$.

而我们知道第二类斯特林数的组合意义是在 $i$ 个物品中选 $j$ 个子集的方案数. 那么后面的 $j!$ 的组合意义相当于让划分的子集带标号, $2^j$ 的组合意义相当于每个子集贡献两种状态.

知道了表达式的组合意义, 我们可以尝试用朴素DP来解它. 我们枚举一下最后一个子集中的元素个数为 $i$ (因为子集带标号, 我们相当于枚举最后一个), 那么我们就有递推式:

$$g(n)=[n=0]+sum_{i=1}^n {nchoose i} g(n-i) imes 2$$

最后乘 $2$ 是因为每个集合会贡献两种状态, $[n=0]$ 是边界条件.

观察这个式子, 发现是二项卷积的形式, 我们果断上EGF来解决. 设 $G(x)$ 是 $g(n)$ 的EGF, $H(x)$ 是数列 $langle 0,2,2,2,dots angle$ 的EGF, 那么我们有:

$$ G(x)=G(x)H(x)+1 $$

注意到 $H(x)$ 代表的数列的第 $0$ 项是 $0$, 因为 $g(n)$ 的递推式中求和下指标是 $1$.

那么我们移项之后就可以开心地得到:

$$ G(x)=frac 1 {1-H(x)} $$

NTT爆算一发就可以了

代码实现

记得EGF要除以阶乘...出答案的时候还得乘回来...

#include <bits/stdc++.h>

const int G=3;
const int DFT=1;
const int IDFT=-1;
const int MAXN=550000;
const int MOD=998244353;
const int INV2=(MOD+1)>>1;
const int PHI=MOD-1;

typedef std::vector<int> Poly;

Poly Sqrt(Poly);
void Read(Poly&);
Poly Inverse(Poly);
Poly Ln(const Poly&);
Poly Exp(const Poly&);
void Print(const Poly&);
void NTT(Poly&,int,int);
Poly Pow(const Poly&,int);
Poly Integral(const Poly&);
Poly Derivative(const Poly&);
Poly operator*(Poly,Poly);
Poly operator/(Poly,Poly);
Poly operator%(Poly,Poly);
Poly operator+(const Poly&,const Poly&);
Poly operator-(const Poly&,const Poly&);

int rev[MAXN];

int NTTPre(int);
int Pow(int,int,int);

int main(){
    int n;
    scanf("%d",&n);
    Poly h(n+1),one(1,1);
    h[1]=2;
    for(int i=2;i<=n;i++)
        h[i]=1ll*h[i-1]*Pow(i,MOD-2,MOD)%MOD;
    Poly g(Inverse(one-h));
    int ans=0;
    int fact=1;
    for(int i=0;i<=n;i++){
        (ans+=1ll*g[i]*fact%MOD)%=MOD;
        fact=1ll*fact*(i+1)%MOD;
    }
    printf("%d
",ans);
    return 0;
}

void Read(Poly& a){
    for(auto& i:a)
        scanf("%d",&i);
}

void Print(const Poly& a){
    for(auto i:a)
        printf("%d ",i);
    puts("");
}

Poly Pow(const Poly& a,int k){
    Poly log=Ln(a);
    for(auto& i:log)
        i=1ll*i*k%MOD;
    return Exp(log);
}

Poly Sqrt(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,int(sqrt(a[0])));
    else{
        Poly b=a;
        b.resize((len+1)>>1);
        b=Sqrt(b);
        b.resize(len);
        Poly inv=Inverse(b);
        int bln=NTTPre(inv.size()+a.size());
        NTT(a,bln,DFT);
        NTT(inv,bln,DFT);
        for(int i=0;i<bln;i++)
            a[i]=1ll*a[i]*INV2%MOD*inv[i]%MOD;
        NTT(a,bln,IDFT);
        for(int i=0;i<len;i++)
            b[i]=(1ll*b[i]*INV2%MOD+a[i])%MOD;
        return b;
    }
}

Poly Exp(const Poly& a){
    size_t len=1;
    Poly ans(1,1),one(1,1);
    while(len<(a.size()<<1)){
        len<<=1;
        Poly b=a;
        b.resize(len);
        ans=ans*(one-Ln(ans)+b);
        ans.resize(len);
    }
    ans.resize(a.size());
    return ans;
}

Poly Ln(const Poly& a){
    Poly ans=Integral(Derivative(a)*Inverse(a));
    ans.resize(a.size());
    return ans;
}

Poly Integral(const Poly& a){
    int len=a.size();
    Poly ans(len+1);
    for(int i=1;i<len;i++)
        ans[i]=1ll*a[i-1]*Pow(i,MOD-2,MOD)%MOD;
    return ans;
}

Poly Derivative(const Poly& a){
    int len=a.size();
    Poly ans(len-1);
    for(int i=1;i<len;i++)
        ans[i-1]=1ll*a[i]*i%MOD;
    return ans;
}

Poly operator/(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans(1);
    if(n>=m){
        std::reverse(a.begin(),a.end());
        std::reverse(b.begin(),b.end());
        b.resize(n-m+1);
        ans=Inverse(b)*a;
        ans.resize(n-m+1);
        std::reverse(ans.begin(),ans.end());
    }
    return ans;
}

Poly operator%(Poly a,Poly b){
    int n=a.size()-1,m=b.size()-1;
    Poly ans;
    if(n<m)
        ans=a;
    else
        ans=a-(a/b)*b;
    ans.resize(m);
    return ans;
}

Poly operator*(Poly a,Poly b){
    int len=a.size()+b.size()-1;
    int bln=NTTPre(len);
    NTT(a,bln,DFT);
    NTT(b,bln,DFT);
    for(int i=0;i<bln;i++)
        a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,bln,IDFT);
    a.resize(len);
    return a;
}

Poly operator+(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+b[i])%MOD;
    return ans;
}

Poly operator-(const Poly& a,const Poly& b){
    Poly ans(std::max(a.size(),b.size()));
    std::copy(a.begin(),a.end(),ans.begin());
    for(size_t i=0;i<b.size();i++)
        ans[i]=(ans[i]+MOD-b[i])%MOD;
    return ans;
}

Poly Inverse(Poly a){
    int len=a.size();
    if(len==1)
        return Poly(1,Pow(a[0],MOD-2,MOD));
    else{
        Poly b(a);
        b.resize((len+1)>>1);
        b=Inverse(b);
        int bln=NTTPre(b.size()*2+a.size());
        NTT(a,bln,DFT);
        NTT(b,bln,DFT);
        for(int i=0;i<bln;i++)
            b[i]=(2ll*b[i]%MOD-1ll*b[i]*b[i]%MOD*a[i]%MOD+MOD)%MOD;
        NTT(b,bln,IDFT);
        b.resize(len);
        return b;
    }
}

void NTT(Poly& a,int len,int opt){
    a.resize(len);
    for(int i=0;i<len;i++)
        if(rev[i]>i)
            std::swap(a[i],a[rev[i]]);
    for(int i=1;i<len;i<<=1){
        int step=i<<1;
        int wn=Pow(G,(PHI+opt*PHI/step)%PHI,MOD);
        for(int j=0;j<len;j+=step){
            int w=1;
            for(int k=0;k<i;k++,w=1ll*w*wn%MOD){
                int x=a[j+k];
                int y=1ll*a[j+k+i]*w%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
    if(opt==IDFT){
        int inv=Pow(len,MOD-2,MOD);
        for(int i=0;i<len;i++)
            a[i]=1ll*a[i]*inv%MOD;
    }
}

int NTTPre(int n){
    int bln=1,bct=0;
    while(bln<n){
        bln<<=1;
        ++bct;
    }
    for(int i=0;i<bln;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1));
    return bln;
}

inline int Pow(int a,int n,int p){
    int ans=1;
    while(n>0){
        if(n&1)
            ans=1ll*a*ans%p;
        a=1ll*a*a%p;
        n>>=1;
    }
    return ans;
}