伯努利数

Bernoulli number

一些性质

Wikipedia 说，

\[ \begin{aligned} B^-_0=1 \\ \sum _{j=0}^{m}\binom{m+1}{j}B^-_{j}=0 \end{aligned} \]

通过上式可以计算伯努利数，计算得到的是第一伯努利数，记做 \(B^-_n\)，与第二伯努利数 \(B^+_n\) 仅有的区别是

\[ \begin{aligned} B^-_1=-\frac{1}{2} \\ B^+_1=+\frac{1}{2} \end{aligned} \]

通过转化可以得到

\[ {\frac {x}{e^{x}-1}}=\sum _{n=0}^{\infty }B^-_{n}{\frac {x^{n}}{n!}} \]

可以用多项式求逆在 \(O(n\log n)\) 的复杂度内求出伯努利数的前 \(n\) 项。

第二伯努利数满足

\[ \begin{aligned} S_{m}(n)&=\sum _{k=1}^{n}k^{m}=1^{m}+2^{m}+\cdots +{n}^{m} \\ &= \frac {1}{m+1}\sum _{k=0}^{m}\binom{m+1}{k} B_{k}^{+}n^{m+1-k} \end{aligned} \]

可以用来求自然数幂和。

例题

Luogu P3711 仓鼠的数学题

题意

设关于 \(x\) 的 \(k+1\) 次多项式 \(S_k(x)=\sum\limits_{i=0}^x i^k\)。

给出 \(n+1\) 个数 \(a_0,a_1,\dotsc,a_n\)，求 \(\sum\limits_{k=0}^n a_k S_k(x)\) 的各项系数。

做法

把上面的式子代入有

\[ \begin{aligned} & \sum_{k=1}^n a_kS_{k}(x) \\ = & \sum_{k=1}^n \frac {a_k}{k+1}\sum _{i=0}^{k}\binom{k+1}{i} B_{i}^{+}x^{k+1-i} \\ = & \sum_{k=1}^n \frac {a_k}{k+1}\sum _{i=1}^{k+1}\binom{k+1}{i} B_{k+1-i}^{+}x^i \\ = & \sum_{k=1}^n a_k k! \sum _{i=1}^{k+1} \frac {B_{k+1-i}^{+}x^{i}}{i!(k+1-i)!} \\ = & \sum _{i=1}^{n+1} \frac{x^i}{i!} \sum_{k=i-1}^n a_k k! \frac {B_{k+1-i}^{+}}{(k+1-i)!} \\ \end{aligned} \]

右边的部分可以求出伯努利数后卷积，总复杂度 \(O(n\log n)\)。

注意通过

\[ {\frac {x}{e^{x}-1}}=\sum _{n=0}^{\infty }B^-_{n}{\frac {x^{n}}{n!}} \]

求出来的是 \(B^-_n\)，后面计算需要改动 \(B_1\) 的值。

代码

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<string.h>
#include<cmath>
#include<vector>

using namespace std;
#define ll long long

inline char read() {
	static const int IN_LEN = 1000000;
	static char buf[IN_LEN], *s, *t;
	return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
	static bool iosig;
	static char c;
	for (iosig=false, c=read(); !isdigit(c); c=read()) {
		if (c == '-') iosig=true;
		if (c == -1) return;
	}
	for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
	if (iosig) x=-x;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
	if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
	*ooh++=c;
}
template<class T>
inline void print(T x) {
	static int buf[30], cnt;
	if (x==0) print('0');
	else {
		if (x<0) print('-'), x=-x;
		for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
		while(cnt) print((char)buf[cnt--]);
	}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

const int N = 250005, M = 1<<19, P = 998244353;
int n, fac[N], ifac[N], w[M];
unsigned ll F[M];
inline int Pow(ll x, int y=P-2){
	int ans=1;
	for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
	return ans;
}
inline int Calc(int x){ int ans=1; while(ans<=x) ans<<=1; return ans;}
inline void DFT(vector<int> &f, int n){
	for(int i=0, j=0; i<n; ++i){
		F[i]=f[j];
		for(int k=n>>1; (j^=k)<k; k>>=1);
	}
	for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1) for(int k=j; k<j+i; ++k){
		int t=F[k+i]*w[i+k-j]%P;
		F[k+i]=F[k]+P-t, F[k]+=t;
	}
	for(int i=0; i<n; ++i) f[i]=F[i]%P;
}
inline void IDFT(vector<int> &f, int n){
	reverse(f.begin()+1, f.end()), DFT(f, n);
	for(int i=0, I=Pow(n); i<n; ++i) f[i]=(ll)f[i]*I%P;
}
inline vector<int> operator *(const vector<int> &a, const vector<int> &b){
	int n=Calc(a.size()+b.size()-2);
	vector<int> A=a, B=b;
	A.resize(n), B.resize(n), DFT(A, n), DFT(B, n);
	for(int i=0; i<n; ++i) A[i]=(ll)A[i]*B[i]%P;
	IDFT(A, n);
	return A.resize(a.size()+b.size()-1), A;
}
vector<int> PolyInv(const vector<int> &a, int n=-1){
	if(n==-1) n=a.size();
	if(n==1) return {Pow(a[0])};
	vector<int> ans=PolyInv(a, (n+1)>>1), A(a.begin(), a.begin()+n);
	int m=Calc(n*2-1);
	ans.resize(m), A.resize(m), DFT(ans, m), DFT(A, m);
	for(int i=0; i<m; ++i) ans[i]=(2+(ll)(P-A[i])*ans[i])%P*ans[i]%P;
	IDFT(ans, m);
	return ans.resize(n), ans;
}
int main() {
	for(int i=1; i<M; i<<=1){
		w[i]=1, w[i+1]=Pow(3, (P-1)/i/2);
		for(int j=2; j<i; ++j) w[i+j]=(ll)w[i+j-1]*w[i+1]%P;
	}
	read(n);
	vector<int> a(n+1), b(n+1);
	fac[0]=1;
	for(int i=1; i<=n+1; ++i) fac[i]=(ll)fac[i-1]*i%P;
	ifac[n+1]=Pow(fac[n+1]);
	for(int i=n+1; i; --i) ifac[i-1]=(ll)ifac[i]*i%P;
	for(int i=0; i<=n; ++i) read(a[i]), a[i]=(ll)a[i]*fac[i]%P;
	print(a[0]);
	for(int i=0; i<=n; ++i) b[i]=ifac[i+1];
	b=PolyInv(b), b[1]=P-b[1];
	reverse(b.begin(), b.end());
	a=a*b;
	for(int i=1; i<=n+1; ++i) print(' '), print((int)((ll)a[i+n-1]*ifac[i]%P));
	return flush(), 0;
}