-- Fast computation of Bernoulli numbers B_n -- http://mathworld.wolfram.com/BernoulliNumber.html -- -- http://www.haskell.org/pipermail/haskell-cafe/2003-March/004063.html -- http://www.haskell.org/pipermail/haskell-cafe/2003-March/004065.html -- http://www.haskell.org/pipermail/haskell-cafe/2003-March/004075.html module Bernoulli where import Ratio -- The 2D table of pre-computed powers = [[r^n | r<-[2..]] | n<-[1..]] -- Thanks to lazy evaluation, the table is automatically `sized' as -- needed. There is no need in any configuration parameters powers = [2..] : map (zipWith (*) (head powers)) powers -- the table of alternating sign powers = [[(-1)^r * r^n | r<-[2..]] | n<-1..] -- Again, the table is sized as needed, with no configuration neg_powers = map (zipWith (\n x -> if n then x else -x) (iterate not True)) powers -- Pre-computed Pascal triangle (again, automatically sized) pascal:: [[Integer]] pascal = [1,2,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal -- A slow algorithm: an optimized version of the code by Damien R. Sullivan -- It is based on the equality: Sum[ (n; k) B_k | k=0..n-1] = 0 b' 0 = fromIntegral 1 b' 2 = 1%6 b' n = -(sumbn n)/(fromIntegral (n+1)) sumbn:: Int -> Rational sumbn n = 1 - (fromIntegral (n+1)%(fromIntegral 2)) + sum [ (b' i) * fromIntegral(comb (n+1) i) | i <- [2,4 .. n-1] ] where comb n i = pascal!!(n-2)!!i -- Binomial coefficient -- B_(22) = (854513)/(138) -- The fast algorithm: it is faster than the above by at least 3 orders -- of magnitude. It is based on the classical double-sum formula -- over binomial coefficients -- see the Mathworld page cited above. -- The internal sum below computes with integers rather than rationals bernoulli 0 = 1 bernoulli 1 = -(1%2) bernoulli n | odd n = 0 bernoulli n = (-1)%2 + sum [ fromIntegral ((sum \$ zipWith (*) powers (tail \$ tail combs)) - fromIntegral k) % fromIntegral (k+1) | (k,combs)<- zip [2..n] pascal] where powers = (neg_powers!!(n-1)) -- To check that the denominator of the computed Bernoulli number is correct, -- we use von Staudt-Clausen theorem -- denom(B_2k) = -- product [ p | p <- takeWhile (< 2k+1) primes, 2k `rem` (p-1) == 0] primes = 2:map head (iterate sieve [3,5..]) sieve (p:xs) = [ x | x<-xs, x `rem` p /= 0 ] b_denom twok = product [ p | p <- takeWhile (<= twok1) primes, twok `rem` (p-1) == 0] where twok1 = twok + 1 -- denom(B_n)==n for 1806, but for no other n (Kellner 2005) {- "Bernoulli of 82 is " 1677014149185145836823154509786269900207736027570253414881613 % 498 "Bernoulli of 2000 is " (-6773 ...4145 more digits.. 117) % 2338224387510 "Bernoulli of 3000 is " (-2891939216292500 ... 6740 more digits.. 179981) % 12072109463901626300591430 -}