# Bernoulli numbers (Haskell)

### From LiteratePrograms

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.

This version is probably better though.

## Memoize

We use a lazy list to memoize calculations. This takes linear time to lookup, instead of constant time, but oh well; it's not a big deal for small numbers. We use "genericIndex" instead of the (!!) operator because (!!) only takes Ints, which will force the rest of our program to use Ints, which will overflow when numbers get sufficiently big, producing incorrect answers; whereas if arbitrary-precision Integers are used, it will not overflow.

<<memoize>>=bernoulli :: Integer -> Ratio Integer bernoulli = genericIndex bernoulli_table bernoulli_table = map bernoulli'[0..]

## Implementation

Now to implement the actual function. The computation can be performed more efficiently according to the following identities due to Srinivasa Ramanujan:

where

<<impl>>=bernoulli' :: Integer -> Ratio Integer bernoulli' 0 = 1 bernoulli' 1 = -1 % 2 bernoulli' n | odd n = 0 | n `mod` 6 == 0 =letval =(n+3)% 3 - bernoulliA n(n `div` 6)inval / fromIntegral(choose(n+3)n)| n `mod` 6 == 2 =letval =(n+3)% 3 - bernoulliA n((n-2)`div` 6)inval / fromIntegral(choose(n+3)n)| n `mod` 6 == 4 =letval = -(n+3)% 6 - bernoulliA n((n-4)`div` 6)inval / fromIntegral(choose(n+3)n)

Now we implement A:

<<a>>=--This can be sped up by not re-calculating the binary coefficient every time or --memoizing that function. bernoulliA :: Integer -> Integer -> Ratio Integer bernoulliA n m = sum(map foo[1..m])wherefoo j = fromIntegral(choose(n+3)(n-6*j))* bernoulli(n-6*j)

The last piece we need is to calculate the binary coefficient:

<<bincoef>>=choose ::(Integral a)=> a -> a -> a choose n r = product[n-r'+1..n]`div` product[1..r']wherer' = min r(n-r)

Now we put it all together; remembering to import Data.Ratio to use the "%" ratio operator, and Data.List for genericIndex:

<<bernoulli.hs>>=importData.RatioimportData.List memoize impl a bincoef

## Running

Now the code can be run from the Haskell shell:

Prelude> :l bernoulli.hs [1 of 1] Compiling Main ( bernoulli.hs, interpreted ) Ok, modules loaded: Main. *Main> bernoulli 6 1%42 *Main> bernoulli 12 (-691)%2730 *Main> bernoulli 16 (-3617)%510 *Main> bernoulli 18 43867%798