speedup help

[oleg@pobox.com]

| comb is only called from here:
| sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i 
| <- [0 .. n-1] ]

Probably I misunderstand what "bernoulli i" stands for. If it is meant
Bernoulli number B_i,
	http://mathworld.wolfram.com/BernoulliNumber.html
then the above expression is quite inefficient. Bernoulli numbers with
odd indices are all zero, except B_1, which is -1/2. Therefore, the above
expression ought to be written as

sumbn n = 1 - (fromIntegral (n+1)%(fromIntegral 2)) + 
    sum [ (b' i) * fromIntegral(comb (n+1) i) | i  <- [2,4 .. n-1] ]

It appears that you compute the sumbn to obtain B_n from the equality

	sum [bernoulli i * (comb (n+1) i) | i<-[0..n]] == 0

Have you tried 
	bernoulli n = sum [ (sumib k n) % (k+1) | k<-[1..n]]
where
       sumib k n = sum [ (comb k r) * r^n | r<-[2..k]] 
	         - sum [ (comb k r) * r^n | r<-[1,3..k]]
the advantage of the latter series is that sumib is an Integer, rather
than a rational. The powers of r can be memoized.


Here's the code

-- powers = [[r^n | r<-[1..]] | n<-[1..]]
powers = [1..] : map (zipWith (*) (head powers)) powers

-- neg_powers = [[(-1)^r * r^n | r<-[1..]] | n<-[1..]]
neg_powers = 
  map (zipWith (\n x -> if n then x else -x) (iterate not False)) powers

pascal:: [[Integer]]
pascal = [1,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal

bernoulli n = sum [ fromIntegral (sumib k n) % fromIntegral (k+1) | k<-[1..n]]

sumib:: Int -> Int -> Integer
sumib k n = sum $ zipWith (*) (neg_powers!!(n-1)) (tail $ pascal!!(k-1))

This code seems to compute 'bernoulli 82' in less then a second, in
Hugs (on a 2GHz Pentium IV).