跳转至

伯努利数

伯努利数 Bn 是一个与数论有密切关联的有理数序列.前几项被发现的伯努利数分别为:

B0=1,B1=−12,B2=16,B3=0,B4=−130,…

等幂求和

伯努利数是由雅各布·伯努利的名字命名的,他在研究 m 次幂和的公式时发现了奇妙的关系.我们记

Sm(n)=∑k=0n−1km=0m+1m+⋯+(n−1)m

伯努利观察了如下一列公式,勾画出一种模式:

S0(n)=nS1(n)=12n2−12nS2(n)=13n3−12n2+16nS3(n)=14n4−12n3+14n2S4(n)=15n5−12n4+13n3−130n

可以发现,在 Sm(n) 中 nm+1 的系数总是 1m+1,nm 的系数总是 −12,nm−1 的系数总是 m12,nm−3 的系数是 −m(m−1)(m−2)720,nm−4 的系数总是零等.

而 nm−k 的系数总是某个常数乘以 mk―,mk― 表示下降阶乘幂,即 m!(m−k)!.

递推公式

Sm(n)=1m+1(B0nm+1+(m+11)B1nm+⋯+(m+1m)Bmn)=1m+1∑k=0m(m+1k)Bknm+1−k

伯努利数由隐含的递推关系定义:

∑j=0m(m+1j)Bj=0,(m>0)B0=1

例如,(20)B0+(21)B1=0,前几个值显然是

n012345678…
Bn1−12160−13001420−130…

证明

利用归纳法证明

这个证明方法来自 Concrete Mathematics 6.5 BERNOULLI NUMBER.

运用二项式系数的恒等变换和归纳法进行证明:

Sm+1(n)+nm+1=∑k=0n−1(k+1)m+1=∑k=0n−1∑j=0m+1(m+1j)kj=∑j=0m+1(m+1j)Sj(n)

令 S^m(n)=1m+1∑k=0m(m+1k)Bknm+1−k,我们希望证明 Sm(n)=S^m(n),假设对 j∈[0,m),有 Sj(n)=S^j(n).

将原式中两边都减去 Sm+1(n) 后可以得到:

Sm+1(n)+nm+1=∑j=0m+1(m+1j)Sj(n)nm+1=∑j=0m(m+1j)Sj(n)=∑j=0m−1(m+1j)S^j(n)+(m+1m)Sm(n)

尝试在式子的右边加上 (m+1m)S^m(n)−(m+1m)S^m(n) 再进行化简,可以得到:

nm+1=∑j=0m(m+1j)S^j(n)+(m+1)(Sm(n)−S^m(n))

不妨设 Δ=Sm(n)−S^m(n),并且将 S^j(n) 展开,那么有

nm+1=∑j=0m(m+1j)S^j(n)+(m+1)Δ=∑j=0m(m+1j)1j+1∑k=0j(j+1k)Bknj+1−k+(m+1)Δ

将第二个 ∑ 中的求和顺序改为逆向,再将组合数的写法恒等变换可以得到:

nm+1=∑j=0m(m+1j)1j+1∑k=0j(j+1j−k)Bj−knk+1+(m+1)Δ=∑j=0m(m+1j)1j+1∑k=0j(j+1k+1)Bj−knk+1+(m+1)Δ=∑j=0m(m+1j)1j+1∑k=0jj+1k+1(jk)Bj−knk+1+(m+1)Δ=∑j=0m(m+1j)∑k=0j(jk)Bj−kk+1nk+1+(m+1)Δ

对两个求和符号进行交换,可以得到:

nm+1=∑k=0mnk+1k+1∑j=km(m+1j)(jk)Bj−k+(m+1)Δ

对 (m+1j)(jk) 进行恒等变换:

(m+1j)(jk)=(m+1k)(m−k+1j−k)

那么式子就变成了:

nm+1=∑k=0mnk+1k+1∑j=km(m+1k)(m−k+1j−k)Bj−k+(m+1)Δ=∑k=0mnk+1k+1(m+1k)∑j=km(m−k+1j−k)Bj−k+(m+1)Δ

将所有的 j−k 用 j 代替,那么就可以得到:

nm+1=∑k=0mnk+1k+1(m+1k)∑j=0m−k(m−k+1j)Bj+(m+1)Δ

考虑我们前面提到过的递归关系

∑j=0m(m+1j)Bj=0,(m>0)B0=1∑j=0m(m+1j)Bj=[m=0]

代入后可以得到:

nm+1=∑k=0mnk+1k+1(m+1k)[m−k=0]+(m+1)Δ=nm+1m+1(m+1m)+(m+1)Δ=nm+1+(m+1)Δ

于是 Δ=0,且有 Sm(n)=S^m(n).

利用指数生成函数证明

对递推式 ∑j=0m(m+1j)Bj=[m=0]

两边都加上 Bm+1,即得到:

∑j=0m+1(m+1j)Bj=[m=0]+Bm+1∑j=0m(mj)Bj=[m=1]+Bm∑j=0mBjj!⋅1(m−j)!=[m=1]+Bmm!

设 B(z)=∑i≥0Bii!zi,注意到左边为卷积形式,故:

B(z)ez=z+B(z)B(z)=zez−1

设 Fn(z)=∑m≥0Sm(n)m!zm,则:

Fn(z)=∑m≥0Sm(n)m!zm=∑m≥0∑i=0n−1imzmm!

调换求和顺序:

Fn(z)=∑i=0n−1∑m≥0imzmm!=∑i=0n−1eiz=enz−1ez−1=zez−1⋅enz−1z

代入 B(z)=zez−1:

Fn(z)=B(z)⋅enz−1z=(∑i≥0Bii!)(∑i≥1nizi−1i!)=(∑i≥0Bii!)(∑i≥0ni+1zi(i+1)!)

由于 Fn(z)=∑m≥0Sm(n)m!zm,即 Sm(n)=m![zm]Fn(z):

S×m(n)=m![zm]Fn(z)=m!∑i=0mB×ii!⋅nm−i+1(m−i+1)!=1m+1∑i=0m(m+1i)Binm−i+1

故得证.

参考实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
using ll = long long;
constexpr int MAXN = 10000;
constexpr int mod = 1e9 + 7;
ll B[MAXN];        // 伯努利数
ll C[MAXN][MAXN];  // 组合数
ll inv[MAXN];      // 逆元(计算伯努利数)

void init() {
  // 预处理组合数
  for (int i = 0; i < MAXN; i++) {
    C[i][0] = C[i][i] = 1;
    for (int k = 1; k < i; k++) {
      C[i][k] = (C[i - 1][k] % mod + C[i - 1][k - 1] % mod) % mod;
    }
  }
  // 预处理逆元
  inv[1] = 1;
  for (int i = 2; i < MAXN; i++) {
    inv[i] = (mod - mod / i) * inv[mod % i] % mod;
  }
  // 预处理伯努利数
  B[0] = 1;
  for (int i = 1; i < MAXN; i++) {
    ll ans = 0;
    if (i == MAXN - 1) break;
    for (int k = 0; k < i; k++) {
      ans += C[i + 1][k] * B[k];
      ans %= mod;
    }
    ans = (ans * (-inv[i + 1]) % mod + mod) % mod;
    B[i] = ans;
  }
}