I am a big fan of the Akiyama–Tanigawa algorithm¹² for computing Bernoulli numbers. It works by generating a structure resembling Pascal’s triangle, with rows whose entries are determined by a simple calculation on the two entries above them, from which we can read off the sequence of Bernoulli numbers.

What I like about this algorithm, besides its simplicity and elegance, is that it admits a concordantly simple and elgant implementation in Haskell using infinite data structures.

Let’s describe the structure the algorithm emerges from. Like Pascal’s triangle it has infinitely many rows. Unlike Pascal’s triangle, its rows are also infinite. Specifically, it is a rightwards- and downwards-infinite grid of rational numbers generated as follows:

1. The first row consists of the reciprocals of the positive integers.
2. In any subsequent row, the n-th element is given by the difference of the previous row’s n-th and (n+1)-st elements, all multiplied by n.

The first entries of each row—the leftmost vertical column—are the Bernoulli numbers.

If you want to generate the first n Bernoulli numbers, then you can start with just n entries in the initial row. This is illustrated below for n=5; our first five Bernoulli numbers are 1, 1/2, 1/6, 0, -1/30, as we expect.

  1    1/2   1/3   1/4   1/5
  ↓  ↙  ↓  ↙  ↓  ↙  ↓  ↙  
 1/2   1/3   1/4   1/5
  ↓  ↙  ↓  ↙  ↓  ↙  
 1/6   1/6   3/20  
  ↓  ↙  ↓  ↙ 
  0    1/30  
  ↓  ↙    
-1/30   

-- for example, the third element of the third row is
-- 3*(1/4 - 1/5) = 3/20

However, Haskell’s infinite data structures and lazy evaluation lets us define the whole bi-infinte array—hence all the Bernoulli numbers—in one go.³ Let’s see how this is done.⁴

First we import Data.Ratio to handle rationals (representing p/q with p % q) and initialize an infinite array of the reciprocals of natural numbers:

import Data.Ratio

recips :: [Rational]
recips = map (\n -> 1 % n) [1..]

Next we write our row update rule and iterate it infinitely on the initial row to get the bi-infinite grid:

nextRow :: [Rational] -> [Rational]
nextRow xs = zipWith (*) [1..] (zipWith (-) xs (drop 1 xs))

bernoulliGrid :: [[Rational]]
bernoulliGrid = iterate nextRow recips

Taking the first entries of each row in the grid gives the series of Bernoulli numbers:

bernoulliNums :: [Rational]
bernoulliNums = map head bernoulliGrid

Now we can get as many Bernoulli numbers as we want with take <num> bernoulliNums. And indeed we have:

Prelude Data.Ratio> take 5 bernoulliNums 
[1 % 1,1 % 2,1 % 6,0 % 1,(-1) % 30]

Neat!