题解 P3711 【仓鼠的数学题】
shadowice1984
2018-12-12 09:19:23
不知道为什么这题把自然数幂和函数的求和上限额外加了一个1,这样的话我们可能需要用二项式定理重新把原来的式子展开一遍……
### 前置芝士:伯努利数
如果不会伯努利数的话可以看一下[这篇博客](https://www.luogu.org/blog/ShadowassIIXVIIIIV/guan-yu-bo-nu-li-shuo-zhuai-hua-zi-ran-shuo-mi-hu-gong-shi-di-zheng-mi)学习一下
### 前置芝士:多项式求逆
蛤?不会多项式求逆,可以出门左转你站模板区去学习一下
### 前置芝士:快速处理伯努利数
我们使用伯努利数的指数型生成函数来快速的预处理它
具体来讲我们有这样一个和式
$$\sum_{i}\frac{B(i)}{i!}z^{i}=\frac{z}{e^{z}-1}$$
等会显然左边是一个无限和式,这东西显然没有办法计算啊
没关系啊,我们把$e^{z}$大力泰勒展开一下
$$e^{z}=\sum_{i}\frac{z^{i}}{i!}$$
那么我们要计算的式子就是
$$\frac{z}{\sum_{i}\frac{z^{i}}{i!}-1}$$
然后我们发现一个很尴尬的事实,这东西分母常数项是0求不了逆……
还好分子是$z$所以我们上下同时除一个$z$
$$\frac{1}{\sum_{i}\frac{z^{i}}{(i+1)!}}$$
那这样的话我们就可以愉快的多项式求逆然后把伯努利数的前$n$项刷出来了~
## 本题题解
众所周知,标准的自然数幂和函数$S(n,k)$是等于$\sum_{i=0}^{n-1}i^{k}$的但是这题我们的循环上限变成了$n$这就十分的难受了
所以设答案函数是$f(n)$,我们先计算比较好算的$f(n-1)$然后再用一些科技把答案的系数展开就好了
那么在以下的讨论当中我们认为$S(n,k)$代表$\sum_{i=0}^{n-1}i^{k}$而不是$\sum_{i=0}^{n}i^{k}$
题目中让我们求
$$\sum_{k=0}^{n}S(n,k)a(k)$$
强行用伯努利数展开一波可以得到
$$\sum_{k=0}^{n}a(k)\frac{1}{k+1}\sum_{i=0}^{k}{k+1 \choose i}B(i)n^{k+1-i}$$
接下来我们翻转后半部分的求和顺序
$$\sum_{k=0}^{n}a(k)\frac{1}{k+1}\sum_{i=1}^{k+1}{k+1 \choose k+1-i}B(k+1-i)n^{i}$$
让我们把组合数拆了搞成阶乘的形式
$$\sum_{k=0}^{n}a(k)k!\sum_{i=1}^{k+1}\frac{B(k+1-i)}{(k+1-i)!}\frac{n^{i}}{i!}$$
接下来又是喜闻乐见的交换求和号环节了
$$\sum_{i=1}^{n+1}\frac{n^{i}}{i!}\sum_{k=i-1}^{n}a(k)k!\frac{B(k+1-i)}{(k+1-i)!}$$
如果我们令$f,g$分别表示这两个式子的话
$$f(i)=\frac{B(i)}{i!}$$
$$g(i)=a(i)i!$$
式子会被写的好看一点
$$\sum_{i=1}^{n+1}\frac{n^i}{i!}\sum_{p-q=i-1}f(p)g(q)$$
似乎是个差卷积,众所周知我们的ntt只能处理和的卷积,所以让我们把g函数翻转一下
$$g_{r}(x)=g(n-x)$$
那么原来的式子就会变成
$$\sum_{i=1}^{n+1}\frac{n^i}{i!}\sum_{p+q=n+(i-1)}f(p)g_{r}(q)$$
这样的话我们一波ntt然后除上一个阶乘就可以求出$f(n-1)$第i项系数了,让我们设这个东西为$c_{i}$接下来我们要做的是求出这个多项式的系数
$$\sum_{i=1}^{n+1}c_{i}(x+1)^{i}$$
做法简单粗暴,直接二项式定理展开一下
$$\sum_{i=1}^{n+1}c_{i}\sum_{j=0}^{i}{i\choose j}x^{j}$$
让我们继续交换求和号,由于$c_{0}=0$所以我们无需担心求和下界的问题
$$\sum_{j=0}^{n+1}x^{j}\sum_{i=j}^{n+1}{i\choose j}c_{i}$$
把组合数用阶乘拆了可以得到
$$\sum_{j=0}^{n+1}x^{j}\frac{1}{j!}\sum_{i=j}^{n+1}c_{i}i!×\frac{1}{(i-j)!}$$
设$f,g$分别为
$$f(i)=c_{i}i!$$
$$g(i)=\frac{1}{i!}$$
那么我们会得到这样一个式子
$$\sum_{j=0}^{n+1}x^{j}\frac{1}{j!}\sum_{p-q=j}f(p)g(q)$$
看起来又是烦人的差卷积,没关系接着翻转数列,令
$$g_{r}(x)=g(n+1-x)$$
则原式为
$$\sum_{j=0}^{n+1}x^{j}\frac{1}{j!}\sum_{p+q=n+1+j}f(p)g(q)$$
然后我们大力ntt一波除个阶乘就行了
上代码~
```C
// luogu-judger-enable-o2
#include<cstdio>
#include<algorithm>
using namespace std;const int N=524288+10;typedef unsigned long long ll;const ll mod=998244353;
inline ll po(ll a,ll p){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}
# define md(x) (x=(x>mod)?x-mod:x)
int rv[22][N];ll rt[2][22];ll tr[N];ll ifac[N];ll bu[N];ll fac[N];ll a[N];int n;
ll f1[N];ll ans[N];
inline void pre()//预处理
{
for(int i=1;i<=19;i++)
for(int j=1;j<(1<<i);j++)rv[i][j]=(rv[i][j>>1]>>1)|((j&1)<<(i-1));
for(int i=1,t=(mod-1)/2;i<=20;i++,t>>=1)rt[0][i]=po(3,t);
for(int i=1,t=(mod-1)/2;i<=20;i++,t>>=1)rt[1][i]=po(332748118,t);
}
inline void ntt(ll* a,int len,int d,int o)//简易ntt板子
{
for(int i=0;i<len;i++)if(i<rv[d][i])swap(a[i],a[rv[d][i]]);
for(int k=1,j=1;k<len;k<<=1,j++)
for(int s=0;s<len;s+=(k<<1))
for(int i=s,w=1;i<s+k;i++,w=w*rt[o][j]%mod)
{ll a1=a[i+k]*w%mod;a[i+k]=a[i]+mod-a1;md(a[i+k]);a[i]+=a1;md(a[i]);}
if(o){ll iv=po(len,mod-2);for(int i=0;i<len;i++)(a[i]*=iv)%=mod;}
}
inline void poly_inv(ll* a,ll* b,int len)//简易多项式求逆
{
b[0]=po(a[0],mod-2);
for(int k=2,d=1;k<=len;k<<=1,d++)
{
for(int i=0;i<k;i++)tr[i]=a[i];ntt(tr,(k<<1),d+1,0);ntt(b,(k<<1),d+1,0);
for(int i=0;i<(k<<1);i++)b[i]=b[i]*(2+mod-tr[i]*b[i]%mod)%mod;
ntt(b,(k<<1),d+1,1);for(int i=k;i<(k<<1);i++)b[i]=0;
}
}
int main()
{
pre();ifac[0]=1;
ifac[1]=1;for(int i=2;i<=250005;i++)ifac[i]=(mod-mod/i)*ifac[mod%i]%mod;
for(int i=2;i<=250005;i++)(ifac[i]*=ifac[i-1])%=mod;
poly_inv(ifac+1,bu,262144);//这里懒了直接处理的所有的伯努利数
fac[0]=1;for(int i=1;i<=250005;i++)fac[i]=fac[i-1]*i%mod;
scanf("%d",&n);
for(int i=0;i<=n;i++)scanf("%lld",&a[i]);int len=524288;int d=19;
for(int i=0;i<=n;i++)(a[i]*=fac[i])%=mod;ntt(a,len,d,0);
for(int i=0;i<=n;i++)f1[i]=bu[n-i];ntt(f1,len,d,0);
for(int i=0;i<len;i++)(f1[i]*=a[i])%=mod;ntt(f1,len,d,1);
for(int i=1;i<=n+1;i++)f1[i]=f1[n+(i-1)]%mod;f1[0]=0;//预处理f(n-1)的系数
for(int i=n+2;i<len;i++)f1[i]=0;ntt(f1,len,d,0);
for(int i=0;i<=n+1;i++)ans[i]=ifac[n+1-i];ntt(ans,len,d,0);
for(int i=0;i<len;i++)(ans[i]*=f1[i])%=mod;ntt(ans,len,d,1);//二项式定理展开
for(int i=0;i<=n+1;i++)printf("%lld ",ans[n+1+i]*ifac[i]%mod);return 0;//拜拜程序~
}
```