Bernoulli numbers (Haskell)

From LiteratePrograms

Jump to: navigation, search
Other implementations: Erlang | Haskell | Python

The Bernoulli numbers Bm are a sequence of rational numbers defined by the identity

with the initial condition B0 = 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 Bm 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 = let val = (n+3) % 3 - bernoulliA n (n `div` 6) in
                                  val / fromIntegral (choose (n+3) n)
             | n `mod` 6 == 2 = let val = (n+3) % 3 - bernoulliA n ((n-2) `div` 6) in
                                  val / fromIntegral (choose (n+3) n)
             | n `mod` 6 == 4 = let val = -(n+3) % 6 - bernoulliA n ((n-4) `div` 6) in
                                  val / 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])
  where foo 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']
  where r' = 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>>=
import Data.Ratio
import Data.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
Retrieved from "http://literateprograms.org/Bernoulli_numbers_%28Haskell%29"
Views