题解 P3711 【仓鼠的数学题】

· · 题解

题意

设 S_k(x)=\sum\limits_{i=0}^{x}i^k，给定正整数 n 和序列 a，求：

\sum\limits_{k=0}^{n}S_k(x)a_k \texttt{Data Range:}1\leq n\leq 2\times 10^5

题解

伯努利数板子题。

假设 S_k(x)=\sum\limits_{i=0}^{x-1}i^k，之后只要把所有的 x 换成 x+1 即可，这么做是方便用伯努利数表示，那么有

S_k(x)=\frac{1}{k+1}\sum\limits_{i=0}^{k}\binom{k+1}{i}B_{i}x^{k+1-i}

这个并不是很好看，所以考虑枚举 k-i：

S_k(x)=\frac{1}{k+1}\sum\limits_{i=0}^{k}\binom{k+1}{k-i}B_{k-i}x^{i+1}

利用二项式的对称性：

S_k(x)=\frac{1}{k+1}\sum\limits_{i=0}^{k}\binom{k+1}{i+1}B_{k-i}x^{i+1}

于是答案为

\sum\limits_{k=0}^{n}\frac{a_k}{k+1}\sum\limits_{i=0}^{k}\binom{k+1}{i+1}B_{k-i}x^{i+1}

显然交换求和顺序

\sum\limits_{i=0}^{n}x^{i+1}\sum\limits_{k=i}^{n}\frac{a_k}{k+1}\binom{k+1}{i+1}B_{k-i}

拆掉二项式

\sum\limits_{i=0}^{n}\frac{x^{i+1}}{(i+1)!}\sum\limits_{k=i}^{n}a_kk!\frac{B_{k-i}}{(k-i)!}

直接减法卷积即可，但是这里的 x 实际上是 x+1，于是我们考虑下面的问题，求

\sum\limits_{i=0}^{n}f_i(x+1)^i

直接利用二项式定理：

\sum\limits_{i=0}^{n}f_i\sum\limits_{j=0}^{i}\binom{i}{j}x^j

交换求和顺序之后减法卷积即可。

代码

#include<bits/stdc++.h>
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=524291,MOD=998244353;
ll n;
ll omgs[MAXN],fact[MAXN],finv[MAXN],f[MAXN],g[MAXN],h[MAXN];
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
}
inline ll qpow(ll base,ll exponent)
{
    ll res=1;
    while(exponent)
    {
        if(exponent&1)
        {
            res=(li)res*base%MOD;
        }
        base=(li)base*base%MOD,exponent>>=1;
    }
    return res;
}
inline void setupOmg(ll cnt)
{
    ll limit=__lg(cnt)-1;
    omgs[0]=1,omgs[1<<limit]=qpow(31,1<<(21-limit));
    for(register int i=limit;i;i--)
    {
        omgs[1<<(i-1)]=(li)omgs[1<<i]*omgs[1<<i]%MOD;
    }
    for(register int i=1;i<cnt;i++)
    {
        omgs[i]=(li)omgs[i&(i-1)]*omgs[i&-i]%MOD;
    }
}
inline void setup(ll cnt)
{
    fact[0]=fact[1]=finv[0]=1;
    for(register int i=1;i<=cnt;i++)
    {
        fact[i]=(li)fact[i-1]*i%MOD;
    }
    finv[cnt]=qpow(fact[cnt],MOD-2);
    for(register int i=cnt-1;i;i--)
    {
        finv[i]=(li)finv[i+1]*(i+1)%MOD;
    }
}
inline ll& reduce(ll &x)
{
    return x+=x>>31&MOD;
}
inline void DIF(ll *cp,ll cnt)
{
    ll lim=__lg(cnt),len=cnt>>1,x;
    for(register int i=0;i<lim;i++,len>>=1)
    {
        for(register int *j=cp,*omg=omgs;j!=cp+cnt;j+=len<<1,omg++)
        {
            for(register int *k=j;k!=j+len;k++)
            {
                x=(li)*omg*k[len]%MOD,reduce(k[len]=*k-x),reduce(*k+=x-MOD);
            }
        }
    }
}
inline void DIT(ll *cp,ll cnt)
{
    ll lim=__lg(cnt),len=1,x,invl;
    for(register int i=0;i<lim;i++,len<<=1)
    {
        for(register int *j=cp,*omg=omgs;j!=cp+cnt;j+=len<<1,omg++)
        {
            for(register int *k=j;k!=j+len;k++)
            {
                reduce(x=*k+k[len]-MOD);
                k[len]=(li)(*k-k[len]+MOD)**omg%MOD,*k=x;
            }
        }
    }
    reverse(cp+1,cp+cnt),invl=MOD-(MOD-1)/cnt;
    for(register int i=0;i<cnt;i++)
    {
        cp[i]=(li)cp[i]*invl%MOD;
    }
}
inline void NTT(ll *cp,ll cnt,ll inv)
{
    inv==1?DIF(cp,cnt):DIT(cp,cnt);
}
inline void conv(ll fd,ll *f,ll *g,ll *res)
{
    ll cnt=1,limit=-1;
    while(cnt<(fd<<1))
    {
        cnt<<=1,limit++;
    }
    for(register int i=fd;i<cnt;i++)
    {
        f[i]=g[i]=0;
    }
    NTT(f,cnt,1),NTT(g,cnt,1);
    for(register int i=0;i<cnt;i++)
    {
        res[i]=(li)f[i]*g[i]%MOD;
    }
    NTT(res,cnt,-1);
}
inline void inv(ll fd,ll *f,ll *res)
{
    static ll tmp[MAXN],tmpr[MAXN];
    ll cnt=2,limit=0;
    tmpr[0]=res[0]=qpow(f[0],MOD-2);
    while(cnt<(fd<<1))
    {
        for(register int i=0;i<cnt;i++)
        {
            tmp[i]=f[i],tmpr[i]=i<(cnt>>1)?res[i]:0;
        }
        NTT(tmp,cnt,1),NTT(tmpr,cnt,1);
        for(register int i=0;i<cnt;i++)
        {
            tmp[i]=(li)tmp[i]*tmpr[i]%MOD;
        }
        NTT(tmp,cnt,-1),memset(tmp,0,cnt<<1),NTT(tmp,cnt,1);
        for(register int i=0;i<cnt;i++)
        {
            tmp[i]=(li)tmp[i]*tmpr[i]%MOD;
        }
        NTT(tmp,cnt,-1);
        for(register int i=(cnt>>1);i<cnt;i++)
        {
            res[i]=!tmp[i]?0:MOD-tmp[i];
        }
        cnt<<=1,limit++;
    }
    for(register int i=fd;i<cnt;i++)
    {
        res[i]=0;
    }
}
int main()
{
    setup((n=read())+5),setupOmg(524288);
    for(register int i=0;i<=n;i++)
    {
        f[i]=(li)read()*fact[i]%MOD,h[i]=finv[i+1];
    }
    inv(n+1,h,g),reverse(f,f+n+1),conv(n+1,f,g,f),reverse(f,f+n+1),g[0]=0;
    for(register int i=0;i<=n+1;i++)
    {
        g[i+1]=f[i],h[i]=finv[i];
    }
    reverse(g,g+n+2),conv(n+2,g,h,f),reverse(f,f+n+2);
    for(register int i=0;i<=n+1;i++)
    {
        printf("%d ",(li)f[i]*finv[i]%MOD);
    }
}