题解 P3711 【仓鼠的数学题】

shadowice1984

2018-12-12 09:19:23

Solution

不知道为什么这题把自然数幂和函数的求和上限额外加了一个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;//拜拜程序~ } ```