伯努利数学习笔记&&Luogu P3711 仓鼠的数学题

新科技

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=;char zf=;char ch=getchar();
while(ch!='-'&&!isdigit(ch))ch=getchar();
if(ch=='-')zf=-,ch=getchar();
while(isdigit(ch))x=x*+ch-'',ch=getchar();return x*zf;
}
void write(ll y){if(y<)putchar('-'),y=-y;if(y>)write(y/);putchar(y%+);}
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-){
int ans=;
for(;y;y>>=,x=1ll*x*x%p)if(y&)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=,j=;i<n;i++){
F[i]=A[j];
for(rt k=n>>;(j^=k)<k;k>>=);
}
for(rt i=;i<n;i<<=){
const int w=W[]=ksm(,(p-)//i);W[]=;
for(rt k=;k<i;k++)W[k]=1ll*W[k-]*w%p;
for(rt j=;j<n;j+=i<<){
for(rt k=;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=;i<n;i++)A[i]=F[i]%p;
if(fla==-){
const int invn=ksm(n);
reverse(A.begin()+,A.end());
for(rt i=;i<n;i++)A[i]=1ll*A[i]*invn%p;
}
}
poly Mul(poly x,poly y){
int sz=x.size()+y.size()-,lim=;
while(lim<=sz)lim<<=;
NTT(lim,x,);NTT(lim,y,);
for(rt i=;i<lim;i++)x[i]=1ll*x[i]*y[i]%p;
NTT(lim,x,-);x.resize(sz);return x;
}
poly Inv(poly a,int n=-){
if(n==-)n=a.size();
if(n==)return {ksm(a[])};
poly c=Inv(a,n+>>),d(&a[],&a[n]);
int lim=;while(lim<=n*)lim<<=;
NTT(lim,c,);NTT(lim,d,);
for(rt i=;i<lim;i++)c[i]=1ll*c[i]*(2ll+p-1ll*d[i]*c[i]%p)%p;
NTT(lim,c,-);c.resize(n);return c;
}
}
using namespace Poly;
int inv[],jc[],njc[],a[];
poly B;
void init(int k){
for(rt i=;i<=;i++)inv[i]=jc[i]=njc[i]=;
for(rt i=;i<=k+;i++){
inv[i]=1ll*inv[p%i]*(p-p/i)%p;
jc[i]=1ll*jc[i-]*i%p;
njc[i]=1ll*njc[i-]*inv[i]%p;
}
B.resize(k+);
for(rt i=;i<=k;i++)B[i]=njc[i+];
B=Inv(B);
for(rt i=;i<=k;i++)B[i]=1ll*B[i]*jc[i]%p;
}
int main(){
n=read();init(n+);
for(rt i=;i<=n;i++)a[i]=read();
poly ans(n+),C(n+);
for(rt i=;i<=n;i++)B[i]=1ll*B[i]*njc[i]%p;
for(rt i=;i<=n;i++)C[i]=1ll*jc[i]*a[i]%p;
reverse(&B[],&B[n+]);B.resize(n+);C.resize(n+);
ans=Mul(B,C);
for(rt i=;i<=n+;i++)ans[n+i-]=1ll*ans[n+i-]*njc[i]%p;
for(rt i=;i<=n;i++)(ans[n+i-]+=a[i])%=p;write(a[]),putchar(' ');
for(rt i=;i<=n+;i++)write(ans[n+i-]),putchar(' ');
return ;
}
上一篇:linux中历史命令的一点发现


下一篇:HTML页面表单输入框去掉鼠标选中后边框变色的效果