伯努利数仓鼠的数学题

原创于 2021-03-30 07:24:38 发布 · 241 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#数学

数学专栏收录该内容

15 篇文章

订阅专栏

题面
设 $B_i$ 为伯努利数的第 $i$ 项。伯努利数有许多定义方法，这里定义 $\sum_{i=0}\frac{B_i}{i!}x^i=\frac{x}{1-e^{-x}}$ 。设 $G_n(x)=\sum_{i=0}\frac{S_i(n)}{n!}x^i$ ，带入 $S_i(n)$ 得到：
$G_n(x)=\sum_{i=0}\sum_{j=0}^n\frac{(jx)^i}{i!}$
利用 $e^x$ 的泰勒展开：
$G_n(x)=\sum_{i=0}^n e^{ix}$
利用等比数列求和公式：
$G_n(x)=\sum_{i=0}^n\frac{e^{-x}-e^{nx}}{e^{-x}-1}$
根据伯努利数的定义：
$G_n(x)=1-\sum_{i=0}\frac{B_i}{i!}x^{i-1}(1-e^{(n+1)x})$
再次泰勒展开：
$G_n(x)=1+\sum_{i=0}\frac{B_i}{i!}x^{i-1}\sum_{j=1}\frac{(nx)^j}{j!}$
做多项式乘法：
$G_n(x)=1+\sum_{i=0}\sum_{j=0}^i \frac{B_j n^{p+1-j}}{j!(i+1-j)!} x^i$
阶乘转化为组合数：
$G_n(x)=1+\sum_{i=0}\frac{x^i}{i!}\frac{1}{p+1}\sum_{j=0}^i \binom{p+1}{j}B_j n^{p+1-j}$
代回原式：
$\sum_{i=0}^n i^j=[j=0]+\frac{1}{j+1}\sum_{k=0}^j \binom{p+1}{j}B_j n^{p+1-j}$
因此伯努利数可以用于求自然数幂和。
回归原题：
$c_i=\sum_{j=i-1}^n\frac{a_j}{j+1}\binom{j+1}{j+1-i}B_{j+1-i}$
拆掉组合数：
$c_i i!=\sum_{j=i-1}^n\frac{a_j j! B_{j+1-i}}{(j+1-i)!}$
显然可以多项式卷积求出 $c$ ，其中 $c_0$ 需要特殊处理，即 $c_0=a_0$ 。
根据伯努利数的生成函数，可以利用多项式求逆求出。
时间复杂度 $O(n \log_2n)$ ，空间复杂度 $O (n)$

#include<stdio.h>
#define R register int
#define L long long
#define I inline
#define N 524288
#define M 250003
#define P 998244353
int a[N],b[N],c[N],fac[M],invf[M];
I void Swap(int&x,int&y){
	int tem=x;
	x=y;
	y=tem;
}
I int Add(int x,int y){
	x+=y;
	if(x<0){
		return x+P;
	}
	return x>=P?x-P:x;
}
I int PowMod(int x,int y){
	int s=1;
	while(y!=0){
		if((y&1)==1){
			s=(L)s*x%P;
		}
		x=(L)x*x%P;
		y>>=1;
	}
	return s;
}
I void NTT(int*A,const int len,const short type){
	int tem=0;
	for(R i=0;i!=len;i++){
		if(i<tem){
			Swap(A[i],A[tem]);
		}
		R j=len;
		do{
			j>>=1;
			tem^=j;
		}while(tem<j);
	}
	static int w[N];
	w[0]=1;
	for(R i=1;i!=len;i<<=1){
		tem=i<<1;
		int omg=PowMod(3,P-1+(P-1)*type/tem);
		for(R j=1;j!=i;j++){
			w[j]=(L)w[j-1]*omg%P;
		}
		for(R j=0;j!=len;j+=tem){
			for(R k=j;k!=i+j;k++){
				int t1=A[k],t2=(L)A[i+k]*w[k-j]%P;
				A[k]=Add(t1,t2);
				A[i+k]=Add(t1,-t2);
			}
		}
	}
	if(type==-1){
		tem=PowMod(len,P-2);
		for(R i=0;i!=len;i++){
			A[i]=(L)A[i]*tem%P;
		}
	}
}
I void PolyInverse(int*A,int*B,const int len){
	static int tem[N];
	B[0]=PowMod(A[0],P-2);
	for(R i=1;i!=len;i<<=1){
		int p=i<<2;
		NTT(B,p,1);
		for(R j=0;j!=i<<1;j++){
			tem[j]=A[j];
		}
		NTT(tem,p,1);
		for(R j=0;j!=p;j++){
			B[j]=(L)B[j]*(P+2-(L)tem[j]*B[j]%P)%P;
		}
		NTT(B,p,-1);
		for(R j=i<<1;j!=p;j++){
			B[j]=tem[j]=0;
		}
	}
	for(R i=0;i!=len<<1;i++){
		tem[i]=0;
	}
}
int main(){
	fac[0]=1;
	int n,x,len=1;
	scanf("%d",&n);
	for(R i=1;i!=n+3;i++){
		fac[i]=(L)fac[i-1]*i%P;
	}
	for(R i=0;i<=n;i++){
		scanf("%d",&x);
		a[n-i]=(L)x*fac[i]%P;
	}
	printf("%d",a[n]);
	n+=2;
	invf[n]=PowMod(fac[n],P-2);
	for(R i=n;i!=0;i--){
		invf[i-1]=(L)invf[i]*i%P;
	}
	while(len<n){
		len<<=1;
	}
	for(R i=0;i<=n;i++){
		c[i]=invf[i+1];
		if((i&1)==1){
			c[i]=P-c[i];
		}
	}
	PolyInverse(c,b,len);
	len<<=1;
	NTT(a,len,1);
	NTT(b,len,1);
	for(R i=0;i!=len;i++){
		a[i]=(L)a[i]*b[i]%P;
	}
	NTT(a,len,-1);
	for(R i=1;i!=n;i++){
		printf(" %d",(L)invf[i]*a[n-i-1]%P);
	}
	return 0;
}