# 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: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).

```