[sage-devel] Re: computing the number of partitions of an integer
Since you've extended this thread some more, I'll mention that the current version of the code that I am working on runs a little bit more than twice as fast as that version. It's posted at http://www.math.lsa.umich.edu/~bober/partitions_c.cc if anyone is interested in testing it. (But don't waste too much time testing right now, because it is likely to change.) It's not quite ready for inclusion, but, eventually, maybe Bill Hart's prediction about it only taking 10 seconds to compute p(10^9) on sage.math will come true. On Fri, 2007-08-03 at 23:21 -0500, Alec Mihailovs wrote: From: William Stein [EMAIL PROTECTED] Thanks. I've fixed this for sage = 2.7.3 in the attached patch. Just upgraded SAGE to 2.7.3 and number_of_partitions works twice as fast as in 2.7.2 with algorithm='bober' for 10^9 - it took 112 sec in 2.7.2 and only 56 sec in 2.7.3. This is in Ubuntu 7.04 on Dual Core 64-bit 2.6GHz Athlon with 2 GB RAM. By the way, to compare with other systems, SAGE 2.7.2 installation took about 53 minutes. I noticed the following interesting thing with time measuring: assigning number_of_partitions to a variable (that shouldn't take time at all), took about 10% or more of total time. For example, sage: time number_of_partitions(10^7) took about 0.8 sec while sage: time a=number_of_partitions(10^7) took slightly more than 1 sec. For 10^8 the corresponding times were 5 sec and 5.45 sec. That seems odd, isn't it? Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From: Jonathan Bober [EMAIL PROTECTED] It's not quite ready for inclusion, but, eventually, maybe Bill Hart's prediction about it only taking 10 seconds to compute p(10^9) on sage.math will come true. That's fantastic! Great work! Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From: William Stein [EMAIL PROTECTED] Thanks. I've fixed this for sage = 2.7.3 in the attached patch. Just upgraded SAGE to 2.7.3 and number_of_partitions works twice as fast as in 2.7.2 with algorithm='bober' for 10^9 - it took 112 sec in 2.7.2 and only 56 sec in 2.7.3. This is in Ubuntu 7.04 on Dual Core 64-bit 2.6GHz Athlon with 2 GB RAM. By the way, to compare with other systems, SAGE 2.7.2 installation took about 53 minutes. I noticed the following interesting thing with time measuring: assigning number_of_partitions to a variable (that shouldn't take time at all), took about 10% or more of total time. For example, sage: time number_of_partitions(10^7) took about 0.8 sec while sage: time a=number_of_partitions(10^7) took slightly more than 1 sec. For 10^8 the corresponding times were 5 sec and 5.45 sec. That seems odd, isn't it? Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From: Jonathan Bober [EMAIL PROTECTED] Incidentally, I should have mentioned here that I submitted a patch for version .4, and also updated it at http://www.math.lsa.umich.edu/~bober/partitions_c.cc It uses long doubles now when then precision is small enough (and then, later, just doubles like before), and the speedup is significant. By the way, I just installed SAGE 2.7.2 (from source) on my Ubuntu 7.04 system, with AMD Athlon 64 X2 5200+ and tested the version included there (using algorithm='bober'). It has the following feature - for input 2^32 it gave the ValueError saying that the input should be less than 2^64-2 (evaluated). Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 8/2/07, Alec Mihailovs [EMAIL PROTECTED] wrote: It uses long doubles now when then precision is small enough (and then, later, just doubles like before), and the speedup is significant. By the way, I just installed SAGE 2.7.2 (from source) on my Ubuntu 7.04 system, with AMD Athlon 64 X2 5200+ and tested the version included there (using algorithm='bober'). It has the following feature - for input 2^32 it gave the ValueError saying that the input should be less than 2^64-2 (evaluated). Thanks. I've fixed this for sage = 2.7.3 in the attached patch. William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~--- 5561.patch Description: Binary data
[sage-devel] Re: computing the number of partitions of an integer
On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote: [...] Checking partition count computation: On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9: Mathematica 5.2 (PartitionsP[10^9]:95.5115 s Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s Jon Bober's code (i386: 'jb 10'): 160 s I wonder why our Core 2 Duo timings are so different? All of the executables are i386 binaries (32-bit x86). One oddity: Jon Bober's code produced a value ending in '30457526857797923685688339', while Mathematica and SAGE produced one ending in '30457526831036667979062760', so they differ. Addendum to last email: It looks like the 2.7.2.1 is using the new code by default, so the 'sage' time above should be the time for the newest version of the code, so I'm going to assume that the other timing is for an old version of my code. (I running with algorithm='pari' right now and it is taking a really long time.) While I'm assuming things, I'm also going to assume that you must have mixed up the output values, and thus Mathematica must be producing the wrong answer, in which case its no longer fair to compare running times to Mathematica. (Of course, Mathematica does have the number ending in 339 in the Mathematica book.) Justin -- Justin C. Walker, Curmudgeon at Large Institute for the Absorption of Federal Funds --- If it weren't for carbon-14, I wouldn't date at all. --- --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Jul 31, 2007, at 22:54 , Jonathan Bober wrote: On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote: On Jul 31, 2007, at 18:36 , William Stein wrote: [snip] On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9: Mathematica 5.2 (PartitionsP[10^9]:95.5115 s Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s Jon Bober's code (i386: 'jb 10'): 160 s I wonder why our Core 2 Duo timings are so different? That is puzzling. Are you sure that you have the latest version of the code? The source file should have 'Version .3' at the top. Yup. Double-checked the file; and re-ran the test. Essentially the same time. I'll check this tomorrow, along with the Mathematica oddity. You can get the latest version, or at least that latest version minus some memory leak fixes, via hg_sage.pull(), or from http://www.math.lsa.umich.edu/~bober/partitions_c.cc The last version I posted to the list was significantly slower than what I have now (especially since I think I accidently had significant improvements commented out in the last code I posted to the list.) That was the the Jon Bober timing above. I'll try the latest one later. Justin -- Justin C. Walker, Curmudgeon-At-Large Institute for the Absorption of Federal Funds Men are from Earth. Women are from Earth. Deal with it. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 8/1/07, Justin C. Walker [EMAIL PROTECTED] wrote:
The last version I posted to the list was significantly slower than
what
I have now (especially since I think I accidently had significant
improvements commented out in the last code I posted to the list.)
That was the the Jon Bober timing above. I'll try the latest one
later.
Two comments for this thread:
(1) The PARI in SAGE is the latest stable release, which means
number_of_partitions(n, algorithm='pari') is *wrong* even on
some fairly small values of n.
(2) If you do
sage: hg_sage.apply('http://sage.math.washington.edu/home/was/sage.hg')
you'll get an update to my latest code (not 100% doc tested!) which
has the latest version of jon's code, and for which
number_of_partitions(n)
uses Jon's code by default.
William
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Hi, Regarding the discrepancy in timings between me and Justin Walker using OS X *intel* Mathematica, it turns out I was running Mathematica on my laptop under OS X via Rosetta, so those times should be ignored. The timings I've posted under Linux are all fine. SO, currently SAGE and Mathematica on OSX Intel MacBook Pro 2.33Ghz take almost exactly the same amount of time (95 seconds) to compute p(10^9). On 64-bit Linux Mathematica is still faster than SAGE, so probably there is a 64-bit only precision trick that is relevant... -- William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Tue, 2007-07-31 at 14:24 -0700, Bill Hart wrote: I do highly recommend this quad double library by the way. And they've implemented all manor of transcendental functions too!! The quad- doubles would give you 206 bits, even on your machine. Bill. URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~--- Have you found the quad double library to be faster than mpfr, or just more convenient? I tried using it in the partition counting code, and it actually slowed things down when I used it for all computations between 200 and 64 bits. Alternately, if I just use it between 200 and, say, 180 bits, it gives almost no change in speed. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Incidentally, I should have mentioned here that I submitted a patch for version .4, and also updated it at http://www.math.lsa.umich.edu/~bober/partitions_c.cc It uses long doubles now when then precision is small enough (and then, later, just doubles like before), and the speedup is significant. On Wed, 2007-08-01 at 12:34 -0700, Justin C. Walker wrote: On Aug 1, 2007, at 12:29 PM, Justin C. Walker wrote: On Jul 31, 2007, at 10:54 PM, Jonathan Bober wrote: On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote: On Jul 31, 2007, at 18:36 , William Stein wrote: [snip] That is puzzling. Are you sure that you have the latest version of the code? I downloaded the .3 source and ran it on my 2.33 GHz Core 2 Duo. I compiled the program with several -O settings, and ran them with the argument '10': Adding to the background noise: On a 2.33GHz Core 2 Duo (Mac OS X, 10.4.10): Mathematica 6.0: 83s Mathematica 6.0.1: 77s Justin -- Justin C. Walker, Curmudgeon-At-Large Institute for the Enhancement of the Director's Income Experience is what you get when you don't get what you want. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Aug 1, 2007, at 1:26 PM, Jonathan Bober wrote:
On Wed, 2007-08-01 at 12:29 -0700, Justin C. Walker wrote:
[snip]
Since the time with the self-compiled code is so similar to the time
with the sage-compiled code, I would guess that you are linking to the
sage mpfr library, which wasn't compiled with -O2 (or at least linking
to another version of mpfr that wasn't compiled with optimizations.)
Oops. You're right. I forgot about this part of the discussion.
I'll get an 'mpfr' library compiled with more optimization.
Has anyone looked at the differences between the various optimization
choices? On Mac OS X, there is -Ox, for x in {0,1,2,3,s,z}, -O,
and -fast.
Also, did you double check the answers? If my code is giving the wrong
answer on your computer, there could be some processor specific bug
somewhere.
I'm verifying this now, but I think that my claim (that Mathematica
and sage agreed and computed a different value than your code) is a
typo. It should have said 'sage', which at the time was using Pari
(and according to Stein, was generating incorrect results).
Yup: all versions of Mathematica (5.2, 6.0, 6.0.1) agree with you :-}
Justin
--
Justin C. Walker, Curmudgeon-At-Large
Director
Institute for the Enhancement of the Director's Income
Weaseling out of things is what separates us from the animals.
Well, except the weasel.
- Homer J Simpson
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Certainly no, not everything from a public institution is in the public domain. This should be analyzed case by case. In case of doubt it would be better to ask the author. Pablo On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
In this case, the license does not says that it is in the public domain (but that it is a copyrighted work!), but you can use it as AS IS. I think that the only condition that is imposed to us is to include the declaimer. Pablo On 7/31/07, Pablo De Napoli [EMAIL PROTECTED] wrote: Certainly no, not everything from a public institution is in the public domain. This should be analyzed case by case. In case of doubt it would be better to ask the author. Pablo On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Oh well, I don't understand all this licensing stuff. So do you understand this license? What about derived works. Does that mean it is not possible to modify this library and redistribute the modified version? In particular, as a C++ library it is no use to me unmodified. There are also some functions I'd like to add, and the calling convention is the opposite to what I use. So basically I'd have to modify it to make it useful. Bill. On 31 Jul, 20:35, Pablo De Napoli [EMAIL PROTECTED] wrote: In this case, the license does not says that it is in the public domain (but that it is a copyrighted work!), but you can use it as AS IS. I think that the only condition that is imposed to us is to include the declaimer. Pablo On 7/31/07, Pablo De Napoli [EMAIL PROTECTED] wrote: Certainly no, not everything from a public institution is in the public domain. This should be analyzed case by case. In case of doubt it would be better to ask the author. Pablo On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Mon, 2007-07-30 at 17:24 -0700, Bill Hart wrote: Wow!! Excellent work indeed. In fact on 64 bit X86 systems you could actually use the 128 bit long doubles to give you a little bit more precision (I believe it only gives you 80 bits including exponent and sign, so probably 64 bit mantissa). Actually, even on my 32 bit core duo, the long double type seems to give 64 bits of precision, so using it might help a little. Do you have any idea how to check at run/compile time what the precision of a double or a long double is? It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Certainly I think you are right that any remaining optimization would be in making sure it uses no unnecessary precision. Here is a page giving information about the remainder of the series after N terms: http://mathworld.wolfram.com/PartitionFunctionP.html(see eqn 26). That's what I'm using, but I took it straight from Rademacher's paper. I think that what is needed is a careful analysis of how many terms (N) will need to be computed to get the remainder to be less than, say, r .5. Then we know that the error in each term will have to be less than (.5 - r)/N to guarantee that when we round, we will get the right answer. Also in Pari, I noted that the computation of the Psi function could dramatically slow the whole computation. I was surprised to find it figured in the runtime. It may be worthwhile checking if this is slowing things down at all. It should be computed relatively quickly if implemented correctly. The main issue was again using the minimum possible precision. I don't think that this is slowing things down much, but I could be wrong. I think that most time is spent in the actual computation of A(n,k), (this is the notation I use for the sum of exponentials -- I think that is the notation that Mathworld uses also.) If I set s(h,k) to return 1 every time, without computing anything, I only save about 20 seconds on a 2m 30s computation of p(10^9), I think. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From COPYING, Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met --Mike On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Oh well, I don't understand all this licensing stuff. So do you understand this license? What about derived works. Does that mean it is not possible to modify this library and redistribute the modified version? In particular, as a C++ library it is no use to me unmodified. There are also some functions I'd like to add, and the calling convention is the opposite to what I use. So basically I'd have to modify it to make it useful. Bill. On 31 Jul, 20:35, Pablo De Napoli [EMAIL PROTECTED] wrote: In this case, the license does not says that it is in the public domain (but that it is a copyrighted work!), but you can use it as AS IS. I think that the only condition that is imposed to us is to include the declaimer. Pablo On 7/31/07, Pablo De Napoli [EMAIL PROTECTED] wrote: Certainly no, not everything from a public institution is in the public domain. This should be analyzed case by case. In case of doubt it would be better to ask the author. Pablo On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
I believe that the IEEE standard guarantees you 80 bits (though it's only 64 bits of mantissa or something like that). The trouble is, you aren't guaranteed the IEEE standard. I've spent much time researching this, but either I didn't look at the right websites, or this stuff isn't documented well. At any rate, does autoconf allow you to check this? Unfortunately I don't know about such things. I'm an admitted ignoramous when it comes to anything to do with runtime/compile time checking. Bill. On 31 Jul, 21:33, Jonathan Bober [EMAIL PROTECTED] wrote: On Mon, 2007-07-30 at 17:24 -0700, Bill Hart wrote: Wow!! Excellent work indeed. In fact on 64 bit X86 systems you could actually use the 128 bit long doubles to give you a little bit more precision (I believe it only gives you 80 bits including exponent and sign, so probably 64 bit mantissa). Actually, even on my 32 bit core duo, the long double type seems to give 64 bits of precision, so using it might help a little. Do you have any idea how to check at run/compile time what the precision of a double or a long double is? It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Certainly I think you are right that any remaining optimization would be in making sure it uses no unnecessary precision. Here is a page giving information about the remainder of the series after N terms: http://mathworld.wolfram.com/PartitionFunctionP.html (see eqn 26). That's what I'm using, but I took it straight from Rademacher's paper. I think that what is needed is a careful analysis of how many terms (N) will need to be computed to get the remainder to be less than, say, r .5. Then we know that the error in each term will have to be less than (.5 - r)/N to guarantee that when we round, we will get the right answer. Also in Pari, I noted that the computation of the Psi function could dramatically slow the whole computation. I was surprised to find it figured in the runtime. It may be worthwhile checking if this is slowing things down at all. It should be computed relatively quickly if implemented correctly. The main issue was again using the minimum possible precision. I don't think that this is slowing things down much, but I could be wrong. I think that most time is spent in the actual computation of A(n,k), (this is the notation I use for the sum of exponentials -- I think that is the notation that Mathworld uses also.) If I set s(h,k) to return 1 every time, without computing anything, I only save about 20 seconds on a 2m 30s computation of p(10^9), I think. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Ah, that's better. Excellent. I feel much happier with this library now. Thanks. Bill. On 31 Jul, 21:38, Mike Hansen [EMAIL PROTECTED] wrote: From COPYING, Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met --Mike On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Oh well, I don't understand all this licensing stuff. So do you understand this license? What about derived works. Does that mean it is not possible to modify this library and redistribute the modified version? In particular, as a C++ library it is no use to me unmodified. There are also some functions I'd like to add, and the calling convention is the opposite to what I use. So basically I'd have to modify it to make it useful. Bill. On 31 Jul, 20:35, Pablo De Napoli [EMAIL PROTECTED] wrote: In this case, the license does not says that it is in the public domain (but that it is a copyrighted work!), but you can use it as AS IS. I think that the only condition that is imposed to us is to include the declaimer. Pablo On 7/31/07, Pablo De Napoli [EMAIL PROTECTED] wrote: Certainly no, not everything from a public institution is in the public domain. This should be analyzed case by case. In case of doubt it would be better to ask the author. Pablo On 7/31/07, Bill Hart [EMAIL PROTECTED] wrote: Yes, I've ended up using Hida and Bailey's quad-double package. Very cool. But the license just says not to use the LB name to promote any derived product. Am I right in assuming this is GPL compatible, i.e. because they are a public institution everything is automatically public domain? Bill. On 31 Jul, 06:03, [EMAIL PROTECTED] wrote: Last update, 2005? On Mon, 30 Jul 2007, Bill Hart wrote: Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
I do highly recommend this quad double library by the way. And they've implemented all manor of transcendental functions too!! The quad- doubles would give you 206 bits, even on your machine. Bill. On 31 Jul, 21:33, Jonathan Bober [EMAIL PROTECTED] wrote: On Mon, 2007-07-30 at 17:24 -0700, Bill Hart wrote: Wow!! Excellent work indeed. In fact on 64 bit X86 systems you could actually use the 128 bit long doubles to give you a little bit more precision (I believe it only gives you 80 bits including exponent and sign, so probably 64 bit mantissa). Actually, even on my 32 bit core duo, the long double type seems to give 64 bits of precision, so using it might help a little. Do you have any idea how to check at run/compile time what the precision of a double or a long double is? It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Certainly I think you are right that any remaining optimization would be in making sure it uses no unnecessary precision. Here is a page giving information about the remainder of the series after N terms: http://mathworld.wolfram.com/PartitionFunctionP.html (see eqn 26). That's what I'm using, but I took it straight from Rademacher's paper. I think that what is needed is a careful analysis of how many terms (N) will need to be computed to get the remainder to be less than, say, r .5. Then we know that the error in each term will have to be less than (.5 - r)/N to guarantee that when we round, we will get the right answer. Also in Pari, I noted that the computation of the Psi function could dramatically slow the whole computation. I was surprised to find it figured in the runtime. It may be worthwhile checking if this is slowing things down at all. It should be computed relatively quickly if implemented correctly. The main issue was again using the minimum possible precision. I don't think that this is slowing things down much, but I could be wrong. I think that most time is spent in the actual computation of A(n,k), (this is the notation I use for the sum of exponentials -- I think that is the notation that Mathworld uses also.) If I set s(h,k) to return 1 every time, without computing anything, I only save about 20 seconds on a 2m 30s computation of p(10^9), I think. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Or probably 212 actually. :-) On 31 Jul, 22:24, Bill Hart [EMAIL PROTECTED] wrote: I do highly recommend this quad double library by the way. And they've implemented all manor of transcendental functions too!! The quad- doubles would give you 206 bits, even on your machine. Bill. On 31 Jul, 21:33, Jonathan Bober [EMAIL PROTECTED] wrote: On Mon, 2007-07-30 at 17:24 -0700, Bill Hart wrote: Wow!! Excellent work indeed. In fact on 64 bit X86 systems you could actually use the 128 bit long doubles to give you a little bit more precision (I believe it only gives you 80 bits including exponent and sign, so probably 64 bit mantissa). Actually, even on my 32 bit core duo, the long double type seems to give 64 bits of precision, so using it might help a little. Do you have any idea how to check at run/compile time what the precision of a double or a long double is? It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Certainly I think you are right that any remaining optimization would be in making sure it uses no unnecessary precision. Here is a page giving information about the remainder of the series after N terms: http://mathworld.wolfram.com/PartitionFunctionP.html (see eqn 26). That's what I'm using, but I took it straight from Rademacher's paper. I think that what is needed is a careful analysis of how many terms (N) will need to be computed to get the remainder to be less than, say, r .5. Then we know that the error in each term will have to be less than (.5 - r)/N to guarantee that when we round, we will get the right answer. Also in Pari, I noted that the computation of the Psi function could dramatically slow the whole computation. I was surprised to find it figured in the runtime. It may be worthwhile checking if this is slowing things down at all. It should be computed relatively quickly if implemented correctly. The main issue was again using the minimum possible precision. I don't think that this is slowing things down much, but I could be wrong. I think that most time is spent in the actual computation of A(n,k), (this is the notation I use for the sum of exponentials -- I think that is the notation that Mathworld uses also.) If I set s(h,k) to return 1 every time, without computing anything, I only save about 20 seconds on a 2m 30s computation of p(10^9), I think. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Actually, even on my 32 bit core duo, the long double type seems to give 64 bits of precision, so using it might help a little. Do you have any idea how to check at run/compile time what the precision of a double or a long double is? Being an assembler programmer, I can say definitely that all FPU registers on x86 are 80-bit and all compilers that I know compile long double as 80-bit numbers. Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From: Alec Mihailovs [EMAIL PROTECTED] Being an assembler programmer, I can say definitely that all FPU registers on x86 are 80-bit and all compilers that I know compile long double as 80-bit numbers. From other point of view, 80-bit real gives 64-bit precision in usual sense (mantissa), while a 64-bit double has 53-bit precision in usual sense. Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Jul 31, 2007, at 18:36 , William Stein wrote:
Hi,
Just to extend this thread some more, a few remarks.
While we're extending this, here's my $0.02 canadian (:-})
On sage.math:
SAGE:
sage: time n = number_of_partitions(10^7)
CPU times: user 0.73 s, sys: 0.00 s, total: 0.73 s
Wall time: 0.73
sage: time n = number_of_partitions(10^8)
CPU times: user 8.42 s, sys: 0.00 s, total: 8.42 s
Wall time: 8.42
sage: time n = number_of_partitions(10^9)
CPU times: user 105.82 s, sys: 0.00 s, total: 105.82 s
Wall time: 105.81
Mathematica:
sage: time s=mathematica.eval('PartitionsP[10^7]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 1.90
sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 9.91
sage: time s=mathematica.eval('PartitionsP[10^9]')
CPU times: user 0.03 s, sys: 0.01 s, total: 0.04 s
Wall time: 70.43
Interesting, SAGE seems better at Mathematica for
smaller input.)
On my OS X core 2 duo 2.33Ghz laptop:
SAGE:
sage: time n=number_of_partitions(10^7)
CPU times: user 0.62 s, sys: 0.00 s, total: 0.62 s
Wall time: 0.62
sage: time n=number_of_partitions(10^8)
CPU times: user 6.99 s, sys: 0.02 s, total: 7.00 s
Wall time: 7.03
sage: time n=number_of_partitions(10^9)
CPU times: user 94.71 s, sys: 0.28 s, total: 94.99 s
Wall time: 95.56
MATHEMATICA:
sage: time s=mathematica.eval('PartitionsP[10^7]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 8.65
sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 48.08
sage: time s=mathematica.eval('PartitionsP[10^9]')
CPU times: user 0.04 s, sys: 0.00 s, total: 0.05 s
Wall time: 350.30
Yep -- Mathematica 5.2 interestingly totally sucks at
computing the number of partitions on an Intel OSX
machine... and SAGE rocks.
Checking partition count computation:
On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9:
Mathematica 5.2 (PartitionsP[10^9]:95.5115 s
Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s
Jon Bober's code (i386: 'jb 10'): 160 s
I wonder why our Core 2 Duo timings are so different?
All of the executables are i386 binaries (32-bit x86).
One oddity: Jon Bober's code produced a value ending in
'30457526857797923685688339', while Mathematica and SAGE produced one
ending in '30457526831036667979062760', so they differ.
Justin
--
Justin C. Walker, Curmudgeon at Large
Institute for the Absorption of Federal Funds
---
If it weren't for carbon-14, I wouldn't date at all.
---
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 7/31/07, Justin C. Walker [EMAIL PROTECTED] wrote: While we're extending this, here's my $0.02 canadian (:-}) Yep -- Mathematica 5.2 interestingly totally sucks at computing the number of partitions on an Intel OSX machine... and SAGE rocks. Checking partition count computation: On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9: Mathematica 5.2 (PartitionsP[10^9]:95.5115 s Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s Jon Bober's code (i386: 'jb 10'): 160 s I wonder why our Core 2 Duo timings are so different? SAGE 2.7.2.1 is different than what I just tested, since I had rebuilt mpfr with optimization turned on, which definitely improves the timings to 95seconds. I am surprised that you get 95.5115s out of mathematica, given that it takes 350 seconds with Mathematica to do the same calculation on my laptop. Let's race Mathematica's at the workshop tomorrow :-). All of the executables are i386 binaries (32-bit x86). One oddity: Jon Bober's code produced a value ending in '30457526857797923685688339', while Mathematica and SAGE produced one ending in '30457526831036667979062760', so they differ. You probably have an old version of his code? I think he emailed newer versions to me, which are what's in SAGE, and which he might not have posted to the list. -- William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Tue, 2007-07-31 at 22:16 -0700, Justin C. Walker wrote: On Jul 31, 2007, at 18:36 , William Stein wrote: Hi, Just to extend this thread some more, a few remarks. While we're extending this, here's my $0.02 canadian (:-}) [...] Checking partition count computation: On a Core 2 Duo 2.33 Mhz, computing the number of partitions of 10^9: Mathematica 5.2 (PartitionsP[10^9]:95.5115 s Sage 2.7.2.1 (number_of_partitions(10^9): 125.2 s Jon Bober's code (i386: 'jb 10'): 160 s I wonder why our Core 2 Duo timings are so different? That is puzzling. Are you sure that you have the latest version of the code? The source file should have 'Version .3' at the top. You can get the latest version, or at least that latest version minus some memory leak fixes, via hg_sage.pull(), or from http://www.math.lsa.umich.edu/~bober/partitions_c.cc The last version I posted to the list was significantly slower than what I have now (especially since I think I accidently had significant improvements commented out in the last code I posted to the list.) All of the executables are i386 binaries (32-bit x86). One oddity: Jon Bober's code produced a value ending in '30457526857797923685688339', while Mathematica and SAGE produced one ending in '30457526831036667979062760', so they differ. I am fairly confident that the answer ends in '339', since I have seen it so many times while testing my code. In fact, this is what is in the graphic on page 5 (page 27 of the pdf) of the Mathematica book. Wikipedia doesn't 'know' any congruences that p(10^9) should satisfy, so I can't use that to check the value. I'll try to figure out what the right answer should be. If your willing, try computing p(10^9 + 139) with all three. The answer should be divisible by 5, 7, and 11. (We should have been using this as a test case all along.) Incidentally, I while I was writing this email, some tests I was running just finished. The latest version of my code and the current version of sage both give the answer ending in '339'. If Mathematica and SAGE (ie PARI?) both give the same wrong answer, then there might be something fishy going on. Justin -- Justin C. Walker, Curmudgeon at Large Institute for the Absorption of Federal Funds --- If it weren't for carbon-14, I wouldn't date at all. --- --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Wow!! Excellent work indeed. In fact on 64 bit X86 systems you could actually use the 128 bit long doubles to give you a little bit more precision (I believe it only gives you 80 bits including exponent and sign, so probably 64 bit mantissa). It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Certainly I think you are right that any remaining optimization would be in making sure it uses no unnecessary precision. Here is a page giving information about the remainder of the series after N terms: http://mathworld.wolfram.com/PartitionFunctionP.html(see eqn 26). Also in Pari, I noted that the computation of the Psi function could dramatically slow the whole computation. I was surprised to find it figured in the runtime. It may be worthwhile checking if this is slowing things down at all. It should be computed relatively quickly if implemented correctly. The main issue was again using the minimum possible precision. Bill. On 28 Jul, 23:51, Jonathan Bober [EMAIL PROTECTED] wrote: I've been working on a from-scratch implementation (attached). Right now it runs a faster than Ralf Stephan's part.c, but not as fast as we would like. (And it seems to work, although I can't guarantee that right now.) On my Core Duo 1.8 ghz , it computes p(10^8) in about 17 seconds, compared to about 70 seconds for the part.c previously posted. However, it took about 270 seconds for p(10^9). (I don't have Mathematica, so I can't give that comparison.) On the other hand, I don't know how much faster sage.math is than my laptop, but since it is 64 bit, it might run the code much faster. I think that there is still a good amount of optimization that can be done to make this faster. Some things that might a lot help include better error estimates for the tail end of the series (if such estimates exist) and, in general, going over everything carefully to try to make sure that no unneeded precision is ever used. (Once the code decides that no more than 53 bits of precision are needed, it switches over to computing with C doubles, and the rest of the computation finishes instantly.) Note that this is C++ code, but it could be switched to pure C quite easily, since is doesn't actually use any real C++. part.cc 22KDownload --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 31 Jul, 01:24, Bill Hart [EMAIL PROTECTED] wrote: It would be interesting to see the time for Mathematica on a 32 bit X86 machine, since this would tell us if that is what they do. Doh! I should have read William's timings more carefully. He gives the times for a 32 bit machine. So I guess Mathematica doesn't use 80 bit long doubles on a 64 bit X86 then. Still it is an option for us. Once there is a stable version of the new code which seems to give correct results, I'll take a closer look and see if I can spot any obvious speed improvements. I can't promise anything. I suspect my fundamental mistake was not realising that you still needed quite a bit of multi-precision code for quite a few terms. In fact now that I think about it, I don't see why I thought you could compute all the s(h,k)'s using single limb arithmetic. It is the multi-precision stuff that is slowing it down, no doubt. mpfr has a 15x overhead over ordinary double precision, even at 53 bits, or so I have read. I guess there is a lot of branching to ensure the accuracy of arithmetic. Whilst that is needed for many applications, it probably isn't here. Sadly there don't seem to be any decent open source alternatives for when that accuracy is not required. I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Bill. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here: http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Hi Didier, Thanks. I also just found: http://www.nongnu.org/hpalib/ which fascinates me. Has anyone used it? Bill. On 31 Jul, 01:46, didier deshommes [EMAIL PROTECTED] wrote: 2007/7/30, Bill Hart [EMAIL PROTECTED]: I have a similar problem in some code I am currently writing. I need precisely quad precision, so mpfr is out of the question. Hi Bill, You might want to consider Yozo Hida's quaddouble C/C++ package here:http://www.cs.berkeley.edu/~yozo/ There is also a wrapper for it in SAGE. didier --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Attached is a somewhat improved version.
(I don't want to flood the mailing list with attachments, so I won't
send this file as an attachment to the list again after this.)
The main improvement is that the code is better commented and more
readable now. The secondary improvement is an increase in speed - I
think it's about 15% faster, or something like that, at computing
p(10^9) than the previous code. (Although I'm not exactly sure what the
phrase 15% faster means, and I was having some weird issues with my
laptop when I timed the old code, so I'm not sure exactly how fast it
ran.) There are no improvements in automated testing, so I'm no more
confident of this correctness of this code than I was of the old code,
but it certainly seems to work.
I'm confident that this can be sped up more, but I'm not sure that we
can reach as low as 10 seconds to compute p(10^9). It certainly might be
possible, though.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
/* Author: Jonathan Bober
* Version:.2
*
* This program computes p(n), the number of integer partitions of n, using Rademacher's
* formula. (See Hans Rademacher, On the Partition Function p(n),
* Proceedings of the London Mathematical Society 1938 s2-43(4):241-254; doi:10.1112/plms/s2-43.4.241,
* currently at
*
* http://plms.oxfordjournals.org/cgi/content/citation/s2-43/4/241
*
* if you have access.)
*
* We use the following notation:
*
* p(n) = lim_{n -- oo} t(n,N)
*
* where
*
* t(n,N) = sum_{k=1}^N a(n,k) f_n(k),
*
* where
*
* a(n,k) = sum_{h=1, (h,k) = 1}^k exp(\pi i s(h,k) - 2 \pi i h n / k)
*
* and
*
* f_n(k) = \pi sqrt{2} cosh(A_n/(sqrt{3}*k))/(B_n*k) - sinh(C_n/k)/D_n;
*
* where
*
* s(h,k) = \sum_{j=1}^{k-1}(j/k)((hj/k))
*
* A_n = sqrt{2} \pi * sqrt{n - 1/24}
* B_n = 2 * sqrt{3} * (n - 1/24)
* C_n = sqrt{2} * \pi * sqrt{n - 1.0/24.0} / sqrt{3}
* D_n = 2 * (n - 1/24) * sqrt{n - 1.0/24.0}
*
* and, finally, where ((x)) is the sawtooth function {x} - 1/2
*
* Some clever tricks are used in the computation of s(h,k), and perhaps at least
* some of these are due to Apostol. (I don't know a reference for this.)
*
* TODO: fix memory leaks
* -- All mpfr_t, mpq_t, and mpz_t variables should be cleared before the function
* mp_t(n) returns, but currently this is not done. This doesn't really matter
* when this is a standalone program, but will matter if this is called as
* a library function from another program.
*
* -- Search source code for other TODO comments.
*
* OTHER CREDITS:
*
* I looked source code written by Ralf Stephan, currently available at
*
* http://www.ark.in-berlin.de/part.c
*
* while writing this code, but didn't really make use of it. More useful
* were notes currently available at
*
* http://www.ark.in-berlin.de/part.pdf
*
* and others at
*
* http://www.math.uwaterloo.ca/~dmjackso/CO630/ptnform1.pdf
*
* Also, it is worth noting that the code for GCD(), while trivial,
* was directly copied from the NTL source code.
*
* Also, Bill Hart made some comments about ways to speed up this computation on the SAGE
* mailing list.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see http://www.gnu.org/licenses/.
*/
#include stdio.h
#include cmath
#include iostream
#include iomanip
#include limits
#include mpfr.h
#include gmp.h
const bool debug = false;
//const bool debug = true;
const bool debugs = false;
//const bool debugf = true;
const bool debugf = false;
//const bool debuga = true;
const bool debuga = false;
//const bool debugt = true;
const bool debugt = false;
using std::cout;
using std::endl;
const unsigned int min_prec = 53;
const unsigned int min_precision = 53;
// The following function can be useful for debugging in come
[sage-devel] Re: computing the number of partitions of an integer
great work, Jonathan! I've tested, and I've found the following problems: 1) part 1 hangs 2) compiling with -Wall gives this warning part.cc:865: warning: unused variable 'temp2' 3) part without arguments returns 42 Pablo On 7/29/07, Jonathan Bober [EMAIL PROTECTED] wrote: Attached is a somewhat improved version. (I don't want to flood the mailing list with attachments, so I won't send this file as an attachment to the list again after this.) The main improvement is that the code is better commented and more readable now. The secondary improvement is an increase in speed - I think it's about 15% faster, or something like that, at computing p(10^9) than the previous code. (Although I'm not exactly sure what the phrase 15% faster means, and I was having some weird issues with my laptop when I timed the old code, so I'm not sure exactly how fast it ran.) There are no improvements in automated testing, so I'm no more confident of this correctness of this code than I was of the old code, but it certainly seems to work. I'm confident that this can be sped up more, but I'm not sure that we can reach as low as 10 seconds to compute p(10^9). It certainly might be possible, though. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
I've been working on a from-scratch implementation (attached). Right now
it runs a faster than Ralf Stephan's part.c, but not as fast as we would
like. (And it seems to work, although I can't guarantee that right
now.)
On my Core Duo 1.8 ghz , it computes p(10^8) in about 17 seconds,
compared to about 70 seconds for the part.c previously posted. However,
it took about 270 seconds for p(10^9). (I don't have Mathematica, so I
can't give that comparison.) On the other hand, I don't know how much
faster sage.math is than my laptop, but since it is 64 bit, it might run
the code much faster.
I think that there is still a good amount of optimization that can be
done to make this faster. Some things that might a lot help include
better error estimates for the tail end of the series (if such estimates
exist) and, in general, going over everything carefully to try to make
sure that no unneeded precision is ever used. (Once the code decides
that no more than 53 bits of precision are needed, it switches over to
computing with C doubles, and the rest of the computation finishes
instantly.)
Note that this is C++ code, but it could be switched to pure C quite
easily, since is doesn't actually use any real C++.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
#include stdio.h
#include cmath
#include iostream
#include iomanip
#include limits
#include mpfr.h
#include gmp.h
const bool debug = false;
//const bool debug = true;
const bool debugs = false;
//const bool debugf = true;
const bool debugf = false;
//const bool debuga = true;
const bool debuga = false;
//const bool debugt = true;
const bool debugt = false;
using std::cout;
using std::endl;
const unsigned int min_prec = 53;
const unsigned int min_precision = 53;
long GCD(long a, long b);
unsigned int calc_precision(unsigned int n, unsigned int N);
//mpfr versions
mpz_t ztemp1, ztemp2, ztemp3;
mpfr_t mp_one_over_12, mp_one_over_24, mp_sqrt2, mp_sqrt3, mp_pi;
mpfr_t mp_A, mp_B, mp_C, mp_D;
mp_rnd_t round_mode = GMP_RNDN;
mpfr_t tempa1, tempa2, tempf1, tempf2, temps1, temps2, tempt1, tempt2; // temp variables for different functions, with precision set and cleared by mp_set_precision
bool mp_vars_initialized = false;
mp_prec_t mp_precision;
void mp_init(unsigned int prec, unsigned int n);
//void mp_f_precompute(unsigned int n);
void mp_f(mpfr_t result, unsigned int k);
void mp_s(mpfr_t result, unsigned int j, unsigned int q);
void mp_a(mpfr_t result, unsigned int n, unsigned int k);
void mp_t(mpfr_t result, unsigned int n);
unsigned int compute_initial_precision(unsigned int n); // computes the precision required to accurately compute p(n)
unsigned int compute_current_precision(unsigned int n, unsigned int N); // computed the precision required to
// accurately compute the tail of the rademacher series
// assuming that N terms have already been computed
double compute_remainder(unsigned int n, unsigned int N); // Gives an upper bound on the error that occurs
// when only N terms of the Rademacher series have been
// computed. NOTE: should only be called when we already know
// that the error is small (ie, compute_current_precision returns
// a small number, eg, something 32)
//low precision (double) versions of the functions:
double df_A, df_B, df_C, df_D;
void d_f_precompute(unsigned int n);
double d_f(unsigned int k);
double d_s(unsigned int h,unsigned int q);
double d_A(unsigned int n, unsigned int k);
double d_t(unsigned int n, unsigned int N);
unsigned int compute_initial_precision(unsigned int n) {
// We just want to know how many bits we will need to
// compute to get an accurate answer.
// We know that
//
// p(n) ~ exp(pi * sqrt(2n/3))/(4n sqrt(3)),
//
// so for now we are assuming that p(n) exp(pi * sqrt(2n/3))/n,
// so we need pi*sqrt(2n/3)/log(2) - log(n)/log(2) + EXTRA bits.
//
// EXTRA should depend on n, and should be something that ensures
// that the TOTAL ERROR in all computations is (something small).
// This needs to be worked out carefully. EXTRA = log(n)/log(2) + 3
// is probably good enough, and is convenient...
//
// but we really need:
//
// p(n) something
//
// to be sure that we
[sage-devel] Re: computing the number of partitions of an integer
On 7/28/07, Jonathan Bober [EMAIL PROTECTED] wrote:
I've been working on a from-scratch implementation (attached). Right now
it runs a faster than Ralf Stephan's part.c, but not as fast as we would
like. (And it seems to work, although I can't guarantee that right
now.)
Jonathan,
I've now included your code in SAGE (for 2.7.2) and made it available
(not by default) so it can be used as follows:
sage: time v=number_of_partitions(10^7, algorithm='bobber')
CPU times: user 0.86 s, sys: 0.00 s, total: 0.87 s
Wall time: 1.21
Your code agrees with Mathematica for 10^7 and 10^8, by the way:
sage: s=mathematica.eval('PartitionsP[10^7]')
sage: s = ''.join(s.replace('','').split())
sage: s == str(v)
True
sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 13.89
sage: time v=number_of_partitions(10^8, algorithm='bobber')
CPU times: user 11.94 s, sys: 0.10 s, total: 12.04 s
Wall time: 16.44
sage: s = ''.join(s.replace('','').split())
sage: s == str(v)
True
I want to release sage-2.7.2 soon, so please let me know if including part.cc
with it is OK.
-- William
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
This is fine, except that my last name is Bober, not Bobber.
(Random question, since it just caused problems for me again - it there
a way to sign up for sage-devel with a non gmail email address?)
On Sat, 2007-07-28 at 16:58 -0700, William Stein wrote:
On 7/28/07, Jonathan Bober [EMAIL PROTECTED] wrote:
I've been working on a from-scratch implementation (attached). Right now
it runs a faster than Ralf Stephan's part.c, but not as fast as we would
like. (And it seems to work, although I can't guarantee that right
now.)
Jonathan,
I've now included your code in SAGE (for 2.7.2) and made it available
(not by default) so it can be used as follows:
sage: time v=number_of_partitions(10^7, algorithm='bobber')
CPU times: user 0.86 s, sys: 0.00 s, total: 0.87 s
Wall time: 1.21
Your code agrees with Mathematica for 10^7 and 10^8, by the way:
sage: s=mathematica.eval('PartitionsP[10^7]')
sage: s = ''.join(s.replace('','').split())
sage: s == str(v)
True
sage: time s=mathematica.eval('PartitionsP[10^8]')
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 13.89
sage: time v=number_of_partitions(10^8, algorithm='bobber')
CPU times: user 11.94 s, sys: 0.10 s, total: 12.04 s
Wall time: 16.44
sage: s = ''.join(s.replace('','').split())
sage: s == str(v)
True
I want to release sage-2.7.2 soon, so please let me know if including part.cc
with it is OK.
-- William
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 7/28/07, Jonathan Bober [EMAIL PROTECTED] wrote: This is fine, except that my last name is Bober, not Bobber. Oops, I'll fix that. Thanks. (Random question, since it just caused problems for me again - it there a way to sign up for sage-devel with a non gmail email address?) I don't know. I don't think so. Let me know if this isn't true. William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer (fixes in latest pari)
This problem seems to be fixed in pari-2.4.2 (not yet realead!, code in CVS) ?numbpart((5*10535+4)) time = 16 ms. %2 = 13277569485341661441070935908285030462857350709714871167223602317834599507871542657403046634712636713000828040255875556477097762416076415269139010200121379629789951459033537585723875697364877024644629549129563676618917774238138926865673007168971398575 The CHANGES file says: Fixed 12- numbpart inaccurate (and slow). E.g numbpart(52602) off by 1 [#557]. This version of pari also includes some new functions that could be useful to Sage. For example, it includes the function polhermite and pollegendre, for computing the Hermite and Lengendre orthogonal polynomials, that currently sage computes using maxima and I think that using libpari could be much faster. Pablo On 7/26/07, Bill Hart [EMAIL PROTECTED] wrote: Pari implements it, but incorrectly: sage: number_of_partitions(5*10535+4,algorithm=pari) 132775694853416614410709359082850304628573507097 148711672236023178345995078715426574030466347126 367130008280402558755564770977624160764152691390 102001213796297899514590335375857238756973648770 246446295491295636766189177742381389268656730071 68971398574 But: sage: number_of_partitions(5*10535+4) 132775694853416614410709359082850304628573507097 148711672236023178345995078715426574030466347126 367130008280402558755564770977624160764152691390 102001213796297899514590335375857238756973648770 246446295491295636766189177742381389268656730071 68971398575 --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
I've tested it but seems to be buggy:: it works up to 156 ./part 156 p(156)= 73232243759 but for 157 gives a floating point exception error (and a gdb tracing says it is in the line 176 of the source code r2 = r0 % r1; in function g I've compiled it using gcc part.c -g -lmpfr -lm -DTEST_CODE -o part and my mpfr library is 2.2.1_p5 and gmp is 4.2.1 The code theres seems indeed to be GPL-ed, it has a copyright notice that says: (C) 2002 Ralf Stephan ([EMAIL PROTECTED]). * See part.pdf on the same path for a summary of the math. * Distributed under GPL (see gnu.org). Version 2002-12. The pdf is: http://www.ark.in-berlin.de/part.pdf Pablo On /26/07, Bill Hart [EMAIL PROTECTED] wrote: Alternatively you could just use the implementation here: http://www.ark.in-berlin.de/part.c which seems to only rely on mpfr and GMP. However, since the Pari version was based on it, I suspect it may also need correction. He does seem to adjust the precision as he goes, but I've no idea how far above or below the required precision this is. However there is no need to use a precision of 55 from a certain point on. That probably forces it to use two doubles instead of one in the floating point computations, which I think would be unnecesary. You can look on mathworld for how well the Rademacher series converges. It would probably be quicker to do a fresh implementation from scratch given the application in mind. The code there may not be GPL'd also. To check it works, one should at least look at the congruences modulo 5, 7, 11 and 13. The author of the mpfr version does seem to check them mod 11 but for very small n and only for a very few values of n. If the gcds prove to be a bottleneck, I have some code that is roughly twice as fast as GMP for computing these. Certainly the code in the mpfr version is suboptimal and needs some serious optimisation when the terms in the dedekind sum fit into a single limb, which on a 64 bit machine is probably always, when n is of a reasonable size. Another suggestion is that for sufficiently small values of h and or k, it may be quicker to compute the dedekind sum directly rather than use the gcd formula. A lookup table would also be useful for the tail end of the gcd/ dedekind sum computation. I'd be shocked if the computation for n = 10^9 couldn't be done in under 10s on sage.math. Bill. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 7/27/07, Pablo De Napoli [EMAIL PROTECTED] wrote: I've tested it but seems to be buggy:: it works up to 156 ./part 156 p(156)= 73232243759 but for 157 gives a floating point exception error (and a gdb tracing says it is in the line 176 of the source code r2 = r0 % r1; in function g I've attached a slightly modified version of part.c that seems to work (it doesn't crash and gives correct answers in all the cases I tested). I just changed the code in g slighlty so if r1 is 0, then r2 is also set to 0. I compile it using gcc part.c -O3 -g -lmpfr -lm -o part -I/home/was/sage/local/include -L/home/was/sage/local/lib/ -lgmp where /home/was/sage is the path to my SAGE install. Interestingly, when I time part.c it is not much better than CVS pari (see below). Bill Hart wrote: This, with the other ideas I gave above, will definitely let us beat Mathematica convincingly, as a simple back of the envelope calculation shows. It will have the added benefit of allowing us to beat Magma at something. Magma is already terrible at computing NumberOfPartitions, though not as bad as GAP (which is even worse): [EMAIL PROTECTED]:~$ magma Magma V2.13-5 Fri Jul 27 2007 10:54:07 on ubuntu [Seed = 1720324423] Type ? for help. Type Ctrl-D to quit. time k := NumberOfPartitions(10^5); Time: 3.070 With the latest CVS pari: ? gettime; n=numbpart(10^5); gettime/1000.0 %1 = 0.01200 ? gettime; n=numbpart(10^6); gettime/1000.0 %2 = 0.2080 ? gettime; n=numbpart(10^7); gettime/1000.0 %3 = 6.973 ? gettime; n=numbpart(10^8); gettime/1000.0 With part.c: I use time ./part 10 out and record system + user time: 10^60.04 seconds 10^74.584 seconds 10^8 47.483 seconds So, to clarify or summarize, Bill is there one thing that we could change in part.c that would speed it up? I realized you sort of answered this before, but the sequence of previous emails yesterday was rather long and may have got confused. Thanks, william The code theres seems indeed to be GPL-ed, it has a copyright notice that says: (C) 2002 Ralf Stephan ([EMAIL PROTECTED]). * See part.pdf on the same path for a summary of the math. * Distributed under GPL (see gnu.org). Version 2002-12. The pdf is: http://www.ark.in-berlin.de/part.pdf Pablo On /26/07, Bill Hart [EMAIL PROTECTED] wrote: Alternatively you could just use the implementation here: http://www.ark.in-berlin.de/part.c which seems to only rely on mpfr and GMP. However, since the Pari version was based on it, I suspect it may also need correction. He does seem to adjust the precision as he goes, but I've no idea how far above or below the required precision this is. However there is no need to use a precision of 55 from a certain point on. That probably forces it to use two doubles instead of one in the floating point computations, which I think would be unnecesary. You can look on mathworld for how well the Rademacher series converges. It would probably be quicker to do a fresh implementation from scratch given the application in mind. The code there may not be GPL'd also. To check it works, one should at least look at the congruences modulo 5, 7, 11 and 13. The author of the mpfr version does seem to check them mod 11 but for very small n and only for a very few values of n. If the gcds prove to be a bottleneck, I have some code that is roughly twice as fast as GMP for computing these. Certainly the code in the mpfr version is suboptimal and needs some serious optimisation when the terms in the dedekind sum fit into a single limb, which on a 64 bit machine is probably always, when n is of a reasonable size. Another suggestion is that for sufficiently small values of h and or k, it may be quicker to compute the dedekind sum directly rather than use the gcd formula. A lookup table would also be useful for the tail end of the gcd/ dedekind sum computation. I'd be shocked if the computation for n = 10^9 couldn't be done in under 10s on sage.math. Bill. -- William Stein Associate Professor of Mathematics University of Washington http://www.williamstein.org --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~--- /* Implementation of the unrestricted partition function using mpfr (see * mpfr.org). (C) 2002 Ralf Stephan ([EMAIL PROTECTED]). * See part.pdf on the same path for a summary of the math. * Distributed under GPL (see gnu.org). Version 2002-12. * * The unrestricted partition function counts the nonnegative integer * solutions to a+2b+3c+...=n.
[sage-devel] Re: computing the number of partitions of an integer
Your version did work for me! Pablo On 7/27/07, William Stein [EMAIL PROTECTED] wrote: On 7/27/07, Pablo De Napoli [EMAIL PROTECTED] wrote: I've tested it but seems to be buggy:: it works up to 156 ./part 156 p(156)= 73232243759 but for 157 gives a floating point exception error (and a gdb tracing says it is in the line 176 of the source code r2 = r0 % r1; in function g I've attached a slightly modified version of part.c that seems to work (it doesn't crash and gives correct answers in all the cases I tested). I just changed the code in g slighlty so if r1 is 0, then r2 is also set to 0. I compile it using gcc part.c -O3 -g -lmpfr -lm -o part -I/home/was/sage/local/include -L/home/was/sage/local/lib/ -lgmp where /home/was/sage is the path to my SAGE install. Interestingly, when I time part.c it is not much better than CVS pari (see below). Bill Hart wrote: This, with the other ideas I gave above, will definitely let us beat Mathematica convincingly, as a simple back of the envelope calculation shows. It will have the added benefit of allowing us to beat Magma at something. Magma is already terrible at computing NumberOfPartitions, though not as bad as GAP (which is even worse): [EMAIL PROTECTED]:~$ magma Magma V2.13-5 Fri Jul 27 2007 10:54:07 on ubuntu [Seed = 1720324423] Type ? for help. Type Ctrl-D to quit. time k := NumberOfPartitions(10^5); Time: 3.070 With the latest CVS pari: ? gettime; n=numbpart(10^5); gettime/1000.0 %1 = 0.01200 ? gettime; n=numbpart(10^6); gettime/1000.0 %2 = 0.2080 ? gettime; n=numbpart(10^7); gettime/1000.0 %3 = 6.973 ? gettime; n=numbpart(10^8); gettime/1000.0 With part.c: I use time ./part 10 out and record system + user time: 10^60.04 seconds 10^74.584 seconds 10^8 47.483 seconds So, to clarify or summarize, Bill is there one thing that we could change in part.c that would speed it up? I realized you sort of answered this before, but the sequence of previous emails yesterday was rather long and may have got confused. Thanks, william The code theres seems indeed to be GPL-ed, it has a copyright notice that says: (C) 2002 Ralf Stephan ([EMAIL PROTECTED]). * See part.pdf on the same path for a summary of the math. * Distributed under GPL (see gnu.org). Version 2002-12. The pdf is: http://www.ark.in-berlin.de/part.pdf Pablo On /26/07, Bill Hart [EMAIL PROTECTED] wrote: Alternatively you could just use the implementation here: http://www.ark.in-berlin.de/part.c which seems to only rely on mpfr and GMP. However, since the Pari version was based on it, I suspect it may also need correction. He does seem to adjust the precision as he goes, but I've no idea how far above or below the required precision this is. However there is no need to use a precision of 55 from a certain point on. That probably forces it to use two doubles instead of one in the floating point computations, which I think would be unnecesary. You can look on mathworld for how well the Rademacher series converges. It would probably be quicker to do a fresh implementation from scratch given the application in mind. The code there may not be GPL'd also. To check it works, one should at least look at the congruences modulo 5, 7, 11 and 13. The author of the mpfr version does seem to check them mod 11 but for very small n and only for a very few values of n. If the gcds prove to be a bottleneck, I have some code that is roughly twice as fast as GMP for computing these. Certainly the code in the mpfr version is suboptimal and needs some serious optimisation when the terms in the dedekind sum fit into a single limb, which on a 64 bit machine is probably always, when n is of a reasonable size. Another suggestion is that for sufficiently small values of h and or k, it may be quicker to compute the dedekind sum directly rather than use the gcd formula. A lookup table would also be useful for the tail end of the gcd/ dedekind sum computation. I'd be shocked if the computation for n = 10^9 couldn't be done in under 10s on sage.math. Bill. -- William Stein Associate Professor of Mathematics University of Washington http://www.williamstein.org --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 7/27/07, Pablo De Napoli [EMAIL PROTECTED] wrote: Your version did work for me! Pablo Great. Is there any chance you could work on any of the optimization or timing ideas Bill Hart suggested in the previous post? I won't be able to work on this myself for a while. -- William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On Thu, 2007-07-26 at 13:10 -0700, Bill Hart wrote:
Perhaps if one had a fast way of evaluating the Dedekind sum, one
might have a chance.
Bill.
I think this gives a faster way to compute it:
Write the sum as
s(h,k) = sum_{j=1}^{k-1} j/k [hj/k - floor(hj/k) - 1/2]
(This isn't strictly correct in general, but in our case hj/k will never
be an integer, so we are ok.)
Then if we separate this into three different sums and use some simple
summation formulas, we get
s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).
(To compute the floor in the inner sum, you can just use integer
division.)
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Here we go. Apparently g(h,q) = q*s(h,q) where s(h,q) is a dedekind
sum.
Then for h q if r_0, r_1, ..., r_{n+1} are the remainders occurring
in the Euclidean algorithm when computing gcd(h,q) (of which there
should be about log q) then:
s(h,q) = 1/12 * sum_{i=1}^{n+1}((-1)^(i+1) (r_i^2 + r_{i-1}^2 + 1) /
(r_i * r_{i-1}) - ((-1)^n +1)/8)
This, with the other ideas I gave above, will definitely let us beat
Mathematica convincingly, as a simple back of the envelope calculation
shows. It will have the added benefit of allowing us to beat Magma at
something.
Bill.
On 27 Jul, 00:12, Bill Hart [EMAIL PROTECTED] wrote:
OK, apparently the Dedekind sums should be done using an algorithm due
to Apostol. I haven't tracked it down yet, but this is clearly what we
need.
Bill.
On 26 Jul, 23:37, Bill Hart [EMAIL PROTECTED] wrote:
On 26 Jul, 23:18, Jonathan Bober [EMAIL PROTECTED] wrote:
s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).
(To compute the floor in the inner sum, you can just use integer
division.)
Yes, this will be faster. Of course integer division is not faster
than floating point division, but since we are doing a sum from j = 1
to k-1 we only need to know when floor(hj/k) increases by 1. This is
simply a matter of adding h each iteration and seeing if we have gone
over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
and so on.
This gets the inner computation down to about 8 cycles or so on
average. Not enough to beat Mathematica though, unless I made a
mistake in my back of the envelope computation somewhere.
Bill.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
On 26 Jul, 23:18, Jonathan Bober [EMAIL PROTECTED] wrote:
s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).
(To compute the floor in the inner sum, you can just use integer
division.)
Yes, this will be faster. Of course integer division is not faster
than floating point division, but since we are doing a sum from j = 1
to k-1 we only need to know when floor(hj/k) increases by 1. This is
simply a matter of adding h each iteration and seeing if we have gone
over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
and so on.
This gets the inner computation down to about 8 cycles or so on
average. Not enough to beat Mathematica though, unless I made a
mistake in my back of the envelope computation somewhere.
Bill.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
There are two tricks I can see that are required to make this faster. Firstly notice that L(n,q)*Psi(n,q) is relatively small once q gets beyond a certain point. Thus L(n,q)*Psi(n,q) can be computed using ordinary double precision for q greater than some very small limit (depending on n). If one knew how fast this series converges (which must have been worked out somewhere), one could even progressively reduce the precision as the computation went, for even greater speed. For example the following Pari script should compute all but the last 150 or so digits of P(10^9) quite quickly: Psi(n, q) = local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b); (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c)) L(n, q) = if(q==1,1,sum(h=1,q-1,if(gcd(h,q)1,0,cos((g(h,q)-2*h*n)*Pi/ q g(h, q) = if(q3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2))) n=14 \p36000 a1 = L(n,1)*Psi(n,1); \p18000 a2 = L(n,2)*Psi(n,2); \p12000 a3 = L(n,3)*Psi(n,3); \p9000 a4 = L(n,4)*Psi(n,4); \p8000 a5 = L(n,5)*Psi(n,5); \p6000 a6 = L(n,6)*Psi(n,6); \p6000 a7 = L(n,7)*Psi(n,7); \p5000 a8 = L(n,8)*Psi(n,8); \p4000 a9 = sum(y=9,14,L(n,y)*Psi(n,y)); \p3000 a10 = sum(y=15,17,L(n,y)*Psi(n,y)); \p2000 a11 = sum(y=18,34,L(n,y)*Psi(n,y)); \p1000 a12 = sum(y=35,300,L(n,y)*Psi(n,y)); round(a1+a2+a3+a4+a5+a6+a7+a8+a9+a10+a11+a12) The second trick is that computing gcd(h,q) is costly. One should factor q and then sieve for all h not coprime to q in short cache friendly segments, especially for very large n. But I don't believe Mathematica uses this Rademacher series thing anyway. If you look at the inner most loop, that computation involving frac must take at least 80 cycles in double precision. But for n = 10^9, that expression must get evaluated about 2*10^11 times. That's about 1.6*10^13 cycles, which on sage.math has go to take about 1s. Even if I'm out by a factor of 10 with the cycles, mathematica clearly doesn't do this. Perhaps if one had a fast way of evaluating the Dedekind sum, one might have a chance. Bill. --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
William Stein [EMAIL PROTECTED] writes: QUESTIONS: Why is Mathematica about 10 times faster than PARI at this? What are the best ways to compute the number of partitions of n? Is it a calculation involving fast arithmetic with dense polynomials of large degree, which would be best done using the upcoming FLINT library or NTL? Please correct me if I'm crazy, but isn't there an asymptotic formula due to Hardy and Rademacher that can evaluate $P(n)$ to a very high accuracy very quickly? Surely both of these packages implement such a convergent series approach, and it is possible that SAGE has faster real arithmetic and could do this faster. Again, I may be completely incorrect -- please let me know. Nick --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
The fastest way I know of to compute the number of integer partitions is to use the Hardy-Ramanujan-Rademacher formula ( see Rademacher series on http://en.wikipedia.org/wiki/Integer_partition ), and I'm pretty sure this is what PARI uses. From modules/part.c: * This program is basically the implementation of the script * * Psi(n, q) = local(a, b, c); a=sqrt(2/3)*Pi/q; b=n-1/24; c=sqrt(b); * (sqrt(q)/(2*sqrt(2)*b*Pi))*(a*cosh(a*c)-(sinh(a*c)/c)) * L(n, q) = if(q==1,1,sum(h=1,q-1,if(gcd(h,q)1,0,cos((g(h,q)-2*h*n)*Pi/q * g(h, q) = if(q3,0,sum(k=1,q-1,k*(frac(h*k/q)-1/2))) * part(n) = round(sum(q=1,5 + 0.24*sqrt(n),L(n,q)*Psi(n,q))) I'm not sure where the bottleneck in PARI is since I can't imagine Mathematica uses a different method to compute the number of partitions. --Mike P.S. I'm going to push hard this weekend to get the combinatorics stuff I've been working on to you. On 7/26/07, William Stein [EMAIL PROTECTED] wrote: On 7/25/07, Timothy Clemans [EMAIL PROTECTED] wrote: How do you in general find out how say Magma or Mathematica does something? Short answer: it is often virtually impossible -- that's one of the main reasons for the existence of SAGE. Read this page for the official Mathematica stance on this sort of question: http://reference.wolfram.com/mathematica/tutorial/WhyYouDoNotUsuallyNeedToKnowAboutInternals.html -- William --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Alternatively you could just use the implementation here:
http://www.ark.in-berlin.de/part.c
which seems to only rely on mpfr and GMP. However, since the Pari
version was based on it, I suspect it may also need correction.
He does seem to adjust the precision as he goes, but I've no idea how
far above or below the required precision this is. However there is no
need to use a precision of 55 from a certain point on. That probably
forces it to use two doubles instead of one in the floating point
computations, which I think would be unnecesary. You can look on
mathworld for how well the Rademacher series converges.
It would probably be quicker to do a fresh implementation from scratch
given the application in mind. The code there may not be GPL'd also.
To check it works, one should at least look at the congruences modulo
5, 7, 11 and 13. The author of the mpfr version does seem to check
them mod 11 but for very small n and only for a very few values of n.
If the gcds prove to be a bottleneck, I have some code that is roughly
twice as fast as GMP for computing these. Certainly the code in the
mpfr version is suboptimal and needs some serious optimisation when
the terms in the dedekind sum fit into a single limb, which on a 64
bit machine is probably always, when n is of a reasonable size.
Another suggestion is that for sufficiently small values of h and or
k, it may be quicker to compute the dedekind sum directly rather than
use the gcd formula.
A lookup table would also be useful for the tail end of the gcd/
dedekind sum computation.
I'd be shocked if the computation for n = 10^9 couldn't be done in
under 10s on sage.math.
Bill.
On 27 Jul, 00:39, Bill Hart [EMAIL PROTECTED] wrote:
Here we go. Apparently g(h,q) = q*s(h,q) where s(h,q) is a dedekind
sum.
Then for h q if r_0, r_1, ..., r_{n+1} are the remainders occurring
in the Euclidean algorithm when computing gcd(h,q) (of which there
should be about log q) then:
s(h,q) = 1/12 * sum_{i=1}^{n+1}((-1)^(i+1) (r_i^2 + r_{i-1}^2 + 1) /
(r_i * r_{i-1}) - ((-1)^n +1)/8)
This, with the other ideas I gave above, will definitely let us beat
Mathematica convincingly, as a simple back of the envelope calculation
shows. It will have the added benefit of allowing us to beat Magma at
something.
Bill.
On 27 Jul, 00:12, Bill Hart [EMAIL PROTECTED] wrote:
OK, apparently the Dedekind sums should be done using an algorithm due
to Apostol. I haven't tracked it down yet, but this is clearly what we
need.
Bill.
On 26 Jul, 23:37, Bill Hart [EMAIL PROTECTED] wrote:
On 26 Jul, 23:18, Jonathan Bober [EMAIL PROTECTED] wrote:
s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).
(To compute the floor in the inner sum, you can just use integer
division.)
Yes, this will be faster. Of course integer division is not faster
than floating point division, but since we are doing a sum from j = 1
to k-1 we only need to know when floor(hj/k) increases by 1. This is
simply a matter of adding h each iteration and seeing if we have gone
over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
and so on.
This gets the inner computation down to about 8 cycles or so on
average. Not enough to beat Mathematica though, unless I made a
mistake in my back of the envelope computation somewhere.
Bill.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
I just found this on http://reference.wolfram.com/mathematica/note/SomeNotesOnInternalImplementation.html: PartitionsP[n] uses Euler's pentagonal formula for small n, and the non-recursive Hardy-Ramanujan-Rademacher method for larger n. --Mike On 7/26/07, Alec Mihailovs [EMAIL PROTECTED] wrote: From: Mike Hansen [EMAIL PROTECTED] I'm not sure where the bottleneck in PARI is since I can't imagine Mathematica uses a different method to compute the number of partitions. I don't know what is used in the latest Mathematica version, but originally NumberOfPartitions function in the Combinatorica package used the recursion with pentagonal numbers, see http://www.cs.uiowa.edu/~sriram/Combinatorica/NewCombinatorica.m Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
Pari implements it, but incorrectly:
sage: number_of_partitions(5*10535+4,algorithm=pari)
132775694853416614410709359082850304628573507097
148711672236023178345995078715426574030466347126
367130008280402558755564770977624160764152691390
102001213796297899514590335375857238756973648770
246446295491295636766189177742381389268656730071
68971398574
But:
sage: number_of_partitions(5*10535+4)
132775694853416614410709359082850304628573507097
148711672236023178345995078715426574030466347126
367130008280402558755564770977624160764152691390
102001213796297899514590335375857238756973648770
246446295491295636766189177742381389268656730071
68971398575
At least GAP returns the right answer!! But after some time:
sage: number_of_partitions(5*10535+4)
---
type 'exceptions.RuntimeError' Traceback (most recent call
last)
/home/wbhart/flint/trunk/ipython console in module()
/home/was/s/local/lib/python2.5/site-packages/sage/combinat/
combinat.py in number_of_partitions(n, k, algorithm)
1245 if algorithm == 'gap':
1246 if k==None:
- 1247 ans=gap.eval(NrPartitions(%s)%(ZZ(n)))
1248 else:
1249 ans=gap.eval(NrPartitions(%s,%s)%(ZZ(n),ZZ(k)))
/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/gap.py
in eval(self, x, newlines, strip)
298 if len(x) == 0 or x[len(x) - 1] != ';':
299 x += ';'
-- 300 s = Expect.eval(self, x)
301 if newlines:
302 return s
/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/
expect.py in eval(self, code, strip, **kwds)
519 code = code.strip()
520 try:
-- 521 return '\n'.join([self._eval_line(L, **kwds) for L
in code.split('\n') if L != ''])
522 except KeyboardInterrupt:
523 self._keyboard_interrupt()
/home/was/s/local/lib/python2.5/site-packages/sage/interfaces/gap.py
in _eval_line(self, line, allow_use_file, wait_for_prompt)
482 return ''
483 else:
-- 484 raise RuntimeError, message
485
486 except KeyboardInterrupt:
type 'exceptions.RuntimeError': Unexpected EOF from Gap executing
NrPartitions(52679);
Bill.
On 26 Jul, 07:38, Nick Alexander [EMAIL PROTECTED] wrote:
William Stein [EMAIL PROTECTED] writes:
QUESTIONS: Why is Mathematica about 10 times faster than PARI
at this? What are the best ways to compute the number
of partitions of n? Is it a calculation involving fast arithmetic
with dense polynomials of large degree, which would be best done
using the upcoming FLINT library or NTL?
Please correct me if I'm crazy, but isn't there an asymptotic formula
due to Hardy and Rademacher that can evaluate $P(n)$ to a very high
accuracy very quickly? Surely both of these packages implement such a
convergent series approach, and it is possible that SAGE has faster
real arithmetic and could do this faster.
Again, I may be completely incorrect -- please let me know.
Nick
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
From: Mike Hansen [EMAIL PROTECTED] I'm not sure where the bottleneck in PARI is since I can't imagine Mathematica uses a different method to compute the number of partitions. I don't know what is used in the latest Mathematica version, but originally NumberOfPartitions function in the Combinatorica package used the recursion with pentagonal numbers, see http://www.cs.uiowa.edu/~sriram/Combinatorica/NewCombinatorica.m Alec --~--~-~--~~~---~--~~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to [EMAIL PROTECTED] For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/ -~--~~~~--~~--~--~---
[sage-devel] Re: computing the number of partitions of an integer
OK, apparently the Dedekind sums should be done using an algorithm due
to Apostol. I haven't tracked it down yet, but this is clearly what we
need.
Bill.
On 26 Jul, 23:37, Bill Hart [EMAIL PROTECTED] wrote:
On 26 Jul, 23:18, Jonathan Bober [EMAIL PROTECTED] wrote:
s(h,k) = h(k - 1)(2k - 1)/(6k) - (k-1)/4 -
(1/k) sum_{j=1}^{k-1} j*floor(hj/k).
(To compute the floor in the inner sum, you can just use integer
division.)
Yes, this will be faster. Of course integer division is not faster
than floating point division, but since we are doing a sum from j = 1
to k-1 we only need to know when floor(hj/k) increases by 1. This is
simply a matter of adding h each iteration and seeing if we have gone
over k (subtracting k if we do). If so, floor(hj/k) has increased by 1
and so on.
This gets the inner computation down to about 8 cycles or so on
average. Not enough to beat Mathematica though, unless I made a
mistake in my back of the envelope computation somewhere.
Bill.
--~--~-~--~~~---~--~~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~--~~~~--~~--~--~---
