speedup help

[oleg@pobox.com]

> 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)