# Bernoulli numbers (Python)

**Other implementations**: Erlang |**Python**

The Bernoulli numbers *B*_{m} are a sequence of rational numbers defined by the identity

with the initial condition *B*_{0} = 1. Here, denotes a binomial coefficient. The sequence, which continues -1/2, 1/6, 0, -1/30, 0, 1/42, ..., obeys no simple pattern, and computing *B _{m}* for

*m*as small as a few hundred is a relatively time-consuming task.

The following programs for calculating Bernoulli numbers will assume the presence of a rational type `mpq` (with support for basic arithmetic) and a function `bincoef` for calculating binomial coefficients, neither of which is available in Python's standard library. The GMPY library provides these, however.

We will have to create a trivial wrapper around GMPY's `bincoef` since for some reason it does not handle the case when *k* < 0, which will be of use later on.

<<dependencies>>=fromgmpyimportmpq, bincoef as _bincoefdefbincoef(n, k):ifk < 0:return0else:return_bincoef(n, k)

## Contents |

## [edit] Basic implementation

The *m*:th Bernoulli number can be obtained by extracting the last term from the sum we used to define the sequence, but this requires that *B*_{1} ... *B*_{m-1} are already known. This is a perfect opportunity to employ recursion: `bernoulli(m)` will call `bernoulli(0) ... bernoulli(m-1)`, which in turn will call `bernoulli(0) ... bernoulli(m-2)` and so on, until all fall back on `bernoulli(0)`. The number of computations needed will grow exponentially with *m*, so a practical implementation must use memoization. We will create a decorator for this purpose:

<<memo>>=classmemo:def __init__(self, f):self.func = f self.cache ={}self.highest = 0def __call__(self, m):ifm < self.highest:returnself.cache[m]else: x = self.func(m) self.cache[m] = x self.highest = mreturnx

Now the main function. As a final optimization, we will exploit the fact that *B _{m}* = 0 if

*m*> 1 is odd.

<<bernoulli_simple.py>>=dependencies memo @memodefbernoulli(m):ifm == 0:return1ifm == 1:return-mpq(1,2)ifm % 2:return0 s = 0forjinrange(0, m): s += bincoef(m+1, j) * bernoulli(j) b = -mpq(s)/bincoef(m+1, m)returnb

The function can be used as follows:

>>> [bernoulli(n) for n in range(10)] [1, mpq(-1,2), mpq(1,6), 0, mpq(-1,30), 0, mpq(1,42), 0, mpq(-1,30), 0]

## [edit] A faster implementation

The computation can be performed more efficiently according to the following identities due to Srinivasa Ramanujan:

where

So instead of summing over *m* Bernoulli numbers and binomial coefficients, a sum over *m*/6 terms suffices.

We implement *A*:

<<A>>=defA(m, M): s = 0forjinrange(1, M+1): s += bincoef(m+3, m-6*j) * bernoulli(m-6*j)returns

We may now implement the three different cases modulo 6 as three separately memoized functions. (Separating them is necessary since the `memo` decorator assumes that if `bernoulli(m)` has been computed, then so have *all* of `bernoulli(0) ... bernoulli(m-1)`. But here we want to rely only on *m* that are congruent modulo 6.)

<<modular cases>>=@memodefb0mod6(m):return(mpq(m+3, 3) - A(m, m//6)) / bincoef(m+3, m) @memodefb2mod6(m):return(mpq(m+3, 3) - A(m, (m-2)//6)) / bincoef(m+3, m) @memodefb4mod6(m):return(-mpq(m+3, 6) - A(m, (m-4)//6)) / bincoef(m+3, m)

Wrapping things up:

<<main faster>>=defbernoulli(m):assertm >= 0ifm == 0:return1ifm == 1:return-mpq(1,2)ifm % 6 == 0:returnb0mod6(m)ifm % 6 == 2:returnb2mod6(m)ifm % 6 == 4:returnb4mod6(m)return0

<<bernoulli_faster.py>>=dependencies memo A modular cases main faster

## [edit] A further enhancement

There is some inefficiency in the implementation of *A*, namely that it calculates a new binomial coefficient "from scratch" for each term of the sum. For large inputs, it is faster to compute each binomial coefficient from the previous one.

<<A faster>>=defproduct(low, high): p = 1forkinrange(low, high+1): p *= kreturnpdefA(m, mtop): s = 0 a = bincoef(m+3, m-6)forjinrange(1, mtop+1): s += a*bernoulli(m-6*j) a *= product(m-6 - 6*j + 1, m-6*j) a //= product(6*j+4, 6*j+9)returns

<<bernoulli_faster_yet.py>>=dependencies memo A faster modular cases main faster

## [edit] Testing

A good way to check that all implementations are correct is to test them against each other. We will also time the computation of *B*_{1000} with each method.

<<test.py>>=frombernoulli_simpleimportbernoulli as bernoulli1frombernoulli_fasterimportbernoulli as bernoulli2frombernoulli_faster_yetimportbernoulli as bernoulli3fromtimeimportclockforiinrange(50):assertbernoulli1(i) == bernoulli2(i) == bernoulli3(i) a = clock() temp = bernoulli1(1000)

The output will probably resemble the following:

simple: 15.7621247952 seconds faster: 1.00002039365 seconds faster yet: 0.568878726207 seconds

The difference will become even more significant with larger input. In fact, the third version was used to compute the 10,000 first Bernoulli numbers in less than an hour on a 1400 MHz Athlon (version two needed four hours). By contrast, *Mathematica* took 1.5 hours to perform the same calculation and the specialized program calcbn needed 2.5 hours. However, it should be noted that the program implemented here is not well suited for calculating individual extremely large Bernoulli numbers, say, B_{100,000}.

Download code |