伯努利（Bernoulli）数学习笔记

伯努利数的指数型生成函数为$B(x)={\sum }_{i=0}\frac{B_i}{i!}{x}^i=\frac{x}{e^x-1}$

由此可得$B_0=1,B_1=-\frac{1}{2},B_3=\frac{1}{6},B_4=0,B_5=\frac{1}{30}...$

如何求伯努利数？显然可以多项式求逆，或者通过一个性质$O(n^2)$求。

（以下$[x^n]B(x)$指的是$B(x)$的$n$次项系数。）

因为$\frac{x}{e^x-1}(e^x-1)=x$

所以$[x^n]B(x)(e^x-1)={\sum }_{i=0}^{n-1}\frac{B_i}{i!}\ast \frac{1}{(n-i)!}=[n=1]$

所以${\sum }_{i=0}^{n-1}{B}_i{C}_n^i=[n=1]n!=[n=1]$

那么已知$B_0=1$，${\sum }_{i=0}^n{B}_i{C}_{n+1}^i=0$，则$B_n=-\frac{{\sum }_{i=0}^{n-1}{B}_i{C}_{n+1}^i}{n+1}$

于是就可以$O(n^2)$方便地求伯努利数，或者$O(n\log n)$常数巨大又不方便地求。

伯努利数还有个重要性质。

记$S_k(n)={\sum }_{i=0}^{n-1}{i}^k$

则$S_k(n)=\frac{1}{k+1}{\sum }_{i=0}^k{C}_{k+1}^i{B}_i{n}^{k+1-i}$

模板题1 51nod 1228 $O(n^2)$求伯努利数

#include<bits/stdc++.h>
using namespace std;
#define RI register int
typedef long long LL;
const int mod=1000000007;
int B[2005],C[2005][2005],mi[2005],K,T;LL n;

int qm(int x) {return x>=mod?x-mod:x;}
int ksm(int x,int y) {
	int re=1;
	for(;y;y>>=1,x=1LL*x*x%mod) if(y&1) re=1LL*re*x%mod;
	return re;
}
void prework() {
	for(RI i=0;i<=2001;++i) {
		C[i][0]=1;
		for(RI j=1;j<=i;++j) C[i][j]=qm(C[i-1][j]+C[i-1][j-1]);
	}
	B[0]=1;
	for(RI i=1;i<=2000;++i) {
		for(RI j=0;j<i;++j) B[i]=qm(B[i]+1LL*B[j]*C[i+1][j]%mod);
		B[i]=qm(mod-1LL*B[i]*ksm(i+1,mod-2)%mod);
	}
}
int main()
{
	prework();
	scanf("%d",&T);
	while(T--) {
		scanf("%lld%d",&n,&K);
		n=qm(n%mod+1);int ans=0;
		mi[0]=1;for(RI i=1;i<=K+1;++i) mi[i]=1LL*mi[i-1]*n%mod;
		for(RI i=0;i<=K;++i) ans=qm(ans+1LL*C[K+1][i]*B[i]%mod*mi[K+1-i]%mod);
		printf("%lld\n",1LL*ans*ksm(K+1,mod-2)%mod);
	}
    return 0;
}

模板题2 51nod 1258 $O(n\log n)$求伯努利数

#include<bits/stdc++.h>
using namespace std;
#define RI register int
typedef long long LL;
typedef long double db;
const int mod=1000000007,N=131080,M=32767;
const db pi=acos(-1);
int K,T;LL n;
int fac[N],inv[N],ifac[N],rev[N],len[N];
int aa[N],bb[N],k1[N],k2[N];
struct com{db r,i;}a[N],b[N],Aa[N],Ab[N],Ba[N],Bb[N];
com operator + (com A,com B) {return (com){A.r+B.r,A.i+B.i};}
com operator - (com A,com B) {return (com){A.r-B.r,A.i-B.i};}
com operator * (com A,com B) {return (com){A.r*B.r-A.i*B.i,A.r*B.i+A.i*B.r};}
com conj(com A) {return (com){A.r,-A.i};}

int qm(int x) {return x>=mod?x-mod:x;}
int ksm(int x,int y) {
	int re=1;
	for(;y;y>>=1,x=1LL*x*x%mod) if(y&1) re=1LL*re*x%mod;
	return re;
}
void FFT(com *a,int n) {
	for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
	for(RI i=1;i<n;i<<=1) {
		com wn=(com){cos(pi/i),sin(pi/i)};
		for(RI j=0;j<n;j+=(i<<1)) {
			com t1,t2,w=(com){1,0};
			for(RI k=0;k<i;++k,w=w*wn)
				t1=a[j+k],t2=a[j+i+k]*w,a[j+k]=t1+t2,a[j+i+k]=t1-t2;
		}
	}
}
void mul(int *ka,int *kb,int *kc,int n) {
	for(RI i=0;i<n;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len[n]-1));
	for(RI i=0;i<(n>>1);++i) {
		a[i]=(com){(db)(ka[i]&M),(db)(ka[i]>>15)};
		b[i]=(com){(db)(kb[i]&M),(db)(kb[i]>>15)};
		a[i+(n>>1)]=b[i+(n>>1)]=(com){0,0};
	}
	FFT(a,n),FFT(b,n);
	for(RI i=0;i<n;++i) {
		int j=(n-i)&(n-1);
		com kAa=(a[i]+conj(a[j]))*(com){0.5,0};
		com kAb=(a[i]-conj(a[j]))*(com){0,-0.5};
		com kBa=(b[i]+conj(b[j]))*(com){0.5,0};
		com kBb=(b[i]-conj(b[j]))*(com){0,-0.5};
		Aa[j]=kAa*kBa,Ab[j]=kAa*kBb,Ba[j]=kAb*kBa,Bb[j]=kAb*kBb;
	}
	for(RI i=0;i<n;++i)
		a[i]=Aa[i]+Ab[i]*(com){0,1},b[i]=Ba[i]+Bb[i]*(com){0,1};
	FFT(a,n),FFT(b,n);
	for(RI i=0;i<n;++i) {
		int kAa=(LL)(a[i].r/n+0.5)%mod;
		int kAb=(LL)(a[i].i/n+0.5)%mod;
		int kBa=(LL)(b[i].r/n+0.5)%mod;
		int kBb=(LL)(b[i].i/n+0.5)%mod;
		kc[i]=qm(((LL)kAa+((LL)(kAb+kBa)<<15)+((LL)kBb<<30))%mod+mod);
	}
}
void getinv(int *a,int *b,int n) {
	if(n==1) {b[0]=ksm(a[0],mod-2),b[1]=0;return;}
	getinv(a,b,n>>1),mul(a,b,k1,n<<1),mul(k1,b,k2,n<<1);
	for(RI i=0;i<n;++i) b[i]=qm(qm(b[i]+b[i])-k2[i]+mod),b[i+n]=0;
}

int C(int d,int u) {return 1LL*fac[d]*ifac[u]%mod*ifac[d-u]%mod;}
void prework() {
	inv[0]=inv[1]=1,ifac[0]=1,fac[0]=1;
	for(RI i=1;i<N;++i) fac[i]=1LL*fac[i-1]*i%mod;
	for(RI i=2;i<N;++i) inv[i]=1LL*(mod-mod/i)*inv[mod%i]%mod;
	for(RI i=1;i<N;++i) ifac[i]=1LL*ifac[i-1]*inv[i]%mod;
	int kn=1;while(kn<=50000) kn<<=1,len[kn]=len[kn>>1]+1;
	len[kn<<1]=len[kn]+1;
	for(RI i=0;i<kn;++i) aa[i]=ifac[i+1];
	getinv(aa,bb,kn);
	for(RI i=0;i<=50000;++i) bb[i]=1LL*bb[i]*fac[i]%mod;
}
int main()
{
	prework();
	scanf("%d",&T);
	while(T--) {
		scanf("%lld%d",&n,&K);
		n=qm(n%mod+1);int ans=0;
		for(RI i=1,j=n;i<=K+1;++i,j=1LL*j*n%mod)
			ans=qm(ans+1LL*C(K+1,K+1-i)*bb[K+1-i]%mod*j%mod);
		printf("%lld\n",1LL*ans*inv[K+1]%mod);
	}
    return 0;
}