BZOJ4555 [Tjoi2016&Heoi2016]求和

不错的题目，这里介绍两种方法。

第一种方法是强行化简这个式子。

首先S(i,j)当j>i时是0，所以原式可写成$\sum_{i=0}^n\sum_{j=0}^nS(i,j)*2^j*j!$

再考虑如何求$S(n,m)*m!$，它的意义是n个不同的球，放进m个不同的盒子里，盒子不允许空的方案数，求它的通项有很多方法，这里介绍比较简便的生成函数求法。

我们固定m，并定义多项式$A(x)=\sum_{i=0}^\infty A_i\frac{x^n}{n!}$的每一项系数$A_i$表示$S(i,m)*m!$

那$A(x)=(e^x-1)^m$（想一想，为什么）。

于是$S(n,m)*m!$即$\frac{x^n}{n!}$的系数就是$\sum_{k=0}^m(-1)^kC_m^k(m-k)^n$

原式即变为$\sum_{i=0}^n\sum_{j=0}^n2^j\sum_{k=0}^j(-1)^kC_j^k(j-k)^i$

变形得$\sum_{j=0}^n2^j*j!\sum_{k=0}^j\frac{(-1)^k}{k!}*\frac{\sum_{i=0}^n\;(j-k)^i}{(j-k)!}$

定义多项式$F(x)$的每一项$F_i=\frac{(-1)^i}{i!}$，定义多项式$G(x)$的每一项$G_i=\frac{\sum_{j=0}^ni^j}{i!}$，定义多项式$H(x)=F(x)*G(x)$

则$ans=\sum_{j=0}^n2^j*j!*H_j$

NTT一发就行了。

第二种方法是考虑原式的意义。

$F_i=\sum_{j=0}^iS(i,j)*2^j*j!$的意义是把n个不同的球，放进若干个不同的盒子了，盒子不允许空，每个盒子有两种状态的方案数。

枚举最后一个盒子的球数可得递推式$F_i=\sum_{j=1}^i2C_i^jF_{i-j}$

变形得$\frac{F_i}{i!}=\sum_{j=1}^i\frac 2{j!}*\frac{F_{i-j}}{(i-j)!}$

这是个卷积的形式，多项式求逆或者分治搞一搞就行了。

第一种做法代码：

#include <cstdio>
#include <algorithm>

typedef long long ll;
const int p=998244353,N=300000;
int n,m,l=-1,r[N];
ll ans,f[N],g[N],fr[N],ni[N];
ll pw(ll a,int b) {
    ll r=1;
    for(;b;b>>=1,a=a*a%p) if(b&1) r=r*a%p;
    return r;
}

void ntt(ll *a,int f) {
    for(int i=0;i<n;i++) if(r[i]>i) std::swap(a[i],a[r[i]]);
    for(int i=2;i<=n;i<<=1) {
        ll wn=pw(3,((p-1)/i*f+p-1)%(p-1)),m=i>>1;
        for(int j=0;j<n;j+=i) {
            ll w=1;
            for(int k=0;k<m;k++,w=w*wn%p) {
                ll x=a[j+k],y=a[j+k+m]*w%p;
                a[j+k]=(x+y)%p,a[j+k+m]=(x-y+p)%p;
            }
        }
    }
    if(f==-1) {
        ll ni=pw(n,p-2);
        for(int i=0;i<n;i++) a[i]=a[i]*ni%p;
    }
}

int main() {
    scanf("%d",&n),fr[0]=ni[0]=1;
    for(int i=1;i<=n;i++) fr[i]=fr[i-1]*i%p,ni[i]=pw(fr[i],p-2);
    for(int i=0;i<=n;i++) {
        if(i&1) f[i]=(p-ni[i])%p; else f[i]=ni[i];
        if(i==1) g[i]=(n+1)*ni[i]%p; else g[i]=(pw(i,n+1)-1+p)*pw(i-1+p,p-2)%p*ni[i]%p;
    }
    for(m=n,n=1;n<=m<<1;n<<=1) l++;
    for(int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|((i&1)<<l);
    ntt(f,1),ntt(g,1);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i]%p;
    ntt(f,-1);
    for(int i=0;i<=m;i++) ans=(ans+pw(2,i)*fr[i]%p*f[i])%p;
    printf("%lld",ans);
    return 0;
}