# speedup help

**[email protected]
**
[email protected]

*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
let n = (read arg)::Int
print $ "Bernoulli of " ++ (show n) ++ " is "
print (bernoulli n)