In this post, I would like to demonstrate how to calculate the Bernoulli numbers. The Bernoulli numbers are a sequence of rational numbers and they’re usually denoted . If you want a list of their values, you can go to the OEIS … oh, wait, you can’t actually look them up directly in the OEIS because they are not integers. But, hey, since they are rational numbers, you can just look up their numerators and denominators separately 😆.

The reason why I am interested in them is that they can be used to evaluate the Riemann Zeta function, as is explained in this paper. This great mathologer video gives another reason why one might be interested in the Bernoulli numbers and also explains how Bernoulli found them in the first place. This will help you understand this post’s featured image. Eventually, the mathologer video leads up to the Euler-Maclaurin formula, which is also the underlying reason why the Bernoulli numbers are useful for evaluating the Riemann Zeta function.

In the following we will aside the genesis and applications of the Bernoulli numbers and we will not discuss the convergence of the series through which we define them. We will assume that this has been taken care of. This post is about the numerical calculation of the Bernoulli numbers.

1. The arithmetic

The Bernoulli numbers are defined via the power series of the following fraction:

To get a handle on the Bernoulli numbers, which are the parts of the coefficients denoted by , it helps to observe the following expansion of the inverse fraction first:

Here we have only used the well-known series expansion of the exponential function and then multiplied out the factor and re-written the result as a sum.

No, if we multiply both fractions, we get the following relation:

and we just imagine that this product has its own power series expansion . We know that this series has to be equal to 1 regardless of the value of . And it is not hard to find out that they obey the following relation:

To obtain this expression of the we have used what is commonly known as the Cauchy product formula. But this just means that if we have to sums and and their prodcut expansion then we just “collect” all terms from the original two series where . For example:

Here, we have two possible combinations of a-b-factors in front of the linear term(s). Therefore, the coefficient in the prodcut’s expansion needs to sum over those. And in the infinite case we can ignore the fact that our toy example has a larger highes term in .

Ok, back to our main result:

Here, we can easily check (if we remember that ) that . Since the total product is identically equal to 1 and is already equal to 1, we know that .

2. The triangle

Ok, so far we have managed to get the first Bernoulli number . But given the fact that all must be zero, we can write down the equations for the next few Bernoulli numbers. Let’s do this and then discuss how to solve that:

…

Now, if we stare at this for a while, we might actually realize that the coefficients of the in this list of equations are actually the entries in Pascal’s triangle. That is, we are laying out all binomial coefficients where the “numerator” of the binomial coefficient increases row-wise and the “denominator” increases within a row.

Therefore, the coefficients of the in any given row are obtained by summing the “neighbouring” coefficients of the previous row, just as one would do it in Pascal’s triangle. The first few Bernoulli numbers are thus:

n	0	1	2	3	4	5	6
	1			0		0

Note that the Bernoulli numbers with odd index greater than 1 are all zero, not just those first two of this kind that can be seen in this table. Wikipedia has a section on why this is the case.

3. The code

For my personal usage I have compiled a small Python function that returns the k-th Bernoulli number. I have not optimized this in any way, for example, I have not added an early exit condition for odd k. I have also not added any caching here, so all numbers will always be calculated from the tip of the triangle.

def my_bernoulli(k):
    coefficients = [1]
    B = [1]
    for i in range(1, k+1):
        new_coefficients = coefficients
        coefficients = coefficients + [1]
        new_coefficients = new_coefficients + [1]
        for j in range(1, i+1):
            new_coefficients[j] = coefficients[j-1] + coefficients[j]
        sum = -B[0]
        for j in range(1, i):
            sum -= new_coefficients[j] * B[j]
        B.append(sum/new_coefficients[i])
        coefficients = new_coefficients

    return B[k]