---------- Forwarded message ---------- From: William Hart <[EMAIL PROTECTED]> Date: Jul 27, 2007 11:58 AM Subject: Re: [sage-devel] Re: computing the number of partitions of an integer To: William Stein <[EMAIL PROTECTED]>
Some things to try (in order of importance): * how long does it take to compute the first few (say 100 odd) terms. It should be nearly instant. If not the precision is held too high for too long. * much of the time is surely taken up by GMP in the s(h,k) computation. One shouldn't use GMP for that. It really needs to be optimally fast. I'd code a C version, though to make any difference you probably need to use GMP's longlong.h to get 128 bit arithmetic. Obviously you still need to use GMP outside that range. Frankly I can't see where the time is going in there. Maybe I've just underestimated how long it takes to do one of these, or underestimated how many of them need to be done. Time how long each s(h, k) is taking when the series has tailed off a bit (say under 150 digits for n= 10^9). My bet is all the time is being taken up by an inefficient implementation of these. * time how much time is spent in gcds. I have code that will double the speed of those if they are a bottleneck. But I sincerely doubt it. Once you know where the bottleneck is, you can look more carefully at that section of code and I might have some more ideas. BTW google groups will not accept postings from this email address. Bill. --- 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.01200000000000000000000000000 > ? gettime; n=numbpart(10^6); gettime/1000.0 > %2 = 0.2080000000000000000000000000 > ? gettime; n=numbpart(10^7); gettime/1000.0 > %3 = 6.973000000000000000000000000 > ? gettime; n=numbpart(10^8); gettime/1000.0 > > > With part.c: > > I use time ./part 100000 > out and record system + > user time: > 10^6 0.04 seconds > 10^7 4.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. > > > > --~--~---------~--~----~------------~-------~--~----~ 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/ -~----------~----~----~----~------~----~------~--~---
