Search code examples
calgorithmgmpgreatest-common-divisor

GCD computation issues with GNU MP Library


I have a question about GNU MP, could you please help me how to proceed with that. I using "The GNU Multiple Precision Arithmetic Library" Edition 5.1.1 on Windows. (MinGW\gcc + MSYS)

There is an mpz_gcd function exist to calculate "gcd" of two integers.

void mpz_gcd (mpz_t rop, mpz_t op1, mpz_t op2);

As far as I get from documentation, there are couple of algorithms implemented in GNU MP for computation of Greatest Common Divisor. Among them:

  • Binary GCD
  • Lehmer’s algorithm
  • Subquadratic GCD

The algorithm that is used seems to be selected automatically, based on the input size of integers.

Currently, the binary algorithm is used for GCD only when N < 3.

For inputs larger than GCD_DC_THRESHOLD, GCD is computed via the HGCD (Half GCD) function as a generalization to Lehmer’s algorithm.

So, I guess there are at least three different approaches to get gcd(a, b). The main problem for me: I want to specify which algorithm to use by myself. I would compare time execution of these algorithms on the random large input (i.e. 10^5 bits) to find out some common trends: what is that point where using "Binary GCD" becomes worse than "Lehmer's method", is "HGCD-Lehmer generalization" is really better than straightforward Lehmer, etc.

Is there any simple way to specify the algorithm you want to use? Any way to pull out this algorithm from the library, any way to modify some "#define" variables. Is it possible to do something like I want without library recompilation? I'm just beginner there and I don't feel able to figure out all sort of things inside the library.

P.S. Probably somebody would be interested what's will come out of that. I've got some code on the github: https://github.com/int000h/gcd_gcc


Solution

  • This is a good time to be reading source code. GMP is open source—take advantage of that!

    In mpn/generic/gcd.c you'll find the function which selects the GCD algorithm (this is actually a public function, it appears in the documentation):

    mp_size_t
    mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
    {
        ...
        if (ABOVE_THRESHOLD (n, GCD_DC_THRESHOLD)) {
            ...
    

    You can see that there are three main branches to the function, each ending with a return statement. Each branch corresponds to a different GCD algorithm. You can copy and paste the code into your own application and modify it so you can specify exactly which algorithm you want. Tips:

    • You can get rid of the #ifdefs. Assume TUNE_GCD_P is not defined.

    • This is an mpn_* function instead of an mpz_* function. It's lower-level: you have to explicitly allocate space for outputs, for example. You may also wish to copy the code from the higher-level function, mpz_gcd().

    • You'll need to extract prototypes for the internal functions, like mpn_hgcd_matrix_adjust(). Just copy the prototypes out of the GMP source code. Don't worry, internal functions appear to be exported from the shared library (they generally shouldn't be, but they are, so you're fine).

    No need to recompile the library, but you will need to do a little bit of work here.