# speedup help

**oleg@pobox.com
**
oleg@pobox.com

*Thu, 6 Mar 2003 16:32:44 -0800 (PST)*

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