伯努利数 仓鼠的数学题

题面
设 $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\left(\genfrac{}{}{0ex}{}{p+1}{j}\right)B_j{n}^{p+1-j}$
代回原式：
${\sum }_{i=0}^n{i}^j=[j=0]+\frac{1}{j+1}{\sum }_{k=0}^j\left(\genfrac{}{}{0ex}{}{p+1}{j}\right)B_j{n}^{p+1-j}$
因此伯努利数可以用于求自然数幂和。
回归原题：
$c_i={\sum }_{j=i-1}^n\frac{a_j}{j+1}\left(\genfrac{}{}{0ex}{}{j+1}{j+1-i}\right)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 }_{2}n)$，空间复杂度 $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;
}