Oleg's blew the heap at 847; mine valiantly struggled on 'til it blew the heap at 1910.
Hmm, I managed to compute bernoulli 2000 and even bernoulli 3000. The code is included. It took 7 minutes (2GHZ Pentium IV, 1GB memory) to compute bernoulli 2000 and 33 minutes for bernoulli 3000. I monitored the memory usage of the compiled application using top and found that the resident set stayed at 30MB, which is a little bit less than the resident set for Emacs. My emacs has a dozen of open windows, and has been running for a month. Just for the record, here's the result of bernoulli 3000: (-2891939 ...6744 other digits... 81) % 12072109463901626300591430 Incidentally, we can show that the denominator is correct, by von Staudt-Clausen theorem:
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
Here's the code (which was compiled with "ghc -O2") import Ratio import System.Environment -- powers = [[r^n | r<-[2..]] | n<-1..] powers = [2..] : map (zipWith (*) (head powers)) powers -- powers = [[(-1)^r * r^n | r<-[2..]] | n<-1..] neg_powers = map (zipWith (\n x -> if n then x else -x) (iterate not True)) powers pascal:: [[Integer]] pascal = [1,2,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal 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)) main = do [arg] <- getArgs let n = (read arg)::Int print $ "Bernoulli of " ++ (show n) ++ " is " print (bernoulli n)