This site is supported by donations to The OEIS Foundation.

The computation and asymptotics of the Bernoulli numbers

This is a sequel of this blog which discusses the Seidel transform.

The difference table of the Bernoulli numbers

The difference table of the Bernoulli numbers is listed on OEIS as A085737/A085738 (apart from signs).

The Bernoulli median numbers

Paul Curtz commented in A181130 : "(-1)^n*a(n) is the numerator on the main diagonal of ... A168516." But much more is true: Vladimir Reshetnikov's A181130/A181131 is the main diagonal of A085737/A085738, in other words, of the difference table of the Bernoulli numbers. And because of this A181130/A181131 are the Bernoulli median numbers. Almost. The signs are missing and the offset should be 0 with value 1. Reshetnikov's definition cannot be extended to n=0 because the defining integral is divergent for n=0.

So let us work out our own variant.

• Name: Bernoulli median numbers.
• Offset: 0
• Numerators:     A212196 = 1, -1, 2, -8, 8, -32, 6112, -3712, ...
• Denominators: A181131 = 1, 3, 15, 105, 105, 231, 15015, 2145, ...
• Definition: The median (central) column in the difference table of the Bernoulli numbers.
• Computation: According to L. Seidel.
def A212196/A181131_list(n) :
def T(S, a) :
R = [a]
for s in S :
a -= s
R.append(a)
return R
def M(A, p) :
R = T(A,0)
S = add(r for r in R)
return -S / (2*p+3)
R = ; A = [1/2, -1/2]
for k in (0..n-2) :
A = T(A, M(A,k))
R.append(A[k+1])
A = T(A,0)
return R

A212196/A181131_list(10)
[1, -1/3, 2/15, -8/105, 8/105, -32/231, 6112/15015, -3712/2145,
362624/36465, -71706112/969969]

Of course all the fun and profit of this algorithm is to build the difference table of the Bernoulli numbers without using the Bernoulli numbers; everything else would be pretty lame. In fact Seidel used this algorithm to compute the Bernoulli numbers themselves (the present author only changed the output of the algorithm). More on this below.

The difference table of the Bernoulli polynomials

Seidel's computing scheme can be extended from the case of integers to the case of polynomials. Recall first the Bernoulli polynomials which we give here also in the more convenient form as polynomials with integer coefficients, A2113615, divided by the denominators of the Bernoulli polynomials, B(n, x) = A2113615(n, x) / A144845(n).

 n B(n, x) A2113615(n, x) = A144845(n)*B(n, x) 0 1 1 1 x-1/2 2*x-1 2 x^2-x+1/6 6*x^2-6*x+1 3 x^3-3/2*x^2+1/2*x 2*x^3-3*x^2+x 4 x^4-2*x^3+x^2-1/30 30*x^4-60*x^3+30*x^2-1

Next let us consider the difference table of the Bernoulli polynomials:

Evaluating the table at x = 1 gives precisely the tableau given by Seidel and displayed above.

(Note by the way that this would of course not work if we choose to evaluate the table at x = 0. This is just one of the many reasons why the definition Bn = Bn(0) does not fit into the landscape of things. For an account on this point see my Bernoulli Manifesto.

The difference table motivates the introduction of the Bernoulli nabla polynomials B(n, x) which are in the first row of the table. This sequence begins:

 n B∇(n, x) A2113616(n, x) = A144845(n)*B∇(n, x) 0 1 1 1 x-3/2 2*x-3 2 x^2-3*x+13/6 6*x^2-18*x+13 3 x^3-9/2*x^2+13/2*x-3 2*x^3-9*x^2+13*x-6 4 x^4-6*x^3+13*x^2-12*x+119/30 30*x^4-180*x^3+390*x^2-360*x+119

From the well known formula DLMF-24.2.5 we immediately see that the sequence of Bernoulli nabla polynomials in its integer form A2113616 can be simply calculated by (Maple):

seq(print(sort(expand(denom(bernoulli(i,t))*bernoulli(i,x-1)))),i=0..4);

The difference table can be computed with this Maple procedure:

DifferenceTableBernoulliPolynomials := proc(n) local A,m,k;

# Clear matrix
A := array(0..n,0..n);
for m from 0 to n do for k from 0 to n do A[m,k] := '~' od od;

# Compute elements
for m from 0 to n do A[m,0] := bernoulli(m,x);
for k from m-1 by -1 to 0 do
A[k,m-k] := A[k+1,m-k-1] - A[k,m-k-1] od od;

# Beautify elements
for m from 0 to n do for k from 0 to n do
A[m,k] := sort(expand(A[m,k])) od od;
convert(A, matrix) end:

DifferenceTableBernoulliPolynomials(5);

With Sage:

def DifferenceTableBernoulliPolynomials(n) :
R.<x> = PolynomialRing(QQ)
T = matrix(R, n)
for m in (0..n-1) :
T[m,0] = bernoulli_polynomial(x,m)
for k in range(m-1,-1,-1) :
T[k,m-k] = T[k+1,m-k-1] - T[k,m-k-1]
return T
DifferenceTableBernoulliPolynomials(5)

The Bernoulli median polynomials

Finally we indicate the Bernoulli median polynomials which are given by the main diagonal of the difference table.

 Bμ(0, x) 1 Bμ(1, x) x^2-2*x+2/3 Bμ(2, x) x^4-4*x^3+5*x^2-2*x+2/15 Bμ(3, x) x^6-6*x^5+13*x^4-12*x^3+4*x^2-8/105 Bμ(4, x) x^8-8*x^7+74/3*x^6-36*x^5+71/3*x^4-4*x^3-4/3*x^2+8/105

Note that evaluating the Bernoulli median polynomials at x = 1 gives the Bernoulli median numbers; evaluating these polynomials at x = 0 to compute the median numbers would be erroneous because Bμ(1, 0) = 2/3  <>  β1 = -1/3.

The connection with the Eulerian polynomials

The definition of the Eulerian polynomials A(n, x) which we will use here is given (via the Eulerian numbers) in A173018. A short introduction can be found in this wiki.

Now let us introduce the rational functions

 $P_{0}(x)={\frac {{x}^{2}}{\left(x+1\right)^{2}}}$ $P_{1}(x)={\frac {-\ {x}^{2}}{\left(x+1\right)^{4}}}$ $P_{2}(x)={\frac {{x}^{2}}{\left(x+1\right)^{6}}}\left(-x+1\right)^{2}$ $P_{3}(x)={\frac {-\ {x}^{2}}{\left(x+1\right)^{8}}}\left({x}^{2}-4\,x+1\right)^{2}$ $P_{4}(x)={\frac {{x}^{2}}{\left(x+1\right)^{10}}}\left(-{x}^{3}+11\,{x}^{2}-11\,x+1\right)^{2}$ More generally these functions can be defined with the Eulerian polynomials A(n, x) as Pn(x) = (-1)^n*x^2*A(n,-x)^2/(1+x)^(2*(n+1)). Then, using Reshetnikov's integral, for n > 0, the integral of Pn(x) over the positive reals is given by the Bernoulli median numbers βn.

$\int _{0}^{\infty }{\frac {\left(-1\right)^{n}{x}^{2}\left(A(n,-x)\right)^{2}}{\left(1+x\right)^{2\,n+2}}}=\beta _{n}$ With Maple:

a := proc(n,m) local k; add((-1)^k*binomial(n+1,k)*(m+1-k)^n,k=0..m) end:
A := proc(n,x) local k; if(n=0,1,add(a(n,k)*x^k,k=0..n-1)) end:
P := proc(n,x) (-1)^n*(x*A(n,-x)/(1+x)^(n+1))^2 end:

Finally we note that (1/4)*Pn(-2) is a signed version of Vladeta Jovovic's A122725.

Bernoulli numbers via Seidel's Treppen Schema

Facsimile from Ludwig Seidel, Ueber eine einfache Entstehungsweise der Bernoulli'schen Zahlen

Seidel's Genocchi triangle, A014781

Here is how Seidel computed the 'Treppen-Schema' in the facsimile. (Note that the parameter n is the number of requested rows in the triangles, which are '0'-based. Thus T(n) returns row(0),..,row(n-1).)

def SeidelGenocchiTriangle(n) :
D = *(n//2+3); D = 1
b = True; h = 1
for i in range(n) :
if b :
for k in range(h,0,-1) : D[k] += D[k+1]
a = 1; h +=1
else :
for k in range(1,h, 1) : D[k] += D[k-1]
a = 0
b = not b
print i, [D[z] for z in (a..h-1)]

SeidelGenocchiTriangle(14)
0 
1 [0, 1]
2 [1, 0]
3 [0, 1, 1]
4 [2, 1, 0]
5 [0, 2, 3, 3]
6 [8, 6, 3, 0]
7 [0, 8, 14, 17, 17]
8 [56, 48, 34, 17, 0]
9 [0, 56, 104, 138, 155, 155]
10 [608, 552, 448, 310, 155, 0]
11 [0, 608, 1160, 1608, 1918, 2073, 2073]
12 [9440, 8832, 7672, 6064, 4146, 2073, 0]
13 [0, 9440, 18272, 25944, 32008, 36154, 38227, 38227]

Unsigned Genocchi numbers (even index), A110501

Two special cases, the unsigned Genocchi numbers and the Genocchi median numbers can be easily extracted from this scheme.

def A110501_list(n) :
D = *(n+2); D = 1
R = ; b = True; h = 2
for i in range(2*n-2) :
if b :
for k in range(h,0,-1) : D[k] += D[k+1]
h += 1
else :
for k in range(1,h, 1) : D[k] += D[k-1]
b = not b
if b : R.append(D[h-2])
return R

A110501_list(13)
[1, 1, 3, 17, 155, 2073, 38227, 929569, 28820619, 1109652905,
51943281731, 2905151042481, 191329672483963]

Genocchi numbers of second kind (or Genocchi medians), A005439

def A005439_list(n) :
D = *(n+2); D = 1
R =  ; b = True; h = 2
for i in range(2*n-3) :
if b :
for k in range(h,0,-1) : D[k] += D[k+1]
h += 1
else :
for k in range(1,h, 1) : D[k] += D[k-1]
if b : R.append(D)
b = not b
return R

A005439_list(13)
[1, 1, 2, 8, 56, 608, 9440, 198272, 5410688, 186043904, 7867739648,
401293838336, 24290513745920]

The notion of the Genocchi median numbers immediately suggest a riddle for the curious reader: "What are the Bernoulli median numbers? Can you spot them in the database? " (Hint: They were only added in 2011 to OEIS! And of course they are also in Seidel's paper.)

The 'interleaved' Genocchi numbers, A099960

def A099960_list(n) :
D = *(n//2+3); D = 1
R = []; b = True; h = 1
for i in (1..n) :
if b :
for k in range(h,0,-1) : D[k] += D[k+1]
R.append(D); h += 1
else :
for k in range(1,h, 1) : D[k] += D[k-1]
R.append(D[h-1])
b = not b
return R

A099960_list(14)
[1, 1, 1, 1, 2, 3, 8, 17, 56, 155, 608, 2073, 9440, 38227]

Bernoulli numbers after Seidel

To make Seidel's basic algorithm more efficient one can bring two ideas into action: to avoid rational arithmetic in the computation and to exploit the symmetry of the difference tableau, i. e. to compute only one half of the antidiagonals. (This explains, by the way, why the 'median' (or 'central') numbers appear at the left hand side of Seidel's staircase scheme.) Together this leads to Seidel's 'Treppen-Schema' for computing B2n for n > 0.

def Bernoulli_list(n) :
G = A110501_list(n)
# first formula in the facsimile above
return [(-1)^i*G[i-1]/(2*(1-4^i)) for i in (1..n)]

Bernoulli_list(13)
[1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798,
-174611/330, 854513/138, -236364091/2730, 8553103/6]

All in one (Maple):

SeidelBernoulli := proc(n) local D, R, i, b, h, k, p, q, s;
D := array(0..n, [seq(0,i=0..n)]); D := 1;
R := 1; b := true; h := 0; p := 1; s := -2;
for i from 0 to 2*n-3 do
if b then h := h + 1; p := 4*p; s := -s; q := s*(p-1);
for k from h by -1 to 1 do D[k] := D[k] + D[k+1] od
else  for k from 1 by  1 to h do D[k] := D[k] + D[k-1] od;
R := R, D[h] / q;
fi;
b := not b;
od;
R end:

SeidelBernoulli(13)
[1, 1/6, -1/30, 1/42, -1/30, 5/66, -691/2730, 7/6, -3617/510, 43867/798,
-174611/330, 854513/138, -236364091/2730]

Zooming into the function we observe in its center the lines

for k from h by -1 to 1 do D[k] := D[k] + D[k+1] od
for k from 1 by  1 to h do D[k] := D[k] + D[k-1] od;

We see here the moment of birth of the boustrophedon transformation.

Indeed Seidel's approach is rather ingenious and Seidel was quite fond of his approach; he closed his paper with the words:

"[Die Bernoulli'schen Zahlen sind] nicht blos privilegiert durch ihre Rolle in der Analysis, sondern auch ausgezeichnet durch ihre arithmetische Natur selbst, vermöge deren sie sich in einfacher und doch charakteristischer Weise sozusagen von selbst aus den Grund-Elementen 1 und 0 aller Zahlenbetrachtung entfalten."

"[The Bernoulli numbers are] not only privileged by their role in the analysis, but also distinguished by their arithmetic nature itself, by virtue of which they unfold in a simple and yet characteristical manner from the basic elements 1 and 0 of all number contemplation."

Let us illustrate this meditation. At the beginning the stars denote the unknown entries .

1
*     0
0     *     *
*      *     *      0
0      *     *     *      *
*     *      *     *      *     0

What Seidel describes in the words "they unfold in a simple and yet characteristical manner" is nowadays known as Seidel's algorithm.

Algorithms for computing the Bernoulli numbers

There are four simple algorithms for computing the Bernoulli numbers which we will compare now.

• Seidel's algorithm (1877) 
• Akiyama & Tanigawa (1999) 
• Luschny (2010) 
• Brent & Harvey (2011) 

The method of Knuth and Buckholtz (1967)  is very similar to the Brent and Harvey method which is a little bit simpler and faster; thus their seems to be no point to consider the K&B method separately. The implementations below all compute the list of the first n even Bernoulli numbers.

Seidel

def SeidelBernoulli(n) :
D = *(n+2); D = 1; b = True;
R = ; h = 1; p = 1; s = -2
for i in range(2*n-2) :
if b :
h += 1; p = 4*p; s = -s; q = s*(p-1)
for k in range(h,0,-1) : D[k] += D[k+1]
else :
for k in range(1,h,1) : D[k] += D[k-1]
R.append(D[h-1]/q)
b = not b
return R

Akiyama-Tanigawa

def AkiyamaTanigawaBernoulli(n) :
A = *(2*n-1)
R = []; ev = True
for m in range(2*n-1) :
A[m] = 1/(m + 1)
for j in range(m, 0, -1) :
A[j - 1] = j * (A[j - 1] - A[j])
if ev : R.append(A)
ev = not ev
return R

Luschny

def LuschnyBernoulli(n) :
A = *(2*n); A = 1;
R = ; ev = True
for m in (2..2*n-1) :
A[m] = 0
for j in range(m, 0, -1) :
A[j] = j * (A[j - 1] + A[j])
if ev : R.append(add((-1)^j*A[j]/(j+1) for j in (1..m)))
ev = not ev
return R

Brent-Harvey

This routine is very similar to the Knuth-Buckholtz method (see formula 6.93 in GKP).

def BrentHarveyBernoulli(n) :
T = *(n+1); T = 1
for k in (2..n) :
T[k] = (k-1)*T[k-1]
for k in (2..n-1) :
for j in (k..n) :
T[j] = (j-k)*T[j-1]+(j-k+2)*T[j]

R = ; E = -2; U = 1
for k in (1..n-1) :
U = 4*U; A = U*(U-1); E = -E
R.append(E*T[k]*k/A)
return R

Benchmarks

 Algorithm n Maple (sec) Sage (cpu/wall sec) Seidel 1000 20 1.7 / 2.5 Brent-Harvey-Knuth-Buckholz 1000 44 1.9 / 2.7 Akiyama-Tanigawa 500 175 8.7 / 12.2 Luschny 500 33 2.6 / 3.6

To sum up: Seidel's computation of the Bernoulli numbers is based on the Genocchi numbers, the computations of Knuth-Buckholtz and Brent-Harvey are based on the Tangent numbers, the Luschny algorithm is based on the Worpitzky numbers. The most striking algorithm is without doubt the Akiyama-Tanigawa algorithm. The poor performance of this algorithm is due to the fact that all operations are rational operations; in contrast the other algorithms avoid rational operations with the exception of a final division.

In any case it is nice to see that the algorithm of some not so unknown 20th and 21th century experts cannot beat Seidel's 19th century staircase scheme. (To be fair: for even larger n  Harvey and Brent-Harvey did develop faster methods. But they are also much more involved.)

Asymptotic approximations of the Bernoulli numbers

Let us look at three asymptotic formulas for the logarithm of the Bernoulli numbers (see also ).

1. LogB1(n)=(1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n
2. LogB2(n)=(1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*
(1-ln((120*n^2+9)/(120*n^2-1)))
3. LogB3(n)=(1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*
(1-1/(12*n^2)*(1-1/(30*n^2)*(1-2/(7*n^2))))

${\text{logB3}}(n)=\left({\frac {1}{2}}+n\right)\ln(n)+\left({\frac {1}{2}}-n\right)\ln(\pi )+\left({\frac {3}{2}}-n\right)\ln(2)-R(n)$ where

$R(n)=n\left(1-{\frac {1}{12}}\,\left(1-{\frac {1}{30}}\,\left(1-{\frac {2}{7}}\,{n}^{-2}\right){n}^{-2}\right){n}^{-2}\right)$ We see that formula 2 as well as formula 3 are refinements of formula 1. We are now focusing on formula 3 in its exponential form exp(LogB3(n)).

What makes this asymptotic approximation especially useful — besides being a much better approximation than those given by the Digital Library of Mathematical Functions (see §24.11) — is that the error of the approximation can be easily estimated.

In fact a convenient way to reason about the validity of an approximation formula is to give a lower bound for the number of exact decimal digits, i. e. to indicate the number of decimal digits which are guaranteed by the formula at least. In the case of exp(LogB3(n)) we have a good and simple way to express this bound:

EddE(n) = floor(3*log(3*n))    (valid for n ≥ 50).

For example for Bernoulli(31622776) this bound says that 55 decimal digits of exp(logB3(n)) are guaranteed to be correct (the true value is 55.7).

To sum up: to compute an approximation to the Bernoulli numbers B(n) with n ≥ 50 and n even compute exp(logB3(n)) and retain floor(3*log(3*n)) decimal digits of the result.

def BernoulliAsympt(n) :
if n < 50 :
print "Value error, n has to be ≥ 50"
return
if is_odd(n): return 0
R = RealField(300)  # to be increased for large values n
nn = n*n
LogB = R((1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n*(1-1/(12*nn)*
(1-1/(30*nn)*(1-2/(7*nn))))
B = (-1)^(1+n//2)*exp(LogB)
Edd = floor(3*ln(3*n))
print SciFormat(B, Edd - 1)
return B

for n in (49..52): print n, BernoulliAsympt(n)
49 → Value error, n has to be ≥ 50
50 → 7.50086674607696e24
51 → 0
52 → -5.0387781014810e26

BernoulliAsympt(31622776)
-7.66922063003368519879820408820523505875084131626143035e198196563

Note that the function displays only valid digits as we used this formatting function: SciFormat. With Sage the asymptotic formulas for the Bernoulli numbers based on these formulas and the exact decimal digits of these approximations can be computed as follows.

Edd1(n) = -log10(abs(1-exp(LogB1(n))/abs(bernoulli(n)))
Edd2(n) = -log10(abs(1-exp(LogB2(n))/abs(bernoulli(n)))
Edd3(n) = -log10(abs(1-exp(LogB3(n))/abs(bernoulli(n)))
EddE(n) = 3*ln(3*n)

ANF = 50; END = 1000; STEP = 20
plot2 = list_plot([[n, Edd2(n)] for n in range(ANF,END,STEP)], color='red')
plot3 = list_plot([[n, Edd3(n)] for n in range(ANF,END,STEP)], color='blue')
plotE = list_plot([[n, EddE(n)] for n in range(ANF,END,STEP)], color='magenta')
show(plot2 + plot3 + plotE)

The plot below compares the number of exact decimal digits of the approximation (formula 3, blue curve) with the number of exact decimal digits guaranteed by this formula, 3*ln(3*n) (magenta curve). For comparison the red curve shows the exact decimal digits of formula 2. Formula 1 (given by DLMF/NIST) gives a poor approximation not shown here.

Restating exp(LogB3(n)) gives the following asymptotic expansion of the Bernoulli number B(n) for n>0 and even:

 $\mid {\text{B}}(n)\mid \ =4\pi \left({\frac {n}{2\pi e}}\right)^{n+1/2}\exp \left({\frac {1}{2}}+{\frac {{n}^{-1}}{12}}\,-{\frac {{n}^{-3}}{360}}\,+{\frac {{n}^{-5}}{1260}}\,+O\left({n}^{-7}\right)\right)$ Thanks to Charles R Greathouse IV for pointing out an error in a previous version.