伯努利数学习笔记&&Luogu P3711 仓鼠的数学题
新科技
题意
设$ S_{k,n}$表示$ \displaystyle\sum_{i=0}^n i^k$
求多项式$\displaystyle\sum_{k=0}^n S_{k,x}a_k$的各项系数
数组$ a$给定,$ n \leq 100000$
伯努利数
伯努利数$B$是一个数列，满足
$$\sum_{i=0}^n B_i\binom{n+1}{i}=0$$
可以用它来求自然数幂和
$$ S_{k,n-1}=\sum_{i=0}^{n-1}i^k=\frac{1}{k+1}\sum_{i=0}^k\binom{k+1}{i}B_in^{k+1-i}$$
如果已经得到了数列$ B$，求自然数幂和$S_{k,n}$是$ O(k)$的
直接根据定义可以$ O(n^2)$递推伯努利数，考虑更快速的推法
$$
\begin{aligned}
\sum_{i=0}^n B_i\binom{n+1}{i}&=0\\
\sum_{i=0}^{n-1} B_i\binom{n}{i}&=0 \ (n>1)\\
B_n+\sum_{i=0}^{n-1} B_i\binom{n}{i}&=B_n \ (n>1)\\
B_n&=\sum_{i=0}^nB_i\binom{n}{i} \ (n>1)\\
\frac{B_n}{n!}&=\sum_{i=0}^n\frac{B_i}{i!(n-i)!}\\
\end{aligned}
$$
设伯努利数的指数型生成函数为$ B$，伯努利数的第一项$ B_1=-\frac{1}{2}$
则有$B*e^x=B+x$
整理得$B=\frac{x}{e^x-1}=(\frac{e^x-1}{x})^{-1}$
直接多项式求逆即可
时间复杂度$ O(n \log n)$
回到原题
用伯努利数展开得
$$
\begin{aligned}
ans&=\sum_{k=0}^na_k S_{k,x}\\
&=\sum_{k=0}^na_k(x^k+\frac{1}{k+1}\sum_{i=0}^k\binom{k+1}{i}B_ix^{k+1-i})\\
&=(\sum_{k=0}^na_kx^k)+(\sum_{k=0}^nk!\sum_{i=0}^k\frac{B_i}{i!(k+1-i)!}x^{k+1-i})\\
ans[x^d]&=a_d+\sum_{i=0}^{n+1}\frac{B_i}{d!i!}(d+i-1)!\\
\frac{ans[x^d]}{d!}&=a_d+\sum_{i=0}^{n+1}\frac{B_i}{i!}(d+i-1)!
\end{aligned}
$$
发现这是一个差卷积的形式
按套路反转之后$ NTT$即可
总复杂度仍是$ O(n \log n)$
代码
#include<ctime> #include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #include<queue> #include<vector> #define p 998244353 #define rt register int #define ll long long #define ull unsigned long long using namespace std; inline ll read(){ ll x=0;char zf=1;char ch=getchar(); while(ch!='-'&&!isdigit(ch))ch=getchar(); if(ch=='-')zf=-1,ch=getchar(); while(isdigit(ch))x=x*10+ch-'0',ch=getchar();return x*zf; } void write(ll y){if(y<0)putchar('-'),y=-y;if(y>9)write(y/10);putchar(y%10+48);} void writeln(const ll y){write(y);putchar('\n');} int k,m,n,x,y,z,cnt,ans; namespace Poly{ #define poly vector<int> #define MAXN 524288 int ksm(int x,int y=p-2){ int ans=1; for(;y;y>>=1,x=1ll*x*x%p)if(y&1)ans=1ll*ans*x%p; return ans; } void NTT(int n,poly &A,int fla){ static ull F[MAXN],W[MAXN];A.resize(n); for(rt i=0,j=0;i<n;i++){ F[i]=A[j]; for(rt k=n>>1;(j^=k)<k;k>>=1); } for(rt i=1;i<n;i<<=1){ const int w=W[1]=ksm(3,(p-1)/2/i);W[0]=1; for(rt k=2;k<i;k++)W[k]=1ll*W[k-1]*w%p; for(rt j=0;j<n;j+=i<<1){ for(rt k=0;k<i;k++){ const ull x=F[j+k],y=F[i+j+k]*W[k]%p; F[j+k]=x+y,F[i+j+k]=x+p-y; } } } for(rt i=0;i<n;i++)A[i]=F[i]%p; if(fla==-1){ const int invn=ksm(n); reverse(A.begin()+1,A.end()); for(rt i=0;i<n;i++)A[i]=1ll*A[i]*invn%p; } } poly Mul(poly x,poly y){ int sz=x.size()+y.size()-1,lim=1; while(lim<=sz)lim<<=1; NTT(lim,x,1);NTT(lim,y,1); for(rt i=0;i<lim;i++)x[i]=1ll*x[i]*y[i]%p; NTT(lim,x,-1);x.resize(sz);return x; } poly Inv(poly a,int n=-1){ if(n==-1)n=a.size(); if(n==1)return {ksm(a[0])}; poly c=Inv(a,n+1>>1),d(&a[0],&a[n]); int lim=1;while(lim<=n*2)lim<<=1; NTT(lim,c,1);NTT(lim,d,1); for(rt i=0;i<lim;i++)c[i]=1ll*c[i]*(2ll+p-1ll*d[i]*c[i]%p)%p; NTT(lim,c,-1);c.resize(n);return c; } } using namespace Poly; int inv[250010],jc[250010],njc[250010],a[250010]; poly B; void init(int k){ for(rt i=0;i<=1;i++)inv[i]=jc[i]=njc[i]=1; for(rt i=2;i<=k+2;i++){ inv[i]=1ll*inv[p%i]*(p-p/i)%p; jc[i]=1ll*jc[i-1]*i%p; njc[i]=1ll*njc[i-1]*inv[i]%p; } B.resize(k+1); for(rt i=0;i<=k;i++)B[i]=njc[i+1]; B=Inv(B); for(rt i=0;i<=k;i++)B[i]=1ll*B[i]*jc[i]%p; } int main(){ n=read();init(n+2); for(rt i=0;i<=n;i++)a[i]=read(); poly ans(n+2),C(n+1); for(rt i=0;i<=n;i++)B[i]=1ll*B[i]*njc[i]%p; for(rt i=0;i<=n;i++)C[i]=1ll*jc[i]*a[i]%p; reverse(&B[0],&B[n+1]);B.resize(n+1);C.resize(n+1); ans=Mul(B,C); for(rt i=0;i<=n+1;i++)ans[n+i-1]=1ll*ans[n+i-1]*njc[i]%p; for(rt i=0;i<=n;i++)(ans[n+i-1]+=a[i])%=p;write(a[0]),putchar(' '); for(rt i=1;i<=n+1;i++)write(ans[n+i-1]),putchar(' '); return 0; }
