Bernoulli Numbers

Bernoulli numbers BnB_n form a sequence of rational numbers closely related to number theory. The first few Bernoulli numbers are as follows:

B0=1,B1=−12,B2=16,B3=0,B4=−130,…B_0=1,B_1=-\frac{1}{2},B_2=\frac{1}{6},B_3=0,B_4=-\frac{1}{30},\dots

Sums of like powers

Bernoulli numbers are named after Jacob Bernoulli, who discovered a remarkable relationship while studying formulas for sums of mm-th powers. Let

Sm(n)=∑k=0n−1km=0m+1m+⋯+(n−1)mS_{m}(n)=\sum_{k=0}^{n-1}k^m=0^m+1^m+\dots+(n-1)^m

Bernoulli observed the following sequence of formulas and sketched a pattern:

S0(n)=nS1(n)=12n2−12nS2(n)=13n3−12n2+16nS3(n)=14n4−12n3+14n2S4(n)=15n5−12n4+13n3−130n\begin{aligned} S_0(n)&=n\\ S_1(n)&=\frac{1}{2}n^2-\frac{1}{2}n\\ S_2(n)&=\frac{1}{3}n^3-\frac{1}{2}n^2+\frac{1}{6}n\\ S_3(n)&=\frac{1}{4}n^4-\frac{1}{2}n^3+\frac{1}{4}n^2\\ S_4(n)&=\frac{1}{5}n^5-\frac{1}{2}n^4+\frac{1}{3}n^3-\frac{1}{30}n \end{aligned}

In Sm(n)S_m(n), the coefficient of nm+1n^{m+1} is always 1m+1\frac{1}{m+1}, the coefficient of nmn^m is always −12-\frac{1}{2}, the coefficient of nm−1n^{m-1} is always m12\frac{m}{12}, the coefficient of nm−3n^{m-3} is −m(m−1)(m−2)720-\frac{m(m-1)(m-2)}{720}, and the coefficient of nm−4n^{m-4} is always zero, among others.

The coefficient of nm−kn^{m-k} is always a constant multiplied by mk‾m^{\underline{k}}, where mk‾m^{\underline{k}} denotes the falling factorial, namely m!(m−k)!\frac{m!}{(m-k)!}.

Recurrence formula

Sm(n)=1m+1(B0nm+1+(m+11)B1nm+⋯+(m+1m)Bmn)=1m+1∑k=0m(m+1k)Bknm+1−k\begin{aligned} S_m{(n)}&=\frac{1}{m+1}(B_0n^{m+1}+\binom{m+1}{1}B_1 n^m+\dots+\binom{m+1}{m}B_m n) \\ &=\frac{1}{m+1}\sum_{k=0}^{m}\binom{m+1}{k}B_kn^{m+1-k} \end{aligned}

Bernoulli numbers are defined by the implicit recurrence:

∑j=0m(m+1j)Bj=0,(m>0)B0=1\begin{aligned} \sum_{j=0}^{m}\binom{m+1}{j}B_j&=0,(m>0)\\ B_0&=1 \end{aligned}

For example, (20)B0+(21)B1=0\binom{2}{0}B_0+\binom{2}{1}B_1=0, and the first few values are clearly

nn001122334455667788…\dots
BnB_n11−12-\frac{1}{2}16\frac{1}{6}00−130-\frac{1}{30}00142\frac{1}{42}00−130-\frac{1}{30}…\dots

Proof

Proof by induction

This proof is adapted from Concrete Mathematics, 6.5, BERNOULLI NUMBER.

Use an identity for binomial coefficients together with induction:

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)\begin{aligned} S_{m+1}(n)+n^{m+1}&= \sum_{k=0}^{n-1}(k+1)^{m+1}\\ &=\sum_{k=0}^{n-1}\sum_{j=0}^{m+1}\binom{m+1}{j}k^j\\ &=\sum_{j=0}^{m+1}\binom{m+1}{j}S_j(n) \end{aligned}

Let S^m(n)=1m+1∑k=0m(m+1k)Bknm+1−k\hat{S}_{m}(n)=\frac{1}{m+1} \sum_{k=0}^{m} \binom{m+1}{k}B_kn^{m+1-k}. We want to prove Sm(n)=S^m(n)S_m(n)=\hat{S}_m(n), assuming that for j∈[0,m)j\in[0,m), Sj(n)=S^j(n)S_j(n)=\hat{S}_j(n).

Subtracting Sm+1(n)S_{m+1}(n) from both sides of the original equation gives:

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)\begin{aligned} S_{m+1}(n)+n^{m+1}&=\sum_{j=0}^{m+1}\binom{m+1}{j}S_j(n)\\ n^{m+1}&=\sum_{j=0}^{m}\binom{m+1}{j}S_j(n)\\ &=\sum_{j=0}^{m-1}\binom{m+1}{j}\hat{S}_j(n)+\binom{m+1}{m}S_m(n) \end{aligned}

Adding and subtracting (m+1m)S^m(n)−(m+1m)S^m(n)\binom{m+1}{m}\hat{S}_m(n)-\binom{m+1}{m}\hat{S}_m(n) on the right and simplifying gives:

nm+1=∑j=0m(m+1j)S^j(n)+(m+1)(Sm(n)−S^m(n))n^{m+1}=\sum_{j=0}^{m}\binom{m+1}{j}\hat{S}_j(n)+(m+1)(S_m(n)-\hat{S}_m(n))

Let Δ=Sm(n)−S^m(n)\Delta = S_m(n)-\hat{S}_m(n) and expand S^j(n)\hat{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)Δ\begin{aligned} n^{m+1}&=\sum_{j=0}^{m}\binom{m+1}{j}\hat{S}_j(n)+(m+1)\Delta\\ &=\sum_{j=0}^{m}\binom{m+1}{j}\frac{1}{j+1}\sum_{k=0}^{j}\binom{j+1}{k}B_kn^{j+1-k}+(m+1)\Delta\\ \end{aligned}

Reverse the order of summation in the second ∑\sum and apply identities for binomial coefficients:

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)Δ\begin{aligned} n^{m+1}&=\sum_{j=0}^{m}\binom{m+1}{j}\frac{1}{j+1}\sum_{k=0}^{j}\binom{j+1}{j-k}B_{j-k}n^{k+1}+(m+1)\Delta\\ &=\sum_{j=0}^{m}\binom{m+1}{j}\frac{1}{j+1}\sum_{k=0}^{j}\binom{j+1}{k+1}B_{j-k}n^{k+1}+(m+1)\Delta\\ &=\sum_{j=0}^{m}\binom{m+1}{j}\frac{1}{j+1}\sum_{k=0}^{j}\frac{j+1}{k+1}\binom{j}{k}B_{j-k}n^{k+1}+(m+1)\Delta\\ &=\sum_{j=0}^{m}\binom{m+1}{j}\sum_{k=0}^{j}\binom{j}{k}\frac{B_{j-k}}{k+1}n^{k+1}+(m+1)\Delta \end{aligned}

Interchanging the two summation signs gives:

nm+1=∑k=0mnk+1k+1∑j=km(m+1j)(jk)Bj−k+(m+1)Δn^{m+1}=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\sum_{j=k}^{m}\binom{m+1}{j}\binom{j}{k}B_{j-k}+(m+1)\Delta

Use the identity (m+1j)(jk)\binom{m+1}{j}\binom{j}{k}

(m+1j)(jk)＝(m+1k)(m−k+1j−k)\binom{m+1}{j}\binom{j}{k}＝\binom{m+1}{k}\binom{m-k+1}{j-k}

to obtain:

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)Δ\begin{aligned} n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\sum_{j=k}^{m}\binom{m+1}{k}\binom{m-k+1}{j-k}B_{j-k}+(m+1)\Delta\\ &=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}\sum_{j=k}^{m}\binom{m-k+1}{j-k}B_{j-k}+(m+1)\Delta\\ \end{aligned}

Replace every j−kj-k by jj:

nm+1=∑k=0mnk+1k+1(m+1k)∑j=0m−k(m−k+1j)Bj+(m+1)Δn^{m+1}=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}\sum_{j=0}^{m-k}\binom{m-k+1}{j}B_{j}+(m+1)\Delta

Recall the recurrence mentioned above:

∑j=0m(m+1j)Bj=0,(m>0)B0=1∑j=0m(m+1j)Bj=[m=0]\begin{aligned} \sum_{j=0}^{m}\binom{m+1}{j}B_j&=0,(m>0)\\ B_0&=1\\ \sum_{j=0}^{m}\binom{m+1}{j}B_j&=[m = 0] \end{aligned}

Substitution yields:

nm+1=∑k=0mnk+1k+1(m+1k)[m−k=0]+(m+1)Δ=nm+1m+1(m+1m)+(m+1)Δ=nm+1+(m+1)Δ\begin{aligned} n^{m+1}&=\sum_{k=0}^{m}\frac{n^{k+1}}{k+1}\binom{m+1}{k}[m - k = 0]+(m+1)\Delta\\ &=\frac{n^{m+1}}{m+1}\binom{m+1}{m}+(m+1)\Delta\\ &=n^{m+1}+(m+1)\Delta \end{aligned}

Thus Δ=0\Delta=0, and Sm(n)=S^m(n)S_m(n)=\hat{S}_m(n).

Proof using exponential generating functions

For the recurrence ∑j=0m(m+1j)Bj=[m=0]\sum_{j=0}^{m}\binom{m+1}{j}B_j=[m=0].

Add Bm+1B_{m + 1} to both sides:

∑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!\begin{aligned} \sum_{j=0}^{m+1}\binom{m+1}{j}B_j&=[m=0]+B_{m+1}\\ \sum_{j=0}^{m}\binom{m}{j}B_j&=[m=1]+B_{m}\\ \sum_{j=0}^{m}\dfrac{B_j}{j!}\cdot\dfrac{1}{(m-j)!}&=[m=1]+\dfrac{B_{m}}{m!} \end{aligned}

Let B(z)=∑i≥0Bii!ziB(z) = \sum\limits_{i\ge 0}\dfrac{B_i}{i!}z^i. The left side is a convolution, so:

B(z)ez=z+B(z)B(z)=zez−1\begin{aligned} B(z)\mathrm{e}^z &= z+B(z)\\ B(z)&=\dfrac{z}{\mathrm{e}^z - 1} \end{aligned}

Let Fn(z)=∑m≥0Sm(n)m!zmF_n(z) = \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m:

Fn(z)=∑m≥0Sm(n)m!zm=∑m≥0∑i=0n−1imzmm!\begin{aligned} F_n(z) &= \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m\\ &= \sum_{m\ge 0}\sum_{i=0}^{n-1}\dfrac{i^mz^m}{m!}\\ \end{aligned}

Interchange the order of summation:

Fn(z)=∑i=0n−1∑m≥0imzmm!=∑i=0n−1eiz=enz−1ez−1=zez−1⋅enz−1z\begin{aligned} F_n(z) &=\sum_{i=0}^{n-1}\sum_{m\ge 0}\dfrac{i^mz^m}{m!}\\ &=\sum_{i=0}^{n-1}\mathrm{e}^{iz}\\ &=\dfrac{\mathrm{e}^{nz} - 1}{\mathrm{e}^z - 1}\\ &=\dfrac{z}{\mathrm{e}^z - 1}\cdot\dfrac{\mathrm{e}^{nz} - 1}{z} \end{aligned}

Substitute B(z)=zez−1B(z)=\dfrac{z}{\mathrm{e}^z - 1}:

Fn(z)=B(z)⋅enz−1z=(∑i≥0Bii!)(∑i≥1nizi−1i!)=(∑i≥0Bii!)(∑i≥0ni+1zi(i+1)!)\begin{aligned} F_n(z) &= B(z)\cdot\dfrac{\mathrm{e}^{nz} - 1}{z}\\ &= \left(\sum_{i\ge 0}\dfrac{B_i}{i!} \right)\left(\sum_{i\ge 1}\dfrac{n^i z^{i - 1}}{i!}\right)\\ &= \left(\sum_{i\ge 0}\dfrac{B_i}{i!} \right)\left(\sum_{i\ge 0}\dfrac{n^{i+1} z^{i}}{(i+1)!}\right) \end{aligned}

Since Fn(z)=∑m≥0Sm(n)m!zmF_n(z) = \sum_{m\ge 0}\dfrac{S_m(n)}{m!}z^m, we have Sm(n)=m![zm]Fn(z)S_m(n)=m![z^m]F_n(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\begin{aligned} S \times m(n)&=m![z^m]F_n(z)\\ &= m!\sum_{i=0}^{m}\dfrac{B \times i}{i!}\cdot\dfrac{n^{m-i+1}}{(m-i+1)!}\\ &=\dfrac{1}{m+1}\sum_{i=0}^{m}\binom{m+1}{i}B_in^{m-i+1} \end{aligned}

This completes the proof.

Reference implementation 
using ll = long long;
constexpr int MAXN = 10000;
constexpr int mod = 1e9 + 7;
ll B[MAXN];        // Bernoulli numbers
ll C[MAXN][MAXN];  // binomial coefficients
ll inv[MAXN];      // modular inverses (for computing Bernoulli numbers)

void init() {
  // Precompute binomial coefficients
  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;
    }
  }
  // Precompute modular inverses
  inv[1] = 1;
  for (int i = 2; i < MAXN; i++) {
    inv[i] = (mod - mod / i) * inv[mod % i] % mod;
  }
  // Precompute Bernoulli numbers
  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;
  }
}
Previous section
230 Bell Numbers
Next section
232 Entringer Number

This content comes from OI Wiki: Bernoulli Numbers. See the original for authors and revision history. AIOJ has adapted the formatting, in-site links, and resource paths.

Content is shared under CC BY-SA 4.0 and SATA license; additional terms are governed by the original notice. Thanks to the OI Wiki authors and contributors; per the SATA terms, please visit the OI Wiki repository and star it to show your support.