| 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).
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 10000 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 wtwjek@winternet.com
On Thu, Mar 06, 2003 at 09:23:44PM -0600, Bill Wood wrote:
I couldn't get getCPUTime (from module CPUTime) to work for me; if
Yeah, I had the same problem -- it would just return ten million or some other number, consistently. Use of getClockTime and diffClockTimes didn't help either.
anyone can enlighten me on how to get timing of function execution
For my purposes I run my program prefixed with the Unix command 'time', and with minimal output because printing 25k of rational digits seems to take a while. Having found parameters which finish in just under a minute I then rerun with output. Then I try to convince my TA this is valid. :)
BTW profiling didn't work; when I tried to compile with profiling
I was having problems with ghc 5.04; I got told to try 5.04.2 and haven't had problems since. Blessed be my sysadmins who responded quickly. To Oleg: I know about odd Bernoullis being zero and use that in the Bernoulli function itself; I hadn't thought about trying to simplify sumbn itself. Interesting... -xx- Damien X-)
. . .
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 wtwjek@winternet.com
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 Incidentally, we can show that the denominator is correct, by von Staudt-Clausen theorem:
primes = 2:map head (iterate sieve [3,5..]) sieve (p:xs) = [ x | x<-xs, x `rem` p /= 0 ]
b_denom twok = product [ p | p <- takeWhile (<= twok1) primes, twok `rem` (p-1) == 0] where twok1 = twok + 1
Here's the code (which was compiled with "ghc -O2") import Ratio import System.Environment -- powers = [[r^n | r<-[2..]] | n<-1..] powers = [2..] : map (zipWith (*) (head powers)) powers -- powers = [[(-1)^r * r^n | r<-[2..]] | n<-1..] neg_powers = map (zipWith (\n x -> if n then x else -x) (iterate not True)) powers pascal:: [[Integer]] pascal = [1,2,1] : map (\line -> zipWith (+) (line++[0]) (0:line)) pascal bernoulli 0 = 1 bernoulli 1 = -(1%2) bernoulli n | odd n = 0 bernoulli n = (-1)%2 + sum [ fromIntegral ((sum $ zipWith (*) powers (tail $ tail combs)) - fromIntegral k) % fromIntegral (k+1) | (k,combs)<- zip [2..n] pascal] where powers = (neg_powers!!(n-1)) main = do [arg] <- getArgs let n = (read arg)::Int print $ "Bernoulli of " ++ (show n) ++ " is " print (bernoulli n)
. . .
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 wtwjek@winternet.com
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
participants (4)
-
Bill Wood -
Damien R. Sullivan -
oleg＠pobox.com -
Tom Moertel