A MODIFIED BERNOULLI NUMBER D. Zagier

The classical Bernoulli numbers B_{n}, defined by the generating function

\frac{x}{e^{x} − 1} = ∑_{n=0}^{∞} B_{n} \frac{x^{n}}{n!} ,  (1)

have many famous and beautiful properties, including the following three:

(i) B_{n} = 0 for odd n > 1. (ii) The fractional part of B_{n} is given by

B_{n} ≡ − ∑_{\begin{matrix}(p−1)|n \\ p prime \end{matrix}} \frac{1}{p} (mod 1) (n > 0, n even) .

(iii) B_{n} is given asymptotically for large even n by the formula

B_{n} ∼ (−1)^{n/2−1} \frac{n!}{2^{n−1}π^{n}} (n → ∞, n even).

In this note we show that the rational numbers defined by

B_{n}^{∗} = ∑_{r=0}^{n} \binom{n + r}{2r} \frac{B_{r}}{n + r} (n > 0) (2)

satisfy the following amusing variants of the above three properties: (I) The value of B_{n}^{∗} for n odd is periodic; more precisely, it is given by

\begin{matrix}n (mod 12) & 1 & 3 & 5 & 7 & 9 & 11 \\ B_{n}^{∗} & 3/4 & −1/4 & −1/4 & 1/4 & 1/4 & −3/4 \end{matrix}

(II) The fractional part of the number B̃_{n} := 2nB_{n}^{∗}−B_{n} is given by

B̃_{n} ≡ ∑_{\begin{matrix}(p+1)|n \\ p prime \end{matrix}} \frac{1}{p} (mod 1) (n > 0, n even) .

(III) B_{n}^{∗} is asymptotically equal to (−1)^{n/2−1} \frac{(n − 1)!}{(2π)^{n}} for n large and even, and is given much more precisely by the approximation

B_{n}^{∗} ≈ (−1)^{n/2} π Y_{n}(4π) (n → ∞, n even),

where Y_{n}(x) denotes the nth Bessel function of the second kind.

The proofs of (I)–(III) will be given in the next three sections, after which we will give the statement and proof of a fourth property, an exact formula for B_{n}^{∗} refining the asymptotic formula (III). The proofs, especially those of (I) and (II), are quite fun and the reader is invited to try to find them him/herself before proceeding.

We end the introduction with a small table of the numbers B_{n}^{∗} and B̃_{n}.

\begin{matrix}n & 2 & 4 & 6 & 8 & 10 & 12 & 14 & 16 & 18 & 20 \\ B_{n}^{∗} & \frac{1}{24} & −\frac{27}{80} & −\frac{29}{1260} & \frac{451}{1120} & −\frac{65}{264} & −\frac{6571}{12012} & \frac{571}{312} & −\frac{181613}{38080} & \frac{23663513}{1220940} & −\frac{10188203}{83600} \\ B̃_{n} & 0 & −\frac{8}{3} & −\frac{3}{10} & \frac{136}{21} & −5 & −\frac{4249}{330} & \frac{651}{13} & −\frac{3056}{21} & \frac{109269}{170} & −\frac{247700}{57} \end{matrix}

Proof of (I). Instead of using the familiar generating function (1), we represent the Bernoulli numbers by the generating function

F (x) = ∑_{r=1}^{∞} \frac{B_{r}}{r} x^{r} ∈ ℚ[[x]] .

This formal power series does not converge anywhere, but occurs in the asymptotic formula

\frac{Γ′(X)}{Γ(X)} ∼ log X − \frac{1}{X} − F (\frac{1}{X} ) (X → ∞)

for the logarithmic derivative of the gamma function, and the functional equation Γ(X + 1) = XΓ(X) of the gamma function implies the functional equation

F (\frac{x}{1 − x}) = F (x) + x + log(1 − x) ∈ ℚ[[x]]  (3)

of the power series F. An elementary proof of (3), or of the equivalent but simpler functional equation G(\frac{x}{1 − x}) = G(x) − x^{2} for the simpler power series G(x) =

∑_{r=0}^{∞} B_{r}x^{r+1} = x + x^{2}F′(x), can be obtained by noting that either one of these functional equations is equivalent to the standard recursion formula

∑_{k=0}^{n−1} \binom{n}{k}B_{k} = 0  (n > 1)  (4)

for the Bernoulli numbers, which is in turn an easy restatement of the definition (1). Now introduce a new power series F_{λ}(x), depending on a parameter λ, by

F_{λ}(x) = F (\frac{x}{1 − λx + x^{2}} ) − log(1 − λx + x^{2}) ∈ ℚ[[x]] .

For λ = 2 this specializes to

F_{2}(x) = ∑_{r=1}^{∞} \frac{B_{r}}{r} \frac{x^{r}}{(1 − x)^{2r}} − 2 log(1 − x) = 2 ∑_{n=1}^{∞} B_{n}^{∗}x^{n} .  (5)

On the other hand, the functional equation (3) together with the symmetry property F (−x) = F (x) + x, which is a restatement of (i), give the functional equations

F_{λ+1}(x) = F_{λ}(x) + \frac{x}{1 − λx + x^{2}} = F_{−λ}(−x)

for the power series F_{λ}. We deduce

F_{2}(x) − F_{2}(−x) = (F_{2}(x) − F_{1}(x)) + (F_{1}(x) − F_{0}(x)) + (F_{0}(x) − F_{−1}(x))

= \frac{x}{1 − x + x^{2}} + \frac{x}{1 + x^{2}} + \frac{x}{1 + x + x^{2}} = \frac{3x − x^{3} − x^{5} + x^{7} + x^{9} − 3x^{11}}{1 − x^{12}} .

Statement (I) follows.

Proof of (II). Rather surprisingly, this property is a consequence of the analogous property (ii) for the usual Bernoulli numbers, the divisibility by p − 1 being metamorphosed into the divisibility by p + 1 by the magic of generating functions. We begin by rewriting the definition of B̃_{n} as

B̃_{n} = ∑_{r=0}^{n−1} \frac{(n + r) + (n − r)}{n + r} \binom{n + r}{2r} B_{r} = β_{n} + β_{n−1} − B_{n} ,  (6)

where

β_{n} = ∑_{r=0}^{n} \binom{n + r}{2r} B_{r} (n ≥ 0) .

Fix a prime p. We want to show that pB̃_{n} is p-integral for all n and is congruent to 1 mod p if p + 1 divides n and to 0 otherwise. We suppose that p > 2. (The case p = 2 is similar but easier.) By (ii) we know that pB_{r} is p-integral and is congruent to −1 mod p if p − 1 divides r > 0 and to 0 otherwise. Equation (6) then immediately gives the p-integrality of pB̃_{n} and the congruence

pB̃_{n} ≡ −γ_{n} − γ_{n−1} + \{ \begin{matrix}1 & if (p − 1)|n \\ 0 & otherwise \end{matrix}  (7)

(here and from now on all congruences are modulo p), where

γ_{n} = ∑_{r>0, (p−1)|r} \binom{n + r}{2r} .

As usual, we use generating functions. From the definition of γ_{n} we have

∑_{n=0}^{∞} γ_{n} x^{n} = ∑_{r>0, (p−1)|r} \frac{x^{r}}{(1 − x)^{2r+1}} = \frac{x^{p−1}}{(1 − x)^{2p−1} − x^{p−1}(1 − x)} .

Hence (7) gives

∑_{n>0} pB̃_{n} x^{n} ≡ −(1 + x) ∑_{n=0}^{∞} γ_{n} x^{n} + \frac{x^{p−1}}{1 − x^{p−1}}

≡ −(1 + x) \frac{x^{p−1}(1 − x)}{(1 − x^{p})^{2} − x^{p−1}(1 − x)^{2}} + \frac{x^{p−1}}{1 − x^{p−1}} .

The expression on the right simplifies to \frac{x^{p+1}}{1 − x^{p+1}} , completing the proof.

Proof of (III). As in the proofs of (I) and (II), we deduce property (III) from its classical Bernoulli analogue. Suppose n is even and large. Then for r = n − k with k fixed (and even) we have

\binom{2n + r}{2r}\frac{B_{r}}{2n + r} = \frac{(−1)^{n/2−1}(n − 1)!}{(2π)^{n}} (\frac{(−1)^{k/2}(4π)^{k}}{k!} + O(\frac{1}{n})),

and since ∑ \frac{(−1)^{k/2}(4π)^{k}}{k!} = cos(4π) = 1 this gives the asymptotic formula in (III). The same argument in conjunction with the binomial coefficient identity

\binom{2n + r}{2r}\frac{1}{2n + r} = \frac{1}{2r!} ∑_{0≤h≤k/2} \frac{(−1)^{h}2^{k−2h}(n − h − 1)!}{h!(k − 2h)!} (r = n − k)

(which holds because both sides express the coefficient of x^{k} in (1 − 2x + x^{2})^{−r}/2r) lets us replace the asymptotic formula for B_{n}^{∗} by the full asymptotic expansion

(−1)^{n/2−1}B_{n}^{∗} ∼ ∑_{h≥0} \frac{(n − h − 1)!}{h!} (2π)^{−n+2h} .

Consulting standard reference works, we discover that the expression on the right is also the asymptotic development of −πY_{n}(4π), since Y_{n} is defined by

−π Y_{n}(2x) = ∑_{h=0}^{n−1} \frac{(n − h − 1)!}{h!}x^{−n+2h} + ∑_{h=0}^{∞}(c_{h} + c_{n+h} − 2 log x) \frac{(−1)^{h}x^{2h+n}}{h!(n + h)!}

with c_{h} = 1 + \frac{1}{2} + ⋯ + \frac{1}{h} − γ = O(log h). (Cf. [1], 7.2(32).) This proves the second assertion in (III), though without any estimate of the error.

Here are a few numerical values to illustrate the accuracy of the two approximations for B_{n}^{∗} in (III) :

\begin{matrix}n & 20 & 30 & 50 \\ (n − 1)! (2π)^{−n} & 13.228 & 10026347.89 & 0.75008667460769577 × 10^{23} \\ −π Y_{n}(4π) & 122.319 & 40532569.11 & 1.69052138468088825 × 10^{23} \\ (−1)^{n/2−1}B_{n}^{∗} & 121.868 & 40532573.81 & 1.69052138468090709 × 10^{23} \end{matrix}

The poorness of the approximations in the first row is explained by the above asymptotic expansion, which shows that the ratio of |B_{n}^{∗}| to (n − 1)!(2π)^{−n} is 1 + Cn^{−1} + O(n^{−2}) with C = 4π^{2} ≈ 40.

An exact formula. In the case of the usual Bernoulli numbers, the rough asymptotic formula (iii) can be replaced by the exact formula

B_{n} = \frac{(−1)^{n/2−1} n!}{2^{n−1} π^{n}} (1 + \frac{1}{2^{n}} + \frac{1}{3^{n}} + ⋯)  (n > 0, n even)  (8)

due to Euler. It is reasonable to look for a corresponding exact formula for B_{n}^{∗}.

We start with numerical data. From the table above we have B_{50}^{∗} + πY_{50}(4π) ≈ 1, 884 × 10^{9}, the exact value being 1884415006.56 . . .. Guessing that this difference might be related to the value −πY_{50}(8π), we compute this latter number, which turns out to be 1884414704.76 …, suggesting that we are on the right track. Going one step further, we find that the difference B_{50}^{∗} + πY_{50}(4π) + πY_{50}(8π) equals 301.79 …, rather close to the value −πY_{50}(12π) = 300.89 . . .. This suggests that B_{n}^{∗} might be very well approximated by, or even equal to, the number U_{n} = ∑_{ℓ=1}^{∞}(−1)^{n/2}πY_{n}(4πℓ). However, this sum diverges, since (−1)^{n/2}πY_{n}(4πℓ) behaves like −\frac{1}{2}ℓ^{−1/2} for ℓ large, so we must renormalize it, setting

U_{n} := ∑_{ℓ=1}^{∞}((−1)^{n/2}πY_{n}(4πℓ) + \frac{1}{2\sqrt{ℓ}}) − \frac{1}{2} ζ(\frac{1}{2}).  (9)

The series converges only like ∑ ℓ^{−3/2}, but can be replaced by the expression

U_{n} = ∑_{ℓ=1}^{∞}((−1)^{n/2}πY_{n}(4πℓ) + ∑_{k=0}^{K} \frac{c_{n,k}}{ℓ^{k+1/2}} ) − ∑_{k=0}^{K} c_{n,k} ζ(k + \frac{1}{2}) ,

c_{n,k} = \frac{(−1)^{[ \frac{k+1}{2} ]}}{2(8π)^{k} k!} (n + k − \frac{1}{2})(n + k − \frac{3}{2}) ⋯ (n − k + \frac{1}{2})

for any K ≥ 0, where the ℓth term is O(ℓ^{−K−3/2}), so the numerical value of U_{n} can be computed easily. Comparing it with B_{n}^{∗}, we find the following table:

\begin{matrix}n & B_{n}^{∗} − U_{n} \\ 2 & −0.886968175 … \\ 4 & −1.988273972 … \\ 10 & −4.999969846 … \\ 30 & −15 + 1.292 … × 10^{−13} \\ 50 & −25 + 5.646 … × 10^{−22} \end{matrix}

suggesting a formula of the form

B_{n}^{∗} + \frac{n}{2} = U_{n} + ε_{n}  (n > 0, n even),  (10)

where ε_{n} is positive and goes to 0 rapidly as n tends to infinity. (Notice that the quantity B_{n}^{∗} + n/2 occurring on the left is given by the same expression as in (2), but with B_{r} replaced by (−1)^{r}B_{r}.) We now prove this and give an exact formula for the error term ε_{n}.

(IV) Define B_{n}^{∗} by (2) and U_{n} by (9). Then equation (10) holds with

ε_{n} = ∑_{k=1}^{∞} \frac{1}{\sqrt{k(k + 4)}} (\frac{\sqrt{k} + \sqrt{k + 4}}{2})^{−2n} ∼ \frac{1}{\sqrt{5}} (\frac{3 + \sqrt{5}}{2})^{−n} . (11)

Proof of (IV). Using the relation (8) between Bernoulli numbers and zeta values, we find

B_{n}^{∗} = \frac{1}{n} − \frac{n}{4} − ∑_{\begin{matrix}2≤r≤n \\ r even \end{matrix}} \binom{n + r − 1}{2r − 1} \frac{(r − 1)!}{(2πi)^{r}} ζ(r)

and hence

B_{n}^{∗} − U_{n} = \frac{1}{n} − \frac{n}{4} + \frac{1}{2}ζ(\frac{1}{2}) + ∑_{ℓ=1}^{∞} b_{n,ℓ} ,  (12)

where

b_{n,ℓ} := (−1)^{n/2−1} π Y_{n}(4πℓ) − ∑_{\begin{matrix}2≤r≤n \\ r even \end{matrix}} \binom{n + r − 1}{2r − 1} \frac{(r − 1)!}{(2πiℓ)^{r}} − \frac{1}{2\sqrt{ℓ}} .

On the other hand, standard formulas for Bessel functions (cf. [1], 7.2(15) and 7.3(16)) give

(−1)^{n/2−1}π Y_{n}(x) = 2 ℜ(K_{n}(ix))

= \frac{\sqrt{π/2}}{Γ(n + \frac{1}{2})} ∫_{0}^{∞} \frac{e^{−ix}(t/2 + ix)^{n− \frac{1}{2}} + e^{ix}(t/2 − ix)^{n− \frac{1}{2}}}{(ix)^{n}} e^{−t} t^{n− \frac{1}{2}} dt

for n even and positive, and hence, after some simple manipulations, the formula

b_{n,ℓ} = \frac{2 Γ(\frac{1}{2})}{Γ(n + \frac{1}{2})} ∫_{0}^{∞} t^{n−1} e^{−t} (f_{n}(\frac{8πℓ}{t}) − \frac{1}{4} \sqrt{\frac{t}{πℓ}} ) dt  (13)

where f_{n}(x) is defined for x > 0 by

f_{n}(x) = \frac{(1 + ix)^{n− \frac{1}{2}} + (1 − ix)^{n− \frac{1}{2}}}{2 (ix)^{n}} − ∑_{\begin{matrix}2≤r≤n \\ r even \end{matrix}} \binom{n − \frac{1}{2}}{n − r}(ix)^{−r} .

Note that f(x) = (2x)^{−1/2} + O(x^{−3/2}) as x → ∞, so that b_{n,ℓ} = O(ℓ^{−3/2}), as we already know. Clearly f_{n}(x) extends to ℂ^{0} := ℂ ∖ (−i∞, −i] ∖ [i, i∞) as an even holomorphic function, and the binomial theorem gives the Taylor expansion

f_{n}(x) = ∑_{r≥0}(−1)^{r} \binom{n − \frac{1}{2}}{n + 2r} x^{2r} (|x| < 1) .

The beta integral identity

\binom{n − \frac{1}{2}}{n + 2r} = \frac{1}{π} \frac{Γ(n + \frac{1}{2})Γ(2r + \frac{1}{2})}{Γ(n + 2r + 1)} = \frac{1}{π} ∫_{0}^{1} u^{2r− \frac{1}{2}} (1 − u)^{n− \frac{1}{2}} du

now gives the integral representation

f_{n}(x) = \frac{1}{π} ∫_{0}^{1} \frac{u^{− \frac{1}{2}} (1 − u)^{n− \frac{1}{2}}}{1 + x^{2}u^{2}} du

for |x| < 1, and by analytic continuation this holds for all x ∈ ℂ^{0}. Consider the more general integral

f_{n}(x, s) = \frac{1}{π} ∫_{0}^{1} \frac{u^{s−1} (1 − u)^{n− \frac{1}{2}}}{1 + x^{2}u^{2}} du (x ∈ ℂ^{0}, s ∈ ℂ, 0 < ℜ(s) < 2) ,

so that f_{n}(x) = f_{n}(x, \frac{1}{2}). It can be estimated for large x by

f_{n}(x, s) = \frac{1}{π} ∫_{0}^{∞} \frac{u^{s−1} du}{1 + x^{2}u^{2}} + O(x^{−s−1}) = \frac{x^{−s}}{2 sin(πs/2)} + O(x^{−s−1}) .

In particular, ∑_{ℓ} f_{n}(ℓx, s) converges for ℜ(s) > 1 and all x > 0, and using the Poisson summation identity

∑_{ℓ=−∞}^{∞} \frac{1}{1 + ℓ^{2} t^{2}} = \frac{π/t}{tanh π/t} = \frac{π}{t} (1 + 2 ∑_{k=1}^{∞} e^{−2πk/t}) (t > 0)

we find

∑_{ℓ=−∞}^{∞} f_{n}(ℓx, s) = \frac{1}{x} (\frac{Γ(s − 1)Γ(n + \frac{1}{2})}{Γ(s + n − \frac{1}{2})} + 2 ∑_{k=1}^{∞} ∫_{0}^{1} u^{s−2} (1 − u)^{n− \frac{1}{2}} e^{−2πk/xu} du)

for ℜ(s) > 1. Writing the left-hand side of this identity as

f_{n}(0, s) + ∑_{ℓ=1}^{∞}(2f_{n}(ℓx, s) − \frac{(ℓx)^{−s}}{sin(πs/2)}) + \frac{ζ(s) x^{−s}}{sin(πs/2)}

gives its analytic continuation to ℜ(s) > 0, and setting s = \frac{1}{2}, x = \frac{8π}{t} we obtain

∑_{ℓ=1}^{∞}(f_{n}(\frac{8πℓ}{t}) − \frac{1}{4}\sqrt{\frac{t}{πℓ}}) = − \frac{Γ(n + \frac{1}{2})}{8Γ(\frac{1}{2})Γ(n)} t − \frac{ζ(\frac{1}{2})}{4\sqrt{π}} t^{\frac{1}{2}} − \frac{Γ(n + \frac{1}{2})}{2Γ(\frac{1}{2})Γ(n + 1)}

+ \frac{t}{8π} ∑_{k=1}^{∞} ∫_{0}^{1} u^{− \frac{3}{2}} (1 − u)^{n− \frac{1}{2}} e^{−kt/4u} du .

Combining this with (13) and performing the integrations over t we get

∑_{ℓ=1}^{∞} b_{n,ℓ} = −\frac{n}{4} − \frac{ζ(\frac{1}{2})}{2} − \frac{1}{n} + \frac{Γ(n + 1)}{4Γ(\frac{1}{2})Γ(n + \frac{1}{2})} ∑_{k=1}^{∞} ∫_{0}^{1} \frac{u^{− \frac{3}{2}} (1 − u)^{n− \frac{1}{2}}}{(1 + k/4u)^{n+1}} du .

The desired result now follows by inserting this into (12) and using the identity

\frac{Γ(n + 1)}{4Γ(\frac{1}{2})Γ(n + \frac{1}{2})} ∫_{0}^{1} \frac{u^{− \frac{3}{2}} (1 − u)^{n− \frac{1}{2}}}{(1 + k/4u)^{n+1}} du = \frac{1}{\sqrt{k(k + 4)}} (\frac{\sqrt{k} + \sqrt{k + 4}}{2})^{−2n} ,

which can be proved either from standard hypergeometric formulas or by expanding both sides as power series in 1/k for k > 4 and then using analytic continuation.

Remarks. We end this paper with a number of remarks.

1. The formula given in (i), which can be rewritten in the form

B_{n}^{∗} = \frac{1}{2} (\frac{−3}{n}) + \frac{1}{4} (\frac{−4}{n}) (n odd) ,

occurred originally in [4] in the context of the Eichler-Selberg trace formula for the traces of the Hecke operator T_{ℓ} acting on modular forms on SL_{2}(ℤ). The method of proof there gave a formula for these traces which had a somewhat different form from the classical one, and in particular involved Bernoulli numbers. The specialization of the formula to the case ℓ = 1 gave the dimension of M_{k}(SL_{2}(ℤ)) in terms of B_{k−1}^{∗}, and the equality of this expression with the standard dimension formula required the periodicity property (I).

2. The same idea as was used to prove (I) leads to another simple expression for the modified Bernoulli numbers B_{n}^{∗}. Specifically, from

F_{2}(x) = F_{0}(x) + \frac{x}{1 − x + x^{2}} + \frac{x}{1 + x^{2}}

and the binomial formula we get the identity

2B_{2n}^{∗} = (\frac{−3}{n}) + ∑_{r=0}^{n}(−1)^{n+r}\binom{n + r}{2r}\frac{B_{2r}}{n + r} ,

whose second term has a pleasing similarity to the original formula (2) defining B_{n}^{∗}.

3. Next, we mention that the definition (2) can also be inverted to express the ordinary Bernoulli numbers in terms of the modified ones, should we for any reason wish to do so. Indeed, from (5) we have

∑_{n=1}^{∞} \frac{B_{n}}{n} x^{n} = 2 log(\frac{\sqrt{1 + 4x} − 1}{2x}) + 2 ∑_{r=1}^{∞} B_{r}^{∗} (\frac{1 + 2x − \sqrt{1 + 4x}}{2x})^{r} ,

and comparing coefficients of B_{n} we get

B_{n} = (−1)^{n}\binom{2n}{n} + 2 ∑_{r=1}^{n}(−1)^{n−r} r\binom{2n}{n − r} B_{r}^{∗} .

4. Finally, we mention that there are (at least) two ways of calculating Bernoulli numbers which are faster than the standard recursion (4). The first is due to M. Kaneko [3]. The second I noticed myself, but Kaneko has informed me that it is in fact a classical identity going back to Kronecker. (See [2] for a historical survey.) Nevertheless, these formulas are both pretty and useful, so for the sake of popularization we reproduce them here. a. One can replace (4) by a recursion of the same type, but with only half as many terms, namely, setting b_{n} = (n + 1)B_{n},

b_{2n} = −\frac{1}{n + 1} ∑_{i=0}^{n−1} \binom{n + 1}{i} b_{n+i}  (14)

together with the conditions b_{1} = −1 and b_{2n+1} = 0 for n > 0. Equation (14) can be seen as a special case of the following fact: Define an involution ^{∗} on the set of sequences {a_{0}, a_{1}, a_{2}, … } by A^{∗}(x) = e^{−x} A(−x), where A(x) := ∑_{n=0}^{∞} a_{n}x^{n}/(n + 1)!, or more explicitly by a_{n}^{∗} = (−1)^{n} ∑_{i=0}^{n} \binom{n+1}{i+1}a_{i}. Then the expression ∑_{i=0}^{n} \binom{n}{i}a_{n+i−1} (n ≥ 1) is anti-invariant under ^{∗} and hence vanishes if A^{∗} = A. (Note that A(x) = x/(e^{x} − 1) = A^{∗}(x) for a_{n} =b_{n}.)

b. The Bernoulli numbers can be calculated directly, rather than recursively, by the closed formula

B_{k} = ∑_{n=0}^{k}(−1)^{n}\binom{k + 1}{n + 1}\frac{0^{k} + ⋯ + n^{k}}{n + 1} .  (15)

To prove this, we apply Bernoulli’s famous formula for 1^{k} + ⋯ + n^{k} to get

\frac{0^{k} + ⋯ + n^{k}}{n + 1} = \frac{B_{k+1}(n + 1) − B_{k+1}}{(k + 1)(n + 1)} = (polynomial of degree k in n) + \frac{B_{k}}{n + 1} ,

where B_{r}(x) denotes the rth Bernoulli polynomial; taking the (k + 1)st difference of both sides kills the polynomial on the right, leaving only an easily computed multiple of B_{k}. Formula (15) is much more convenient for numerical computations than the recursion formula, at least if one wants to compute individual Bernoulli numbers rather than a table up to some limit, since the number of steps needed to compute B_{k} is O(k) rather than O(k^{2}) (each term in the sum can be computed from its predecessor in O(1) steps). Indeed, even for computing a table, (15) is sometimes more useful than the recursion (4), since the time required is about the same but the storage requirements are reduced from O(k) to O(1). Here is a one-line PARI program implementing the formula (15) (for k > 0):

 B(k)=h=0;s=1;c=k+1;for(n=2,k+1,c=c*(n-k-2)/n;h=h+c*s/n;s=s+n^k);h

References.

[1] A. Erdélyi, W. Magnus, F. Oberhettinger and F. Tricomi, Higher Transcendental Functions (Bateman Manuscript Project), Vol. II, McGraw-Hill, New York (1953)

[2] H.W. Gould, Explicit formulas for Bernoulli numbers, Amer. Math. Monthly 79 (1972) 44–51

[3] M. Kaneko, A recurrence formula for the Bernoulli numbers, Proc. J. Acad. 71 (1995) 192–193

[4] D. Zagier, Hecke operators and periods of modular forms, Israel Math. Conf. Proc. 3 (1990) 321–336