Re: speedup help
Bill Wood wrote: I think I got the right results for B_3000: [...] Mathematica 4.1 computes B_3000 as follows: In[1]:= BernoulliB[3000] Out[1]= -28919392162925009628147618267854828678617917853903846822112332719169192942048\ 518130533026045150594816676476469224548430690874860242714680177615276168526041\ [ 83 lines omitted ] 535500476970565917875995082926214969042647564930033701717898024811160065586065\ 5536080415376091806331620179538459843417141322454179981 / 12072109463901626300591430 Cheers, Tom ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
. . . Oleg's blew the heap at 847; mine valiantly struggled on 'til it blew the heap at 1910. Hmm, I managed to compute bernoulli 2000 and even bernoulli 3000. The code is included. It took 7 minutes (2GHZ Pentium IV, 1GB memory) to compute bernoulli 2000 and 33 minutes for bernoulli 3000. I monitored the memory usage of the compiled application using top and found that the resident set stayed at 30MB, which is a little bit less than the resident set for Emacs. My emacs has a dozen of open windows, and has been running for a month. Just for the record, here's the result of bernoulli 3000: (-2891939 ...6744 other digits... 81) % 12072109463901626300591430 Well, kudos to Oleg! I just ran his code from the msg this one replies to and got similar results: On a 650 Mhz Athlon with 128MB RAM I compiled with ghc -O (GHC 5.00.1). Using the time command, bernoulli 2000 took 490 sec. user time while bernoulli 3000 took 2170 sec. user time. Notice there were no heap overflows this time. I hope someone can explain the differences in space behavior between the version in Oleg's mail of March 6 and this version. I'm surprised we are as close as we are in time given that my processor is less than half as fast. I would also expect that my having 1/8-th the memory he has would slow me down more due to page faulting. The current version also fixes a minor glitch: the earlier version gave B_0 = 0 instead of 1. I think I got the right results for B_3000: My print-out had the same denominator along with a 6762-digit numerator with the same initial seven and final two digits. I don't get 6744 digits in the middle, however. I'm impressed by the good performance characteristics of high-level Haskell code. Nice work Oleg, -- Bill Wood [EMAIL PROTECTED] ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
| 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). ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
Oleg has a very interesting approach; in particular, he avoids
explicit recursion and (most) computing with rationals, while
also replacing the factorials in binary coefficients by using
successive rows of Pascal's triangle. He also skips over the
B_{2k+1}, k 0, which are all 0.
I slogged through the standard expansions, deriving a tail
recursive function that rolls along successive Bernoulli numbers,
generating successive rows of Pascal's triangle along the way,
and returning the list of B_n .. B_0. You can think of the list
of Bernoulli numbers as memoizing just-in-time.
Running it in Hugs on a 650Mhz Athlon with 128M RAM, bernoulli 82
took ca. 1 sec. Compiling with ghc -O, bernoulli 1000 took approx.
15 sec. wall time; bernoulli 1 blew the heap.
I couldn't get getCPUTime (from module CPUTime) to work for me; if
anyone can enlighten me on how to get timing of function execution
I'd appreciate it.
BTW profiling didn't work; when I tried to compile with profiling
flags I received error msgs saying that interface files for standard
libraries couldn't be found. Compiling without the flags seems to
work just fine.
Oleg's program brings up another interesting point for all you
mathematicians out there: I found two basically different expansion
formulas for Bernoulli numbers. One is based on the formal
expansion of the equation
(B + 1)^n = B^n
which results in binomial-theorem expansions all of whose terms are
positive. The other is based on formal expansion of the equation
(B - 1)^n = B^n
which results in binomial-theorem expansions whose terms alternate
in sign. The amazing thing is, they two sets of numbers only differ
at one term: the first formula results in B_1 = -1/2 while the
second results in B_1 = +1/2 !
I found the second formula in Conway Guy, _The Book of Numbers_,
p.108.
The second formula, with tiny variations, can be found in Knuth,
_Fundamental Algorithms_, p. 109, Abramowitz Stegun, _Handbook
of Mathematical Functions_, p. 804 and Song Y. Yan, _Number Theory for
Computing_, p. 75
This has gotten a little long, sorry. If you want I can post my Haskell
code or send privately.
-- Bill Wood
[EMAIL PROTECTED]
___
Haskell-Cafe mailing list
[EMAIL PROTECTED]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
. . . This code seems to compute 'bernoulli 82' in less then a second, in Hugs (on a 2GHz Pentium IV). Just a note: I compiled and ran Oleg's and mine for comparison. The two seem to be of the same complexity, with Oleg's a little faster (modulo my using wall time; see previous msg.) Oleg's blew the heap at 847; mine valiantly struggled on 'til it blew the heap at 1910. I must be doing something right, since I'm carrying around the list of all numbers from B_0 through B_n, while Oleg cleverly avoids that. I was also surprised to see Oleg's blow the heap at an *odd* value -- I thought he skipped those. -- Bill Wood [EMAIL PROTECTED] ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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-) ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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 | [EMAIL PROTECTED]
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
[EMAIL PROTECTED]
http://www.haskell.org/mailman/listinfo/haskell-cafe
___
Haskell-Cafe mailing list
[EMAIL PROTECTED]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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 -- [EMAIL PROTECTED] -- http://www.cit.gu.edu.au/~arock/ School of Computing and Information Technology Griffith University -- Nathan, Brisbane, Queensland 4111, Australia ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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 ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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 / [EMAIL PROTECTED] / 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 scarcityI think we should embrace the era of plenty, and work out how to mutually live in it. -- john gilmore ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help
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-) ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: speedup help update
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.50.0 comb Main 21.89.7 sumbn Main 10.7 16.4 Time to look at arrays, finally. -xx- Damien X-) ___ Haskell-Cafe mailing list [EMAIL PROTECTED] http://www.haskell.org/mailman/listinfo/haskell-cafe
