BZOJ 4555 HEOI2016/TJOI2016求和

Problem

BZOJ

Solution

发现我的思路比较清奇？可能是我数学太差了想不到

首先可以用容斥推出第二类斯特林数的通项公式(没见过就是真·推不出来)

$$S_n^m=\frac 1 {m!}\sum_{i=0}^m (-1)^iC_m^i (m-i)^n=\sum_{i=0}^m (-1)^i(i!)\frac {(m-i)^n} {(m-i)!}$$
这是一个卷积的形式，那么我们可以利用NTT在 $O(n\log n)$的时间求出第n行的斯特林数。

我们考虑把$2^j\times (j!)$看作一个系数，根据递推公式，把斯特林数往下一行推。

然后如果你手动模拟一下，就像这样：

i\j	0	1	2
0	1
1	0	1
2	1	3	8

i\j	0	1	2	3
0	1
1	0	1
2	0	0	4
3	1	3	10	48

i\j	0	1	2	3	4
0	1
1	0	1
2	0	0	4
3	0	0	0	30
4	1	3	10	54	384

其中，$S_i^i$的系数会有剩余，但我们知道它恒等于1，所以只需要知道它的系数即可。

又容易发现，这些系数是有规律的，我们可以用递推来模拟做这样的一个事情，复杂度$O(n)$。转移方程可以参考代码，其中a数组为$s_i^i$的系数，b数组为最后一行的系数，c为第i列向第i+1列的转移系数。唯一需要注意的是第n项的系数是没有得到加成的，还是$2^n(n!)$。

时间复杂度$O(n\log n)$。

Code

#include <algorithm>
#include <cstring>
#include <cstdio>
#define rg register
using namespace std;
typedef long long ll;
const int maxn=100010,mod=998244353,G=3;
template <typename Tp> inline void getmin(Tp &x,Tp y){if(y<x) x=y;}
template <typename Tp> inline void getmax(Tp &x,Tp y){if(y>x) x=y;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;char ch=getchar();int f=0;
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int n,m,tmp,N=1,l,ans,f[maxn<<2],g[maxn<<2],rev[maxn<<2];
int fac[maxn],inv[maxn],mi[maxn],a[maxn],b[maxn],c[maxn];
inline int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
int power(int x,int y)
{
    int res=1;
    for(;y;y>>=1,x=(ll)x*x%mod)
      if(y&1)
        res=(ll)res*x%mod;
    return res;
}
void init()
{
    fac[0]=mi[0]=1;
    for(rg int i=1;i<=n;i++)
      fac[i]=(ll)fac[i-1]*i%mod,mi[i]=pls(mi[i-1],mi[i-1]);
    inv[n]=power(fac[n],mod-2);
    for(rg int i=n-1;~i;i--) inv[i]=(ll)inv[i+1]*(i+1)%mod;
    a[0]=1;b[0]=1;c[0]=1;
    for(rg int i=1;i<=n;i++)
    {
        a[i]=dec((ll)fac[i]*mi[i]%mod,(ll)i*c[i-1]%mod);
        b[i]=pls((ll)fac[i]*mi[i]%mod,c[i-1]%mod);
        c[i]=dec(b[i],(ll)i*c[i-1]%mod);
    }
    for(rg int i=0;i<n;i++) ans=pls(ans,a[i]);
}
void ntt(int *a,int f)
{
    for(rg int i=0;i<N;i++)
      if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(rg int i=1;i<N;i<<=1)
    {
        int gn=power(G,(mod-1)/(i<<1));
        for(rg int j=0;j<N;j+=(i<<1))
        {
            int g=1;
            for(rg int k=0;k<i;k++,g=(ll)g*gn%mod)
            {
                int x=a[j+k],y=(ll)g*a[j+k+i]%mod;
                a[j+k]=pls(x,y);a[j+k+i]=dec(x,y);
            }
        }
    }
    if(f==-1)
    {
        int inv=power(N,mod-2);reverse(a+1,a+N);
        for(int i=0;i<N;i++) a[i]=(ll)a[i]*inv%mod;
    }
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif
    read(n);
    init();
    while(N<(n<<1)) N<<=1,l++;
    for(rg int i=1;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for(rg int i=0;i<n;i++)
    {
        f[i]=inv[i];
        if(i&1) f[i]=dec(0,f[i]);
        g[i]=(ll)power(i,n)*inv[i]%mod;
    }
    ntt(f,1);ntt(g,1);
    for(rg int i=0;i<N;i++) f[i]=(ll)f[i]*g[i]%mod;
    ntt(f,-1);
    for(rg int i=0;i<n;i++) ans=pls(ans,(ll)f[i]*b[i]%mod);
    ans=pls(ans,(ll)mi[n]*fac[n]%mod);
    printf("%d\n",ans);
    return 0;
}