Bernoulli numbers (Python)

From LiteratePrograms

Other implementations: Erlang | Haskell | Python

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>>=
from gmpy import mpq, bincoef as _bincoef

def bincoef(n, k):
    if k < 0:
        return 0
    else:
        return _bincoef(n, k)

Contents 1 Basic implementation 2 A faster implementation 3 A further enhancement 4 Testing

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₁ ... 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>>=
class memo:
    def __init__(self, f):
        self.func = f
        self.cache = {}
        self.highest = 0
    def __call__(self, m):
        if m < self.highest:
            return self.cache[m]
        else:
            x = self.func(m)
            self.cache[m] = x
            self.highest = m
            return x

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

@memo
def bernoulli(m):
    if m == 0: return 1
    if m == 1: return -mpq(1,2)
    if m % 2: return 0
    s = 0
    for j in range(0, m):
        s += bincoef(m+1, j) * bernoulli(j)
    b = -mpq(s)/bincoef(m+1, m)
    return b

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]

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>>=
def A(m, M):
    s = 0
    for j in range(1, M+1):
        s += bincoef(m+3, m-6*j) * bernoulli(m-6*j)
    return s

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>>=
@memo
def b0mod6(m):
    return (mpq(m+3, 3) - A(m, m//6)) / bincoef(m+3, m)

@memo
def b2mod6(m):
    return (mpq(m+3, 3) - A(m, (m-2)//6)) / bincoef(m+3, m)

@memo
def b4mod6(m):
    return (-mpq(m+3, 6) - A(m, (m-4)//6)) / bincoef(m+3, m)

Wrapping things up:

<<main faster>>=
def bernoulli(m):
    assert m >= 0
    if m == 0: return 1
    if m == 1: return -mpq(1,2)
    if m % 6 == 0: return b0mod6(m)
    if m % 6 == 2: return b2mod6(m)
    if m % 6 == 4: return b4mod6(m)
    return 0

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

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>>=
def product(low, high):
    p = 1
    for k in range(low, high+1):
        p *= k
    return p

def A(m, mtop):
    s = 0
    a = bincoef(m+3, m-6)
    for j in range(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)
    return s

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

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₁₀₀₀ with each method.

<<test.py>>=
from bernoulli_simple import bernoulli as bernoulli1
from bernoulli_faster import bernoulli as bernoulli2
from bernoulli_faster_yet import bernoulli as bernoulli3
from time import clock

for i in range(50):
    assert bernoulli1(i) == bernoulli2(i) == bernoulli3(i)

a = clock()
temp = bernoulli1(1000)
print "simple:", clock() - a, "seconds"

a = clock()
temp = bernoulli2(1000)
print "faster:", clock() - a, "seconds"

a = clock()
temp = bernoulli3(1000)
print "faster yet:", clock() - a, "seconds"

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