洛谷 P4091 [HEOI2016/TJOI2016]求和 任意模数fft

题目描述

在2016年，佳媛姐姐刚刚学习了第二类斯特林数，非常开心。

现在他想计算这样一个函数的值:

$f(n)=\sum_{i=0}^n\sum_{j=0}^i S(i,j)\times 2^j \times (j!)$

$S(i, j)$表示第二类斯特林数，递推公式为:

$S(i,j)=j×S(i−1,j)+S(i−1,j−1),1≤j≤i−1$ 。

边界条件为：$S(i,i)=1(0≤i),S(i,0)=0(1≤i)$

你能帮帮他吗?

输入输出格式

输入格式：
输入只有一个正整数

输出格式：
输出$f(n)$。由于结果会很大，输出$f(n)$对$998244353(7 × 17 × 2^{23} + 1)$取模的结果即可。

输入输出样例

输入样例#1：
3
输出样例#1：
87
说明

对于50%数据$1≤n≤5000$
对于100%数据$1≤n≤100000$

分析：
显然就是推式子了。
其中，第二类斯特林数的通项公式为

$$S(i,j)=\frac{1}{j!}\sum_{k=0}^{j}(-1)^k*\binom{j}{k}*(j-k)^i$$
证明使用生成函数，至于记的话，可以记住只有$(j-k)$的指数是$i$，其他全是$j$和$k$。
所以，
$$ans=\sum_{i=0}^{n}\sum_{j=0}^{n}s(i,j)*j!*2^j$$
$$=\sum_{j=0}^{n}j!*2^j\sum_{i=0}^{n}s(i,j)$$
$$=\sum_{j=0}^{n}j!*2^j\sum_{i=0}^{n}\frac{1}{j!}\sum_{k=0}^{j}(-1)^k*\binom{j}{k}*(j-k)^i$$
$$=\sum_{j=0}^{n}j!*2^j\sum_{i=0}^{n}\sum_{k=0}^{j}\frac{(-1)^k}{k!}*\frac{(j-k)^i}{(j-k)!}$$
$$=\sum_{j=0}^{n}j!*2^j\sum_{k=0}^{j}\frac{(-1)^k}{k!}*\frac{\sum_{i=0}^{n}(j-k)^i}{(j-k)!}$$
显然是一个卷积的形式，直接任意模数fft即可。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const LL p=998244353;
const int maxn=300007;
const double pi=acos(-1);

using namespace std;

struct rec{
    double x,y;
};

rec operator +(rec a,rec b)
{
    return (rec){a.x+b.x,a.y+b.y};
}

rec operator -(rec a,rec b)
{
    return (rec){a.x-b.x,a.y-b.y};
}

rec operator *(rec a,rec b)
{
    return (rec){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}

rec operator !(rec a)
{
    return (rec){a.x,-a.y};
}

LL n,ans,len;
rec a[maxn],b[maxn],w[maxn],dfta[maxn],dftb[maxn],dftc[maxn],dftd[maxn];
LL bit[maxn],jc[maxn],inv[maxn],f[maxn],g[maxn],r[maxn];

LL power(LL x,LL y)
{
    if (y==1) return x;
    LL c=power(x,y/2);
    c=(c*c)%p;
    if (y%2) c=(c*x)%p;
    return c; 
}

void fft(rec *a,LL f)
{
    for (LL i=0;i<len;i++)
    {
        if (i<r[i]) swap(a[i],a[r[i]]);
    }
    w[0]=(rec){1,0};
    for (LL i=2;i<=len;i*=2)
    {
        rec wn=(rec){cos(2*pi/i),f*sin(2*pi/i)};
        for (LL j=i/2;j>=0;j-=2) w[j]=w[j/2];
        for (LL j=1;j<i/2;j+=2) w[j]=w[j-1]*wn;
        for (LL j=0;j<len;j+=i)
        {
            for (LL k=0;k<i/2;k++)
            {
                rec u=a[j+k],v=a[j+k+i/2]*w[k];
                a[j+k]=u+v;
                a[j+k+i/2]=u-v;
            }
        }
    }
}

void init(LL len)
{
    LL k=trunc(log(len+0.5)/log(2));
    for (LL i=0;i<len;i++)
    {
        r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));
    }
}

void FFT(LL *x,LL *y,LL *z,LL n,LL m)
{   
    len=1;
    while (len<(n+m-1)) len*=2;
    init(len);
    for (LL i=0;i<len;i++)
    {
        LL A,B;
        if (i<n) A=x[i]%p; else A=0;
        if (i<m) B=y[i]%p; else B=0;
        a[i]=(rec){A>>15,A&32767};
        b[i]=(rec){B>>15,B&32767};
    }
    fft(a,1); fft(b,1);
    for (LL i=0;i<len;i++)
    {
        LL j=(len-i)&(len-1);
        rec da,db,dc,dd;
        da=(a[i]+(!a[j]))*(rec){0.5,0};
        db=(a[i]-(!a[j]))*(rec){0,-0.5};
        dc=(b[i]+(!b[j]))*(rec){0.5,0};
        dd=(b[i]-(!b[j]))*(rec){0,-0.5};
        dfta[i]=da*dc;
        dftb[i]=da*dd;
        dftc[i]=db*dc;
        dftd[i]=db*dd;
    }
    for (LL i=0;i<len;i++)
    {
        a[i]=dfta[i]+dftb[i]*(rec){0,1};
        b[i]=dftc[i]+dftd[i]*(rec){0,1};
    }
    fft(a,-1); fft(b,-1);
    for (LL i=0;i<len;i++)
    {
        LL da,db,dc,dd;
        da=(LL)(a[i].x/len+0.5)%p;
        db=(LL)(a[i].y/len+0.5)%p;
        dc=(LL)(b[i].x/len+0.5)%p;
        dd=(LL)(b[i].y/len+0.5)%p;
        z[i]=((da<<30)%p+((db+dc)<<15)%p+dd)%p;
    }
}

int main()
{
    scanf("%lld",&n);
    bit[0]=1; jc[0]=1;  
    for (LL i=1;i<=n;i++)
    {
        bit[i]=(bit[i-1]*2)%p;
        jc[i]=(jc[i-1]*i)%p;
    }   
    inv[n]=power(jc[n],p-2);
    for (LL i=n;i>0;i--) inv[i-1]=(inv[i]*i)%p;
    g[0]=1; f[0]=1;
    for (LL i=1;i<=n;i++)
    {
        if (i%2) f[i]=(p-1)*inv[i]%p;
            else f[i]=inv[i];
        g[i]=(power(i,n+1)+p-1)%p*power(i-1,p-2)%p*inv[i]%p;
    }       
    g[1]=n+1;   
    FFT(f,g,f,n+1,n+1);
    for (LL i=0;i<=n;i++) ans=(ans+bit[i]*jc[i]%p*f[i]%p)%p;
    printf("%lld",ans);
}