{-# LANGUAGE ParallelListComp #-}

import Data.List 
import Data.Ratio

choose :: Integer -> Integer -> Integer
n `choose` k = product [n-k+1..n] `div` product [2..k]

bernoullis = 1 : [sum [(n `choose` k) % (k - n - 1) * b | k <- [0..] | b <- bs]
                    | bs <- tail (inits bernoullis) | n <- [1..]]

bernoulli = genericIndex bernoullis

{-
Alternative solution using array memoisation:

import Data.Array

bernoulli n = bernoulli' n
  where a = listArray (0, n) (map bernoulli' [0..n])
        bernoulli' 0 = 1
        bernoulli' n = sum [(n `choose` k) % (k - n - 1) * a ! k | k <- [0..n-1]]

-}