This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/ComputationAndAsymptoticsOfBernoulliNumbers
From OeisWiki
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 = [1]; 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).
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 the comments on my homepage and my response to Don Knuth (pdf in German).)
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:
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.
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
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.
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 = [0]*(n//2+3); D[1] = 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] 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 = [0]*(n+2); D[1] = 1 R = [1]; 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 = [0]*(n+2); D[1] = 1 R = [1] ; 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[1]) 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 = [0]*(n//2+3); D[1] = 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[1]); 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] := 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) [1]
- Akiyama & Tanigawa (1999) [2]
- Luschny (2010) [3]
- Brent & Harvey (2011) [4]
The method of Knuth and Buckholtz (1967) [5] 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 = [0]*(n+2); D[1] = 1; b = True; R = [1]; 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 = [0]*(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[0]) ev = not ev return R
Luschny
def LuschnyBernoulli(n) : A = [0]*(2*n); A[1] = 1; R = [1]; 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 = [0]*(n+1); T[1] = 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 = [1]; 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
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 [6]).
- LogB1(n)=(1/2+n)*ln(n)+(1/2-n)*ln(pi)+(3/2-n)*ln(2)-n
- 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)))
- 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))))
where shown here.
Restating exp(LogB3(n)) gives the following asymptotic expansion of the Bernoulli number B(n) for n>0 and even:
Thanks to Charles R Greathouse IV for pointing out an error in a previous version.
References
- ↑, volume 7 (1877), 157-187.
- ↑ M. Kaneko, The Akiyama-Tanigawa algorithm for Bernoulli numbers, Journal of Integer Sequences, Vol. 3 (2000).
- ↑ P. Luschny, A sequence transformation and the Bernoulli numbers, OEIS wiki, Apr 2010.
- ↑ R. P. Brent and D. Harvey, Fast computation of Bernoulli, Tangent and Secant numbers, arXiv, Aug 2011.
- ↑ D. E. Knuth and T. J. Buckholtz, Computation of tangent, Euler and Bernoulli numbers. Math. Comp. 21 1967, 663-688.
- ↑ P. Luschny, An approximation of the Bernoulli numbers, 2007, updated 2012.