This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/TheLostBernoulliNumbers

From OeisWiki

Jump to: navigation, search

The Euler-Bernoulli Diamond
the Lost Bernoulli Numbers.

KEYWORDS: Euler zeta numbers, Euler numbers, Bernoulli numbers, Euler tangent numbers, Bernoulli secant numbers, generalized Bernoulli numbers.

Concerned with sequences:

A000111, A000364, A009006, A009843, A099612, A059841, A099617,
A176289, A193472, A193473, A193475, A193476.



Corresponding to the three cases sec, tan, sec + tan there are three kinds of Euler numbers: A000364 (Euler (or secant or "Zig") numbers), A000182 (tangent (or "Zag") numbers), A000111 (Euler or up/down numbers).

Each kind of Euler numbers corresponds to a kind of Bernoulli numbers: The Bernoulli secant numbers, Bernoulli tangent numbers and the Bernoulli secant/tangent numbers.

The classical Euler numbers (HMF, Maple, Mathematica) are the Euler secant numbers. The Euler tangent numbers are also present on OEIS. The classical Bernoulli numbers (HMF, Maple, Mathematica) are the Bernoulli tangent numbers. The Bernoulli secant numbers are missing on OEIS -- until now.

Our development is centered on the Euler Zeta numbers studied by Leonhard Euler in 1735. We give a small collection of formulas, old and new ones, representing these milestone numbers (A099612/A099617).

In passing we will implement the algorithms for fast computation of Bernoulli, Tangent and Secant Numbers as given by Richard P. Brent and David Harvey this month (August 2011).

In the last part we will generalize the numbers to functions using the generalized (Hurwitz) zeta function. These beautiful simple functions seem to be little known. They provide the shortest and most systematic road to the Bernoulli-Euler family of numbers.

The Euler−Bernoulli Diamond

We will discuss the Euler−Bernoulli diamond, a graph showing the various relations between the Bernoulli and the Euler numbers. All sequences start at n = 0. (For better readability right click the picture to open it in a separate window; then (double) click again on the picture in the new window.)


The Bernoulli and Euler numbers are defined in the Digital Library of Mathematical Functions. We will make two small deviations (see the appendix for some motivation).

  • We look at the unsigned numbers.
  • We set the Bernoulli number B1 to 0.

With the use of the generating function of the secant and the tangent functions we can define all the seven functions in the Euler−Bernoulli Diamond.

Name Definition Seq
Euler Zeta EZ := n -> gf(sec + tan, n): Fraction
Euler Generalized EG := n -> `if`(n=0,1,EZ(n)n!): A000111
Euler Tangent ET := n -> `if`(n mod 2 = 0,0,EG(n)): A009006
Euler Classic EC := n -> `if`(n mod 2 = 1,0,EG(n)): A000364
Bernoulli Generalized BG := n -> `if`(n=0,1,EZ(n-1)n!/(4^n-2^n)): Fraction
Bernoulli Classic BC := n -> `if`(n mod 2 = 1,0,BG(n)): Fraction
Bernoulli Secant BS := n -> `if`(n mod 2 = 0,0,BG(n)): Fraction

As can be seen all definitions are based on a single one: the Euler zeta numbers, explored by Leonhard Euler in 1735 in his paper ‘De summis serierum reciprocarum’.

The (generalized) Euler numbers are nothing else than the Euler zeta numbers times n!. Looking at the even or at the odd indexed numbers gives the Euler tangent numbers and the classic Euler numbers respectively.

In a similar way the Bernoulli numbers are obtained. The generalized Bernoulli numbers are the Euler zeta numbers shifted by one times n!/(4n − 2n). Again looking at the even and the odd indexed numbers gives the classical and the secant Bernoulli numbers, respectively.

Clearly the point wise addition of the sequence of the classic Euler numbers and the Euler tangent numbers gives the generalized Euler numbers and the point wise addition of the sequence of the classic Bernoulli numbers and the Bernoulli secant numbers gives the generalized Bernoulli numbers.

  0 1 2 3 4 5 6 7 8
BG 1 1/2 1/6 3/56 1/30 25/992 1/42 427/16256 1/30
BC 1 0 1/6 0 1/30 0 1/42 0 1/30
BS 0 1/2 0 3/56 0 25/992 0 427/16256 0

The Euler zeta numbers and all flavors of the Bernoulli numbers are fractions. Thus we will also look the numerator and the denominator of these numbers separately.

Name Definition Seq
Euler Zeta Numerator EZN := n -> numer(EZ(n)): A099612
Euler Zeta Denominator EZD := n -> denom(EZ(n)): A099617
Bernoulli Generalized Numerator BGN := n -> numer(BG(n)): missing, now
Bernoulli Generalized Denominator BGD := n -> denom(BG(n)): missing, now
Bernoulli Classic Numerator BCN := n -> numer(BC(n)): A059841
Bernoulli Classic Denominator BCD := n -> denom(BC(n)): A176289
Bernoulli Secant Numerator BSN := n -> numer(BS(n)): A009843
Bernoulli Secant Denominator BSD := n -> denom(BS(n)): missing, now

And now the surprise: neither the Bernoulli secant numbers nor the generalized Bernoulli numbers are in the database! Surprise because the database is full of variants of the Bernoulli numbers, several hundreds of them at least. The name lost Bernoulli numbers refers to the Bernoulli secant numbers, the counterpart to the Euler tangent numbers.


Assume n ≥ 1. [ n even] and [n odd] are Iverson brackets. T(n) = (4^n−2^n)/n.

EZ EZ(n)n! [n odd] EZ(n) n! [n even] EZ(n) n! EZ(n-1) (n-1)! /T(n) [n even] EZ(n-1) (n-1)! /T(n) [n odd] EZ(n-1) (n-1)! /T(n)
EG EG(n) /n! [n odd] EG(n) [n even] EG(n) EG(n-1) /T(n) [n even] EG(n−1) /T(n) [n odd] EG(n−1) /T(n)
BG BG(n+1) T(n+1) /n! BG(n+1) T(n+1) [n odd] BG(n+1) T(n+1) [n even] BG(n+1) T(n+1) [n even] BG(n) [n odd] BG(n)

Formulas for Euler's Zeta numbers.

The cornerstones of our setup were the Euler zeta numbers. Together with the extremely simple transformations all sequences can be computed from these numbers.

  • By generating function:
gf := (f, n) -> coeff(series(f(x),x,n+2),x,n):
EZ := n -> gf(sec + tan, n):
  • By recursion:
S := proc(n, k) option remember; 
if k=0 then `if`(n=0,1,0) else S(n,k-1)+S(n-1,n-k) fi end:
EZ := n -> S(n,n)/n!;

This computation can be best understood by considering the triangle of the Euler-Bernoulli numbers A008281. The triangle can be constructed by using only additions.

0 1                
1 0 1              
2 0 1 1            
3 0 1 2 2          
4 0 2 4 5 5        
5 0 5 10 14 16 16      
6 0 16 32 46 56 61 61    
7 0 61 122 178 224 256 272 272  
8 0 272 544 800 1024 1202 1324 1385 1385
n\k 0 1 2 3 4 5 6 7 8
  • Using Elkies' formula (On the sums Sum_(k=-infinity...infinity) (4k+1)^(-n), Amer. Math. Monthly, 110(7) 2003, 561–573):
elksum := n -> 2*(Pi/2)^(-n-1)*sum((4*k+1)^(-n-1),k =-infinity..infinity):
EZ := n -> if n = 0 then 1 else elksum(n) fi:
  • In reference to a much more general view using Lis(z):
EZ := n -> if n = 0 then 1 else 2*I^(n+1)*polylog(-n,-I) / n! fi:
  • With the generalized zeta function:
EZ := s -> if s = 0 then 1 else abs(((2^s-4^s)*Zeta(0,-s+1,1)+

\text{EZ}(s)\Gamma(s)\ =\ \mid (2^s-4^s)\zeta(1-s,1)+\frac12 
4^s\left(\zeta(1-s,\frac14)-\zeta(1-s,\frac34)\right) \mid

P := proc(n, x) local k, j; add(add((-1)^j*2^(-k)*binomial(k,j)
*(k-2*j)^n*x^(n-k), j=0..k), k=0..n) end:
EZ := n -> if n = 0 then 1 else P(n-1, -I) / n! fi:

It is quite instructive to have a look how the EZ numbers are computed by the A193474 polynomials: as the row sums of the triangle below divided by n!.

0 1      
1 1      
2 1      
3 2      
4 6 -1    
5 24 -8    
6 120 -60 1  
7 720 −480 32  
8 5040 -4200 546 -1
9 40320 −40320 8064 -128
n\k 0 2 4 6

\operatorname{EZ}(n) = \frac{1}{n!}\sum_{k=0}^{n-1}\frac{(-1)^{ \lfloor 
(n-1-k)/2 \rfloor }}{2^{k}} \sum_{j=0}^{k} (-1)^j \binom{k}{j}(k-2j)^{n-1}

A combinatorial interpretation of the coefficients of the A193474 polynomials (the unsigned entries of the triangle above) would be very interesting. Clearly the inclusion-exclusion principle is here working behind the scenes.

[Edit, Oct 7] The hint comes from looking at the row sums of the coefficients of the polynomials: this is the number of labeled ordered partitions of an n-set into odd parts A006154. Thus the coefficients are the number of labeled ordered partitions of an n-set into n-k parts of odd size. Peter Bala pointed this out in A196776 (reversing the enumeration of the rows of A193474). Unfortunately he missed to include the empty set; thus his enumeration does not fit with the enumeration of A006154 nor with the enumeration of the not-ordered partitions A136630. [Edit end]

  • Using the Euler polynomials:
ep := n -> `if`(n=0, 1,2^n*(euler(n,1/2) - euler(n,1))*(-1)^iquo(n+1,2));

n! * EZ(n) = ep(n). (Note that this involves rational arithmetics.)

  • Using the Swiss Knife polynomials:

The definition of the Swiss Knife polynomials skp(n,x) will be given below.

n! * EZ(n) = (-1)^floor(n/2) * skp(n, n mod 2) . (Note that this involves only integer arithmetics.)

The connection with Pi.

There is a fascinating connection between Euler's Holy Grail, the Euler zeta numbers, and the number Π. The proposition is:

 \lim_{n \rightarrow \infty} \frac{2 \operatorname{EZ}(n-1)}{\operatorname{EZ}(n)}\ 
=\ \pi

Moreover, two successive terms always enclose the true value of Π.

 \scriptstyle 2, 4, 3, \frac{16}{5}, \frac{25}{8}, \frac{192}{61}, \frac{427}{136}, 
\frac{4352}{1385}, \frac{12465}{3968}, \frac{158720}{50521},\ldots

Compare 132049 and 132050 which start for some unaccountable reason at n = 3. More on this can be found in: J.-P. Delahaye, Le fascinant nombre Pi, Pour la Science, Paris, 1997. German translation: Pi - die Story, Birkhaeuser, 1999 Basel, p. 31. 

A look-alike for the denominators of the Bernoulli secant numbers.

Johannes W. Meijer pointed out the phenomenon of a look-alike of the denominators in the Taylor series for tan(x) (see A156769). Analogically with this we also observe a look-alike for the denominators of the Bernoulli secant numbers. Looking only at the odd case BSD(2*n+1) we find


On the other hand A193475(n) = 4*16^n - 2*4^n gives


and differs at n = 10, 31, 52, 73, 77, 94,... The fraction A193475(n)/BSD(2*n+1) leads to


I have a conjecture with regard to BSD(n):

BSD(2*n+1) = (4*16^n - 2*4^n) / gcd(2*n+1, 4*16^n - 2*4^n).

Are there any takers? I think this conjecture is easy to prove.

Efficient computation

The Euler and Bernoulli numbers can be efficiently computed with simple algorithms up to n ~ 1000. Here we show these algorithms implemented with Maple. Note that we compute lists; this means that the function argument demands that all numbers [0,..,n] are to be computed; however only the non-zero numbers are computed.

Small benchmarks showed that for instance the list of the (even) Bernoulli numbers between 0 and 1000 are computed twice as fast with these scripts than with the build in `bernoulli´ function of Maple VR5.

4^n−2^n and (16^n−4^n)/2 and 4*16^n−2*4^n

A020522 := proc(n) option remember;
if n = 0 then 0 elif n = 1 then 2 else
6*A020522(n-1)-8*A020522(n-2) fi end:

A026337 := proc(n) option remember;
if n = 0 then 0 elif n = 1 then 6 else
20*A026337(n-1)-64*A026337(n-2) fi end:

A193475 := proc(n) option remember;
if n = 0 then 2 elif n = 1 then 56 else
20*A193475(n-1)-64*A193475(n-2) fi end:

Euler classic (even n, signless), Secant numbers.

EulerSecant_list := proc(n) local S,k,j; 
S := array(0..n); S[0] := 1;
for k from 1 to n do
    S[k] := k*S[k-1] od;
for k from 1 to n do
    for j from k to n do
        S[j] := (j-k)*S[j-1]+(j-k+1)*S[j] od od;
convert(S,list) end:

Euler tangent (odd n, signless), Tangent numbers.

EulerTangent_list := proc(n) local T,k,j; 
T := array(1..n); T[1] := 1;
for k from 2 to n do
    T[k] := (k-1)*T[k-1] od;
for k from 2 to n do
    for j from k to n do
       T[j] := (j-k)*T[j-1]+(j-k+2)*T[j] od od;
convert(T,list) end:

Bernoulli classic (even n, signless).

BernoulliTangent_list := proc(n) local T,k;
T := EulerTangent_list(n);

Bernoulli secant (odd n, signless).

BernoulliSecant_list := proc(n) local S,k;
S := EulerSecant_list(n);
[seq(S[k]*(2*k-1)/A193475(k-1),k=1..n)] end:

Euler generalized (signless).

n -> [a(0), a(1), ..., a(n), a(n+[n even])]
Nice property: S and T could be computed parallel!

EulerGeneralized_list := proc(n) local S,T,L;
L := iquo(n+1+irem(n+1,2),2);
S := EulerSecant_list(L);
T := EulerTangent_list(L);
map(op, %) end:

Bernoulli generalized (signless).

n -> [a(0), a(1), ..., a(n)]

BernoulliGeneralized_list := proc(n) local k;
[1,seq(%[k]*k/A020522(k),k=1..n)] end:

Jakob's and Jan-Pierre's good old `archaic´ view

If one prefers the signed version of the numbers one has only to replace a single expression in the above formulas: Instead of setting EZ := n -> gf(sec + tan, n) define EZ := n -> gf(sech + tanh, n).

Basically our suggestion follows the practice of two eminent mathematicians: Jakob Bernoulli himself (Ars Conjectandi, 1713) and Jean-Pierre Serre (A Course in Arithmetic). Both look only at the even Bernoulli numbers starting with 1/6 (i.e. at n = 2 in HMF definition). (On Mathworld this view on the Bernoulli numbers is called the `archaic´ view, doing a little wrong to the first winner of the Abel Prize in modern times.) Thus our proposal just adds the case n = 0 to their setup.

The Bernoulli zeta functions

In this section we will look at the Bernoulli functions which interpolate the Bernoulli numbers in their signed form for n ≥ 2. Both the simplicity of the definition as well as the importance of the generalized zeta function (or Hurwitz zeta function) makes this my favorite route to the Bernoulli numbers.

 \operatorname{ZetaBS}(z) \ = \ \frac{\zeta\left(z,\frac14\right)-\zeta\left(z,\frac34\right)}{2^z 
- 2}

 \operatorname{ZetaBG}(z) \ = \ \zeta\left(z\right) + \operatorname{ZS}(z)

ZetaBS := z -> (Zeta(0,z,1/4) - Zeta(0,z,3/4))/(2^z-2):
ZetaBG := z -> Zeta(z) + ZetaBS(z):

# Generalized Bernoulli function
BGF := z -> -z*ZetaBG(1-z):

# Tangent (classical) Bernoulli function
BTF := z -> -z*Zeta(1-z):

# Secant (lost) Bernoulli function
BSF := z -> -z*ZetaBS(1-z):   

The Euler zeta functions

In this section we will look at the Euler functions which interpolate the Euler numbers in their signed form for n ≥ 1. The definition is again based on their representation via the Hurwitz zeta function. Things are in complete analogy to the Bernoulli case.

 \operatorname{ZetaES}(z) \ = \ \zeta \left(z,\frac14\right) - 

 \operatorname{ZetaET}(z) \ = \ \left(2^z - 2\right) \zeta\left(z\right)

ZetaES := z -> Zeta(0,z,1/4) - Zeta(0,z,3/4):
ZetaET := z -> (2^z-2)*Zeta(z):

# Generalized Euler function
EGF := z -> 2*4^z*(ZetaES(-z) + ZetaET(-z)):

# Secant (classical) Euler function
ESF := z -> 2*4^z*ZetaES(-z):

# Tangent Euler function
ETF := z -> 2*4^z*ZetaET(-z):

The functions in the plot are all scaled by z!. Note how both the classic Euler numbers (blue points) as well as the Euler tangent numbers (red points) lie on the graph of the generalized Euler function (green).

In fact, the points on the green line are exactly the signed Euler Zeta numbers, the absolute values of which constitutes the core sequence of our Euler-Bernoulli diamond graph. So back to the roots, anno 1735 Leonhardo scribit:

The connection with the Swiss-Knife polynomials

As every child knows, the Bernoulli and the Euler numbers are just special cases of the Bernoulli and Euler polynomials. So it would be nice to have also polynomials for all the other members of the Euler−Bernoulli family of numbers, wouldn't it? Well, so we need now the Bernoulli secant polynomials, the Euler zeta polynomials, the Euler tangent polynomials etc. Maybe I should better write this up next month in this blog? Fortunately this is not necessary. For two reasons.

Firstly we do not need Bernoulli secant polynomials, Euler zeta polynomials or Euler tangent polynomials. We even do not need the Bernoulli and the Euler polynomials! What we really want is a single sequence of polynomials which covers all these cases. At least with regard to the numbers involved. Fortunately such polynomials exist. And they are much nicer than the Bernoulli and the Euler polynomials because they have integer coefficients instead of fractions like the E−B polynomials.

And secondly I do not need to write another blog; I have already written a blog on these polynomials last year.

x2 − 1
x3 − 3 x
x4 − 6 x2 + 5
x5 − 10 x3 + 25 x
x6 − 15 x4 + 75 x2 − 61
x7 − 21 x5 + 175 x3 − 427 x
x8 − 28 x6 + 350 x4 − 1708 x2 + 1385

These are Swiss-Knife Polynomials. Let us denote these polynomials by σ(n, x). Then we have:

Euler zeta evenσ(n, 0) / n!
Euler zeta oddσ(n, 1) / n!
Euler zeta(-1)^floor(n/2) σ(n, n mod 2) / n!
Euler classic (secant)σ(n, 0)
Euler tangentσ(n, 1)
Euler generalized(-1)^floor(n/2) σ(n, n mod 2)
Bernoulli classic (tangent)σ(n − 1, 1) n / (4n − 2n)
Bernoulli secantσ(n − 1, 0) n / (4n − 2n)
Bernoulli generalizedσ(n − 1, (n−1) mod 2) n / (4n − 2n)
Numerator Bernoulli secantσ(n, 0) (n+1)
Genocchiσ(n, −1) (n+1) / 2n
Springerσ(n, 1/2) 2n

This table, together with the representation given in my May 2010 blog gives explicit formulas for all the numbers we have considered here. It suffices to set

 \sigma_n(x) = \sum_{k=0}^{n} \frac{\alpha(k+1)}{2^{ \lfloor k/2 
\rfloor}} \sum_{v=0}^{k}(-1)^{v} \binom{k}{v}(x+v+1)^n

where α(k) is the sequence 0,1,1,1,0,-1,-1,-1,0,1,1,1,0,-1,-1,-1,... (A046978).

sigma := proc(n, x) local v, k, pow, A046978;
pow := (a,b) -> if a = 0 and b = 0 then 1 else a^b fi;
A046978 := k -> `if`(k mod 4 = 0,0,(-1)^iquo(k,4));
add((-1)^v*binomial(k,v)*pow(v+x+1,n),v=0..k),k=0..n) end:


Download as pdf files: 

Generalized Bernoulli numbers and Zeta functions,
Generalized Euler numbers and Zeta functions.

Personal tools