Computing Bernoulli and Tangent Numbers
Richard P. Brent
MSI, ANU

17/18 May 2011

Dedicated to Jon Borweinon the occasion of his 60th birthday
Copyright c 2011, R. P. Brent

Richard Brent

Computing Bernoulli and Tangent Numbers

Summary
Bernoulli numbers are rational numbers Bn defined by thegenerating function
Xn≥0

Bnzzn
=.
n!
exp(z) − 1

They are of interest in number theory and are related to special
values of the Riemann zeta function. They also occur as
coefficients in the Euler-Maclaurin formula.
The closely related Tangent numbers Tn , and Secant numbers
Sn , defined byXn>0Tnz 2n−1= tan z,(2n − 1)!

Xn≥0Snz 2n= sec z,(2n)!

are more convenient for computation because they are integers.
Richard Brent

Computing Bernoulli and Tangent Numbers

Summary continued

In this talk I will consider some algorithms for computing
Bernoulli, Secant and Tangent numbers.
Recently, David Harvey [Math. Comp. 2010] showed that thesingle number Bn can be computed inO(n2 (log n)2+o(1) )bit-operations. In fact, the Bernoulli numbers B0 , . . . , Bn can allbe computed with the same complexity bound (and similarly for
Secant and Tangent numbers).

Richard Brent

Computing Bernoulli and Tangent Numbers

Summary continued

I will first give a relatively simple algorithm that achieves theslightly weaker bound O(n2 (log n)3+o(1) ). If time permits, I will
sketch the improvement to O(n2 (log n)2+o(1) ) bit-operations.
I will also give very simple in-place algorithms for computing thefirst n Secant or Tangent numbers using O(n2 ) integer
operations. Although they are not the asymptotically fastest
algorithms, they are extremely simple and convenient formoderate values of n.

Richard Brent

Computing Bernoulli and Tangent Numbers

Advertisement

Much of the material for this talk is drawn from my recent book:
Richard P. Brent and Paul Zimmermann,
Modern Computer Arithmetic,
Cambridge University Press, 2010, 237 pp.
(online version available from my website).
In particular, see §4.7.2 and exercises 4.35–4.41.
Equation numbers such as “(4.58)” are as in the book.

Richard Brent

Computing Bernoulli and Tangent Numbers

Bernoulli numbers
From the generating function
XBnn≥0zzn=n!
exp(z) − 1

it’s easy to see that the Bn are rational numbers,
B2n+1 = 0 if n > 0, and they satisfy the recurrence
B0 = 1,

k Xk +1
jj=0

Bj = 0 for k > 0.

(4.58)

It’s sometimes convenient to consider scaled Bernoulli numbers
Cn = B2n /(2n)!, with generating function
XCn z 2n =n≥0Richard Brentz/2.
tanh(z/2)
Computing Bernoulli and Tangent Numbers

The Von Staudt – Clausen theorem
Computing a few Bn , we findB0 = 1,B6 = 1/42,
B14 = 7/6,B1 = −1/2, B2 = 1/6,
B8 = −1/30, B10 = 5/66,etc.

B4 = −1/30,
B12 = −691/2730,

What are the denominators? This is answered by theVon Staudt – Clausen Theorem (1840), which says that
0B2n:= B2n +X(p−1)|2n1∈ Z,p(where the sum is over primes p for which p − 1 divides 2n).
0Thus, in a program it might be more convenient to store B2nthan B2n .
Richard Brent

Computing Bernoulli and Tangent Numbers

Connection with the Riemann zeta-function
Euler found that the Riemann zeta-function for evennon-negative integer arguments can be expressed in terms ofBernoulli numbers – the relation is(−1)k −1

2ζ(2k )
B2k=.
(2k )!
(2π)2kSince ζ(2k ) = 1 + O(4−k ) as k → +∞, we see that
|B2k | ∼2 (2k )!
.
(2π)2kFrom Stirling’s approximation to (2k )! we see that the numberof bits in the integer part of B2k is 2k lg k + O(k ).
Thus, it takes Ω(n2 log n) space to store B1 , . . . , Bn .

Richard Brent

Computing Bernoulli and Tangent Numbers

Another connection with the zeta-function

From the functional equation for the Riemann zeta-function, wealso haveBn+1
ζ(−n) = −n+1for n ∈ Z, n ≥ 1.

Richard Brent

Computing Bernoulli and Tangent Numbers

Computing Bernoulli numbers
From the generating function z/(exp(z) − 1) we obtain therecurrence
B0 = 1,

k 
Xk +1
j=0j

Bj = 0 for k > 0.

(4.58)

This recurrence has traditionally been used to compute
B0 , . . . , B2k with O(k 2 ) arithmetic operations.
This is unsatisfactory if floating-point numbers are used,because the recurrence is numerically unstable: the relative
error in the computed B2k is of order 4k 2−n if the floating-point
arithmetic has a precision of n bits.

Richard Brent

Computing Bernoulli and Tangent Numbers

A Numerically Stable Recurrence
As before, let Ck = B2k /(2k )!. Thensinh(z/2) XCk z 2k = cosh(z/2),
z/2k ≥0and equating coefficients gives the recurrence
kX
j=0Cj
1=.
k−j(2k + 1 − 2j)! 4(2k )! 4k(4.60)

Using this recurrence to evaluate C0 , C1 , . . . , Ck , the relative
error in the computed Ck is only O(k 2 2−n ), which is satisfactoryfrom a numerical point of view.

Richard Brent

Computing Bernoulli and Tangent Numbers

An Asymptotically Fast Algorithm

Harvey (2010) showed how Bn could be computed exactly,
using a modular algorithm and the Chinese remainder theorem,in O(n2 (log n)2+ε ) bit-operations. (We write ε for terms that areo(1) as n → ∞.)We’ll show how to compute all of B0 , . . . , Bn with almost thesame complexity bound (only larger by a factor O(log n)).

Richard Brent

Computing Bernoulli and Tangent Numbers

Digression – Reciprocals of Power Series

Let A(z) = a0 + a1 z + a2 z 2 + · · · be a power series withcoefficients in some field F (e.g. F = Q or R), with a0 6= 0.
Let B(z) = b0 + b1 z + · · · be the reciprocal power series, soA(z)B(z) = 1.
Suppose we can multiply polynomials of degree n − 1 in F [z]using M(n) = O(n(log n)1+ε ) field operations. Then, usingNewton’s method [Kung and Sieveking], we can computeb0 , . . . , bn−1 with the same complexity O(n(log n)1+ε ), up to aconstant factor.

Richard Brent

Computing Bernoulli and Tangent Numbers

Application – an Asymptotically Fast Algorithm
TakingA(z) = (exp(z) − 1)/zand working with (n lg(n) + O(n))-bit floating-point numbers, weget B0 , . . . , Bn to sufficient accuracy to deduce the exact(rational) result using O(n(log n)1+ε ) floating-point operations,each of which can be done withO(n(log n)2 log log n)
bit-operations by Schönhage–Strassen. Thus, overall we get
B0 , . . . , Bn with
O(n2 (log n)3+ε )bit-operations. Similarly for Secant and Tangent numbers.

Richard Brent

Computing Bernoulli and Tangent Numbers

Tangent and Secant numbers
The Tangent numbers Tn (n > 0) (also called Zag numbers)
are defined by
X

Tnn>0

sin zz 2n−1
= tan z =.
(2n − 1)!
cos z

Similarly, the Secant numbers Sn (n ≥ 0) (also called Euler orZig numbers) are defined by
Xn≥0Snz 2n1= sec z =.
(2n)!
cos z

Unlike the Bernoulli numbers, the Tangent and Secant numbers
are positive integers.

Richard Brent

Computing Bernoulli and Tangent Numbers

Asymptotics
Because tan z and sec z have poles at z = π/2, we expect
Tn to grow roughly like (2n − 1)!(2/π)n and Sn like (2n)!(2/π)n .
More precisely, let
ζ0 (s) = (1 − 2−s )ζ(s) = 1 + 3−s + 5−s + · · ·be the odd zeta-function.ThenTk
22k +1 ζ0 (2k )22k +1=∼.
(2k − 1)!
π 2kπ 2k
We also haveSk22k +2∼ 2k +1 .
(2k )!
πRichard Brent

Computing Bernoulli and Tangent Numbers

Bernoulli numbers via Tangent numbers
From the formulas for Tk and B2k in terms of ζ(2k ), we see that
Tk = (−1)k −1 22k (22k − 1)

B2k2k

(this can be proved directly, without involving the zeta-function).
Since Tk ∈ Z, the odd primes in the denominator of B2k must
divide 22k − 1. This is compatible with the Von Staudt–Clausen
theorem, since (p − 1)|2k =⇒ p|(22k − 1) by Fermat’s littletheorem.
Tk has about 4k more bits than dB2k e, but both have
2k lg k + O(k ) bits, so asymptotically not much difference.
Thus, to compute B2k we don’t lose much by first computing Tk ,and this may be more convenient since Tk ∈ Z, B2k ∈ Q.
Richard Brent

Computing Bernoulli and Tangent Numbers

Getting Rid of a “log” Factor in the Time Bound
To improve the algorithm for Bernoulli numbers, we use theKronecker–Schönhage trick. Here is an outline.
Fix n > 1, choose p = dn lg(n)e, N = 2np, z = 2−p .
Write down N-bit approximations to (2n)! sin(z) and(2n)! cos(z) from the truncated Taylor series.
Perform an N-bit division (using Newton’s method) to get an
N-bit approximation to tan(z) in time O(N log(N) log log(N)).
Multiply by (2n − 1)! and read off the integers
Tk0 = Tk (2n − 1)!/(2k − 1)! from the binary representation.
Now deduce the Tk and B2k , k ≤ n. The overhead involving
factorials can be handled within the overall time bound, which isO(N log(N) log log(N)) = O(n2 (log n)2 log log n).

Richard Brent

Computing Bernoulli and Tangent Numbers

Algorithms based on 3-term Recurrences
Akiyama and Tanigawa gave an algorithm for computing
Bernoulli numbers based on a 3-term recurrence. However, it is
only useful for exact computations, since it is numericallyunstable if applied using floating-point arithmetic.
We’ll give a stable 3-term recurrence and corresponding
in-place algorithm for computing Tangent numbers. It isperfectly stable since all operations are on positive integers andthere is no cancellation. Also, it involves less arithmetic thanthe Akiyama-Tanigawa algorithm.
The three-term recurrence was given by Buckholtz and Knuth(1967), but they did not give the in-place algorithm explicitly.
There is a similar 3-term recurrence and stable in-placealgorithm for computing Secant numbers.

Richard Brent

Computing Bernoulli and Tangent Numbers

A 3-term Recurrence for Computing Tangent Numbers
Write t = tan x, D = d/dx, so Dt = 1 + t 2 andD(t n ) = nt n−1 (1 + t 2 ) for n ≥ 1.
It is clear that D n t is a polynomial in t, say Pn (t).
P
Write Pn (t) = j≥0 pn,j t j . Then deg(Pn ) = n + 1 and,from the formula for D(t n ),pn,j = (j − 1)pn−1,j−1 + (j + 1)pn−1,j+1 .

(4.63)We are interested in Tk = p2k −1,0 , which can be computed fromthe 3-term recurrence in O(k 2 ) operations using the obvious
boundary conditions.
We can save work by noticing that pn,j = 0 if n + j is even.

Richard Brent

Computing Bernoulli and Tangent Numbers

Algorithm TangentNumbers
Input: positive integer n
Output: Tangent numbers T1 , . . . , Tn
T1 ← 1for k from 2 to n
Tk ← (k − 1)Tk −1
for k from 2 to n
for j from k to n
Tj ← (j − k )Tj−1 + (j − k + 2)Tj
return T1 , T2 , . . . , Tn .
The first for loop initializes Tk = pk −1,k = (k − 1)!. The variable
Tk is then used to store pk ,k −1 , pk +1,k −2 , . . ., p2k −2,1 , p2k −1,0 at
successive iterations of the second for loop. Thus, when thealgorithm terminates, Tk = p2k −1,0 (correct).
Richard Brent

Computing Bernoulli and Tangent Numbers

Illustration
The process in the case n = 3 is illustrated in the following(m)diagram, where Tk denotes the value of the variable Tk at
successive iterations m = 1, 2, . . . , n:(1)

T1

= p0,1.

&(1)
T1 = p1,0

(1)

T2
&

= p1,2.

&(2)
T2 = p2,1.
(2)

T2

(1)

T3
&
(2)

= p3,0T3&= p2,3.
= p3,2.
(3)
T3 = p4,1.
(3)
T3 = p5,0Richard BrentComputing Bernoulli and Tangent Numbers

Complexity of Algorithm TangentNumbers

Algorithm TangentNumbers takes Θ(n2 ) operations on positiveintegers. The integers Tn have O(n log n) bits, other integers
have O(log n) bits.
Thus the overall complexity is O(n3 (log n)1+ε ) bit-operations,or O(n3 log n) word-operations if n fits in a single word.
The algorithm is not optimal, but it is good in practice formoderate values of n, and much simpler than asymptoticallyfaster algorithms.

Richard Brent

Computing Bernoulli and Tangent NumbersAcknowledgements
I Jon Borwein, Kurt Mahler and George Szekeres forencouraging my belief that high-precision computations
are useful in “experimental” mathematics.
I David Harvey for discussions on the complexity ofalgorithms for computing Bernoulli numbers.
I Donald Knuth and Thomas Buckholtz for giving thethree-term recurrence (4.63) for Tangent numbers.
I Ben F. “Tex” Logan, Jr., for suggesting the use of Tangent

numbers to compute Bernoulli numbers.
I Christian Reinsch for pointing out the numerical instabilityof the recurrence (4.58) and suggesting the use of thenumerically stable recurrence (4.60).
I Paul Zimmermann for coauthoring our book Modern

Computer Arithmetic.
Richard Brent

Computing Bernoulli and Tangent Numbers

References
1. Jonathan M. Borwein and David H. Bailey, Mathematics byExperiment: Plausible Reasoning in the 21st Century, second
edition, A. K. Peters, 2008.
2. Richard P. Brent, Unrestricted algorithms for elementary andspecial functions, in Information Processing 80, North-Holland,
Amsterdam, 1980, 613–619. arXiv:1004.3621v1
3. Richard P. Brent and Paul Zimmermann, Modern Computer
Arithmetic, Cambridge University Press, 2010, 237 pp.
4. Ronald L. Graham, Donald E. Knuth, and Oren Patashnik,Concrete Mathematics, third edition, Addison-Wesley, 1994.
5. David Harvey, A multimodular algorithm for computing
Bernoulli numbers, Mathematics of Computation 79 (2010),2361–2370.

Richard Brent

Computing Bernoulli and Tangent Numbers

References continued

6. Masanobu Kaneko, The Akiyama-Tanigawa algorithm for
Bernoulli numbers, Journal of Integer Sequences 3 (2000).
Article 00.2.9, 6 pages.
7. Donald E. Knuth, Euler’s constant to 1271 places,
Mathematics of Computation, 16 (1962), 275–281.
8. Donald E. Knuth and Thomas J. Buckholtz, Computation ofTangent, Euler, and Bernoulli numbers, Mathematics ofComputation 21 (1967), 663–688.
9. Christian Reinsch, personal communication, about 1979.

Richard Brent

Computing Bernoulli and Tangent Numbers