I hacked some hooks into gmp.h to export mpz_sgcd and mpz_ngcd. Both of these functions run in the right amount of time.
I think the mpn_gcd is just broken. I don't see why we shouldn't just use ngcd for now. Don't quote me on this, but one possibility is that mpn_gcd still implements Weber's algorithm, but a subordinate part of it is using the ngcd_lehmer algorithm. This appears to make sense as the code runs some of the existing code if a certain threshold is passed, but always runs the ngcd_lehmer code. That would seem to imply that the algorithm was originally designed to break the problem down somewhat (I don't know much about Weber's algorithm, but it does some kind of k-ary reduction and an nmod step, then feeds it to a simpler gcd algorithm as a basecase. It is the binary gcd algorithm, used as a simpler gcd algorithm which has been replaced by ngcd_lehmer in Moller's patches for mpn_gcd. See Weber's paper for details: http://raider.muc.edu/cs/cs492/samples/sce_blog/linked_files/Weber_Presentation_Paper.pdf Of course the Weber algorithm supposedly runs about 5 times faster than a binary gcd algorithm and so the idea may be that for smaller inputs Weber + ngcd will be faster than Weber + binary and even faster than ngcd. But for large inputs ngcd will win due to the subquadratic asymptotics. Anyhow, there appears to be no reason why one can't just call mpn_ngcd directly for large enough inputs, bypassing the Weber reduction. As you've noted though, over a certain range of small inputs, Weber + ngcd is slower than Weber + binary gcd. Of course I may have this all mixed up. I have only taken a brief look. If we only end up using ngcd (quite possible), we should get rid of all the others (sgcd, bgcd, etc) as Moller reports in his paper that they are all inferior to his ngcd. Bill. 2008/12/21 Bill Hart <goodwillh...@googlemail.com>: > There is a graph of expected performance on page 16 of Moller's paper: > > http://www.lysator.liu.se/~nisse/archive/S0025-5718-07-02017-0.pdf > > It looks like the ngcd algorithm is supposed to be the fastest. But it > is still clear from the graph that at 1 million bits ngcd is > dramatically faster than the accelerated gcd, like 10 times faster. > Something is wrong somewhere. > > Perhaps the version of the patches on Moller's site is broken. But > presumably this is the code he used to construct the timings for his > paper. > > Bill. > > 2008/12/21 Bill Hart <goodwillh...@googlemail.com>: >> I see why it is faster. There is a #define in the accelerated gcd >> algorithm switching out some of the old code. Now a function called >> mpn_ngcd_lehmer is called. >> >> This is one of Moller's functions, presumably implementing his "new >> half gcd algorithm". According to his paper this is supposed to be >> faster than the half gcd algorithm of Schoenhage. >> >> Now I have no idea what is going on. >> >> Bill. >> >> 2008/12/21 Bill Hart <goodwillh...@googlemail.com>: >>> The function mpz_gcd is defined in the file /mpz/gcd. In the latest >>> revision of MPIR it invokes a macro GCD_BODY(mpn_gcd). This macro is >>> just a generic wrapper for an mpn_gcd like function. It can be invoked >>> with any function which observes the syntax of mpn_gcd as a parameter. >>> But currently mpz_gcd is invoked with parameter mpn_gcd. >>> >>> mpn_gcd is defined in the file /mpn/generic/gcd.c. Currently it >>> implements the algorithm known as the accelerated gcd algorithm: >>> >>> K. Weber, The accelerated integer GCD algorithm, ACM Transactions on >>> Mathematical Software, v. 21 (March), 1995, pp. 111-122. >>> >>> as implemented by Ken Weber in 2002. It does *NOT* use Moller's >>> patches or a half gcd algorithm. >>> >>> Instead, mpn_gcd should be modified to actually call mpn_hgcd. >>> >>> As you say, it may be possible to have some kind of tuned wrapper >>> which chooses the best algorithm depending on the input. >>> >>> I don't know why mpz_gcd is currently faster than it used to be. >>> Perhaps some other difference in our code base accounts for that. >>> >>> The whole point is, Moller's algorithm is a half gcd algorithm. This >>> is asymptotically faster than the accelerated gcd algorithm (it is >>> subquadratic). This means that as the sizes increase the distance >>> between the before and after timings should become more and more >>> dramatic. >>> >>> One person that might be able to help is Moller. He seemed nice enough >>> when I spoke to him via email. I don't know what the intended usage of >>> his functions is, and I haven't checked to see how well documented >>> they are. Presumably one wouldn't call hgcd when one of the inputs is >>> very small and the other very large. There also appear to be numerous >>> other functions in his files and I haven't looked at his documentation >>> or paper to see what situations they are supposed to be advantageous >>> in. >>> >>> I kind of assumed that you were sorting all this out. >>> >>> Bill. >>> >>> 2008/12/21 Jason Martin <jason.worth.mar...@gmail.com>: >>>> >>>> What about the figures? I see that the new mpn_gcd is faster than the >>>> GMP v. 4.2 mpn_gcd. >>>> >>>> I'm confused about what you want me to do. I thought the assignment >>>> was to put Moller's v2.1 patches into MPIR. I did that. The patches >>>> are there and they DO show a difference with a high enough limb count. >>>> >>>> I can't use the Moller patches on the GMP webpage because they are GPL v. >>>> 3. >>>> >>>> Jason Worth Martin >>>> Asst. Professor of Mathematics >>>> http://www.math.jmu.edu/~martin >>>> >>>> >>>> >>>> On Sat, Dec 20, 2008 at 7:12 PM, Bill Hart <goodwillh...@googlemail.com> >>>> wrote: >>>>> >>>>> You need to look at the figures on the GMP website and look at the >>>>> asymptotics. That has nothing to do with Magma. >>>>> >>>>> Bill. >>>>> >>>>> 2008/12/21 Jason Martin <jason.worth.mar...@gmail.com>: >>>>>> >>>>>>> By the way, the code you posted has a variable missing in the final >>>>>>> printf. It always returns 0.00. >>>>>> >>>>>> Sorry, that was a version control issue. The numbers I posted below >>>>>> are for correct code (with the "result" variable printed. >>>>>> >>>>>> I don't get the error message that you are reporting when building. >>>>>> I'll try building on sage.math and see what happens. >>>>>> >>>>>> I'm also confused about what exactly you're saying the problem is. My >>>>>> testing shows that for big enough numbers, the Moller patches in MPIR >>>>>> are faster than the un-patched MPIR. I'm not comparing it against >>>>>> Magma because I don't have any basis for the comparison. >>>>>> >>>>>> Perhaps you should send me the exact numbers you're talking about and >>>>>> I'll try them with: >>>>>> >>>>>> MPIR rev 1487 (pre-patch) >>>>>> MPIR current >>>>>> and >>>>>> GMP 4.2.4 >>>>>> >>>>>> --jason >>>>>> >>>>>> > >>>>>> >>>>> >>>>> > >>>>> >>>> >>>> >>>> >>>> >>> >> > --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "mpir-devel" group. To post to this group, send email to mpir-devel@googlegroups.com To unsubscribe from this group, send email to mpir-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/mpir-devel?hl=en -~----------~----~----~----~------~----~------~--~---
