BZOJ 4555: [Tjoi2016&Heoi2016]求和

好久之前做的题了，算是我的NTT入门题了，然后我点开题单里的多项式，除了这题全TM是权限题当场去世

废话不多说直接推式子：

\[ans=\sum_{i=0}^n\sum_{j=0}^i \left\{_j^i\right\}\times 2^j\times j!\\=\sum_{i=0}^n\sum_{j=0}^n \left\{_j^i\right\}\times 2^j\times j!\\=\sum_{j=0}^n 2^j\times j!\times \sum_{i=0}^n\left\{_j^i\right\}\\=\sum_{j=0}^n 2^j\times j!\times \sum_{i=0}^n\sum_{k=0}^j (-1)^k\times \frac{1}{k!(j-k)!}\times (j-k)^i\\=\sum_{j=0}^n 2^j\times j!\times \sum_{k=0}^j (-1)^k\times \frac{1}{k!(j-k)!}\times \sum_{i=0}^n(j-k)^i\\=\sum_{j=0}^n 2^j\times j!\times \sum_{k=0}^j \frac{(-1)^k}{k!}\times \frac{\sum_{i=0}^n(j-k)^i}{j-k}\\ \]

注意上面的\(\left\{_m^n\right\}=\sum_{i=0}^m (-1)^i\times \frac{1}{i!(m-i)!}\times (m-i)^n\)是第二类斯特林数的通项公式，可以通过组合意义推出

推到这里我们发现式子已经化成了卷积的形式，设\(A(x)=\frac{(-1)^x}{x!},B(x)=\frac{\sum_{i=0}^n x^i}{x!}=\frac{x^{n+1}-1}{(x-1)x!}\)，并且\(B(1)=n+1\)，然后直接大力NTT就好了，复杂度\(O(n\log n)\)

PS：据说有\(O(n)\)的生成函数的做法但是咱太弱了不会啊QAQ

#include<cstdio>
#include<cmath>
#define RI register int
const int N=100005,mod=998244353;
int n,ans,pw[N],fact[N],inv[N],A[N<<2],B[N<<2],lim;
inline void inc(int &x,int y)
{
    if ((x+=y)>=mod) x-=mod;
}
inline void dec(int &x,int y)
{
    if ((x-=y)<0) x+=mod;
}
inline int quick_pow(int x,int p,int mul=1)
{
    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
class NTT_Solver
{
    private:
        int rev[N<<2],p;
        inline void swap(int &x,int &y)
        {
            int t=x; x=y; y=t;
        }
    public:
        inline void init(int n)
        {
            for (lim=1;lim<=n;lim<<=1,++p);
            for (RI i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<p-1);
        }
        inline void NTT(int *f,int opt)
        {
            RI i,j,k; for (i=0;i<lim;++i) if (i<rev[i]) swap(f[i],f[rev[i]]);
            for (i=1;i<lim;i<<=1)
            {
                int m=i<<1,D=quick_pow(3,opt==1?(mod-1)/m:mod-1-(mod-1)/m);
                for (j=0;j<lim;j+=m)
                {
                    int W=1; for (k=0;k<i;++k,W=1LL*W*D%mod)
                    {
                        int x=f[j+k],y=1LL*f[i+j+k]*W%mod;
                        f[j+k]=f[i+j+k]=x; inc(f[j+k],y); dec(f[i+j+k],y);
                    }
                }
            }
            if (opt==-1)
            {
                int Inv=quick_pow(lim,mod-2);
                for (i=0;i<lim;++i) f[i]=1LL*f[i]*Inv%mod;
            }
        }
}P;
inline void init(int n)
{
    RI i; fact[0]=inv[0]=pw[0]=A[0]=B[0]=1; B[1]=n+1;
    for (i=1;i<=n;++i) pw[i]=2LL*pw[i-1]%mod;
    for (i=1;i<=n;++i) fact[i]=1LL*fact[i-1]*i%mod;
    for (inv[n]=quick_pow(fact[n],mod-2),i=n-1;i;--i)
    inv[i]=1LL*inv[i+1]*(i+1)%mod; for (i=1;i<=n;++i)
    if (i&1) A[i]=mod-inv[i]; else A[i]=inv[i];
    for (i=2;i<=n;++i) B[i]=quick_pow(i,n+1),dec(B[i],1),
    B[i]=1LL*B[i]*quick_pow(i-1,mod-2)%mod*inv[i]%mod;
}
int main()
{
    RI i; scanf("%d",&n); init(n); P.init(n<<1);
    for (P.NTT(A,1),P.NTT(B,1),i=0;i<lim;++i) A[i]=1LL*A[i]*B[i]%mod;
    for (P.NTT(A,-1),i=0;i<=n;++i)
    inc(ans,1LL*pw[i]*fact[i]%mod*A[i]%mod);
    return printf("%d",ans),0;
}