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. comb is only called from here: sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i <- [0 .. n-1] ] Here was one try: fcomb :: Integer -> Integer -> Integer fcomb m 0 = 1 fcomb m n = res where res = last * (m-n+1) `div` n last = res except, obviously, this doesn't work. I hope it's clear what I'm trying to do, or what I would be in a more imperative language -- in C I'd probably have some static variable in fcomb. I figure monads are needed, but I've been unable to figure them out enough to apply them here. Will the monadism propagate all the way up and require changing lots of function types? Bleah. I'm using ghc, can I sneak some mutable in here instead? Any advice? Also on using arrays, where my parameters come off the command line. I imagine in C++ I'd just precompute a bunch of tables and then just use those values in the actual algorithm. Thanks! -xx- Damien X-) (if you're curious, this is for a class, but not a class on using Haskell. I chose to use Haskell for this assignment after ghc -O, to my surprise, outperformed ocaml. I suspect GMP deserves the real credit, but hey.)
I think you would get a big speed-up if you got rid of all the rational stuff and just used: comb m n = fact m `div` (fact n * fact (m-n)) If that doesn't speed it up enouch, you can of course cache fact m in your computation and do something like: sumbn n = sum [ bournoulli i * fm `div` (fn * fact (m-n)) | i <- [0..n-1]] where fm = fact m fn = fact n it is possible that the compiler is inlining the call the comb, in which case this has already been done for you. hard to say for sure. putting '{-# INLINE comb #-}' might help a lot.
From there, you should probably look at arrays if you can bound n.
-- Hal Daume III | hdaume@isi.edu "Arrest this man, he talks in maths." | www.isi.edu/~hdaume On Mon, 3 Mar 2003, Damien R. Sullivan wrote:
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.
comb is only called from here: sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i <- [0 .. n-1] ]
Here was one try:
fcomb :: Integer -> Integer -> Integer fcomb m 0 = 1 fcomb m n = res where res = last * (m-n+1) `div` n last = res
except, obviously, this doesn't work. I hope it's clear what I'm trying to do, or what I would be in a more imperative language -- in C I'd probably have some static variable in fcomb. I figure monads are needed, but I've been unable to figure them out enough to apply them here. Will the monadism propagate all the way up and require changing lots of function types? Bleah. I'm using ghc, can I sneak some mutable in here instead?
Any advice? Also on using arrays, where my parameters come off the command line. I imagine in C++ I'd just precompute a bunch of tables and then just use those values in the actual algorithm.
Thanks!
-xx- Damien X-)
(if you're curious, this is for a class, but not a class on using Haskell. I chose to use Haskell for this assignment after ghc -O, to my surprise, outperformed ocaml. I suspect GMP deserves the real credit, but hey.) _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
G'day all. On Mon, Mar 03, 2003 at 04:59:21PM -0800, Hal Daume III wrote:
I think you would get a big speed-up if you got rid of all the rational stuff and just used:
comb m n = fact m `div` (fact n * fact (m-n))
Or, even better, if you didn't multiply stuff that you're just going to divide out in the first place. choose :: (Integral a) => a -> a -> Integer choose m n | m < 0 = 0 | n < 0 || n > m = 0 | n > m `div` 2 = choose' n (m-n) | otherwise = choose' (m-n) n where choose' i' j' = let i = toInteger i' j = toInteger j' in productRange (i+1) j `div` factorial j factorial :: (Integral a) => a -> Integer factorial n = productRange 1 n productRange :: (Integral a) => Integer -> a -> Integer productRange b n | n < 5 = case n of 0 -> 1 1 -> b 2 -> b*(b+1) 3 -> (b*(b+1))*(b+2) 4 -> (b*(b+3))*((b+1)*(b+2)) _ -> 0 | otherwise = let n2 = n `div` 2 in productRange b n2 * productRange (b+toInteger n2) (n-n2) Note that productRange uses a divide-and-conquer algorithm. The reason for this is that it's more efficient to multiply similarly-sized Integers than dissimilarly-sized Integers using GMP. Cheers, Andrew Bromage
On Tue, Mar 04, 2003 at 12:25:01PM +1100, Andrew J Bromage wrote:
Or, even better, if you didn't multiply stuff that you're just going to divide out in the first place.
I had thought of that before, and used a simple comb m n = product [m, m-1 .. m-n+1] / fact (m-n) but the unmemoized product proved to be slower than the original. This looks rather different, though. :) -xx- Damien X-)
On Mon, Mar 03, 2003 at 04:59:21PM -0800, Hal Daume III wrote:
comb m n = fact m `div` (fact n * fact (m-n))
This was the biggest help, 33 seconds instead of my original 43. fact is the big consumer now, and I think cries out for being arrayed, especially as it gets used a lot elsewhere too.
If that doesn't speed it up enouch, you can of course cache fact m in your computation and do something like:
sumbn n = sum [ bournoulli i * fm `div` (fn * fact (m-n)) | i <- [0..n-1]] where fm = fact m fn = fact n
I'm not sure what this is doing. i has to be in the comb part.
From there, you should probably look at arrays if you can bound n.
Bound at compile time or at some early point in run time? The program's behavior is determined by command line arguments, and filling an array with n factorials would be perfectly appropriate. I'm sorry to report that the other suggestions didn't help much. Andrew Rock's took 80 seconds. Andrew Bromage's did gain a slight improvement -- 41 seconds instead of 43. If I replace the factorial in in productRange (i+1) j `div` factorial j with my own fact then it goes to 37 seconds. But that's still more time than Hal's simple use of Integers. Top 3 functions of the version with Hal's code are: fact Main 28.5 0.0 comb Main 21.8 9.7 sumbn Main 10.7 16.4 Time to look at arrays, finally. -xx- Damien X-)
I would like to say that programming this project (calculate the Euler-Mascheroni constant to as many digits as possible in a minute) in Haskell has been fairly pleasant overall. The startup time was a bit oogy -- ocaml was faster to get a working program -- and ghc's compilation time continues to aggravate me (10x slower than ocaml, or even slower with -O -- but then the code's faster), and I fear getting my code really efficient will be a pain. But the code itself is nice and elegant and it was neat being able to do sums on list comprehensions. My code tends to map pretty straightforwardly to the math on paper, which is really nice for being assured of correctness. -xx- Damien X-)
G'day all. On Mon, Mar 03, 2003 at 09:53:24PM -0500, Damien R. Sullivan wrote:
Time to look at arrays, finally.
This might help: http://haskell.org/wiki/wiki?MemoisingCafs Cheers, Andrew Bromage
On Tuesday, March 4, 2003, at 10:26 AM, Damien R. Sullivan wrote:
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.
comb is only called from here: sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i <- [0 .. n-1] ]
Here was one try:
fcomb :: Integer -> Integer -> Integer fcomb m 0 = 1 fcomb m n = res where res = last * (m-n+1) `div` n last = res
Try this: comb :: Integral a => a -> a -> a comb n r = c n 1 1 where c n' r' p | r' > r = p | otherwise = c (n' - 1) (r' + 1) (p * n' `div` r') Cheers, Rock. -- Andrew Rock -- A.Rock@cit.gu.edu.au -- http://www.cit.gu.edu.au/~arock/ School of Computing and Information Technology Griffith University -- Nathan, Brisbane, Queensland 4111, Australia
I have no idea if the following is faster or not (I suspect not), but it is certainly easier to read: n `choose` k = (n `permute` k) `div` (fact k) n `permute` k = product [(n-k+1) .. n] fact n = product [1 .. n] mike -- mike castleman / mlc67@columbia.edu / http://mlcastle.net aolim: mlcastle / icq: 3520821 / yahoo: mlc000 "we have invented the technology to eliminate scarcity, but we are deliberately throwing it away to be benefit those who profit from scarcity....I think we should embrace the era of plenty, and work out how to mutually live in it." -- john gilmore
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
On Mon, Mar 03, 2003 at 08:46:22PM -0800, Mark P Jones wrote:
pascal :: [[Integer]] pascal = iterate (\row -> zipWith (+) ([0] ++ row) (row ++ [0])) [1]
comb :: Int -> Int -> Integer comb n m = pascal !! n !! m
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
No, the sum for a Bernoulli number is a combination times another Bernoulli number from 0 to n-1. Hard to have B_n depend on B_n. At least in a nice recurrence...
sumbn n = sum [ bernoulli i * fromIntegral c | (i,c) <- zip [0..n-1] (pascal!!(n+1)) ]
This code as is takes about 23 seconds, comparable to the 22 seconds of factorial with array (hardcoded, since I can't get it dynamically in a pretty fashion.) If I turned pascal into arrays it might be even faster. I'd have to change something though, right, zipWith wouldn't work with arrays?
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)
This needed some modification, since bernoulli returns Rationals, so I had zipWith use a special mult function. It just took 25 seconds.
slower, I thought these definitions were interesting and different enough to justify sharing ...
Hey, you're even faster too! At least for messing with comb. Aaron Denney, to his credit, had a pretty similar idea a week ago, but I didn't get what he was talking about then. Newbies like code they can paste in. :) Thanks! -xx- Damien X-)
participants (6)
-
Andrew J Bromage -
Andrew Rock -
Damien R. Sullivan -
Hal Daume III -
Mark P Jones -
mike castleman