# speedup help

**Mark P Jones
**
mpj@cse.ogi.edu

*Mon, 3 Mar 2003 20:46:22 -0800*

Hi Damien,
|* So, I'm having to calculate 'n choose k' an awful lot. At
*|* the moment I've got this:
*|*
*|* comb :: Integer -> Integer -> Integer
*|* comb m 0 = 1
*|* comb m n = (numerator(toRational (fact m) / toRational (fact
*|* n * fact (m-n))))
*|*
*|* where fact is a memoized factorial function. It's not
*|* perfectly memoized, though; I use lists, since that's easier
*|* by default. They should be arrays, and possibly just
*|* changing that would speed comb up a lot. (Comb is currently
*|* 40% of runtime, fact is 23%.) But I think it should be
*|* possible to speed up comb itself, too.
*
If you're looking to calculate memoized n choose k, then I'd suggest
looking at the following. First, where do the n choose k numbers come
from? Forget factorial! Think Pascal's triangle!
pascal :: [[Integer]]
pascal = iterate (\row -> zipWith (+) ([0] ++ row) (row ++ [0])) [1]
Now we can define a variant of your comb function like this:
comb :: Int -> Int -> Integer
comb n m = pascal !! n !! m
(Add an extra line if you want comb to work for values outside the
range: 0 <= m <= n. You could also replace the rows of the triangle
with arrays, if you want. No factorials, multiplies, or divides here,
and natural memoization ...)
|* comb is only called from here:
*|* sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i
*|* <- [0 .. n-1] ]
*
In that case, you can take further advantage of using Pascal's triangle
by recognizing that numbers of the form (comb (n+1) i) are just the
entries in the (n+1)th row. (All but the last two, for reasons I
don't understand ... did you possibly want [0..n+1]?) So we get the
following definition:
sumbn n = sum [ bernoulli i * fromIntegral c
| (i,c) <- zip [0..n-1] (pascal!!(n+1)) ]
Actually, I prefer the following version that introduces an explicit
dot product operator:
sumbn n = take n (map bernoulli [0..]) `dot` (pascal!!(n+1))
dot xs ys = sum (zipWith (*) xs ys)
I don't have your setup to test the performance of these definitions,
but I'd be curious to know how they fare. Even if they turn out to be
slower, I thought these definitions were interesting and different
enough to justify sharing ...
All the best,
Mark