speedup help

oleg@pobox.com oleg@pobox.com
Fri, 7 Mar 2003 17:58:58 -0800 (PST)

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