GCD algorithms for arbitrary-precision arithmetic

Question

I am completely stuck with this question, so I am looking for any help.

I think everybody knows about basic GCD computation algorithms like binary or euclidean GCD. It is not a problem to implement such a method to compute two single-precision numbers. Actually, it's just about couple of strokes.

I need to have this method implemented (in C language) for multiple-precision numbers (more than 10^5 bits). There are couple of GNU libraries available (GNU MP, MPFR, MPIR) and they have means to define multiple-precision numbers and make an actions on them. It looks like one multiple-precision number stored in memory as couple of single-precision parts aka "limbs".

The have some methods implemented for finding gcd(a, b), but they, actually, difficult to use for my needs. Binary method for GCD computation used only when a and b contains exactly two limbs. HGCD method used when min(a,b) contains more than (i.e. 630) limbs, etc. I find difficult to figure out, how any of these methods could be expanded for using with any length of a and b. I've also found that different version of GNU libraries contains different versions and methods of GCD algorithms.

Question: I want to find out is it possible to make the binary GCD algorithm works with multiple-precision integers of any length in terms of "limbs", and if it possible - to get any help or ideas how to implement it in C. Does anyone have any idea or code parts how to implement it?

I'd like to consider any advice or any other solution to solve that problem.

Here is the part of GNU MP binary GCD method for (a = b = 2 limbs) if anyone would take a look:

/* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
   Both U and V must be odd. */
static inline mp_size_t
gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
{
  printf("gcd_2 invoked\n");
  mp_limb_t u0, u1, v0, v1;
  mp_size_t gn;

  u0 = up[0];
  u1 = up[1];
  v0 = vp[0];
  v1 = vp[1];

  ASSERT (u0 & 1);
  ASSERT (v0 & 1);

  /* Check for u0 != v0 needed to ensure that argument to
   * count_trailing_zeros is non-zero. */
  while (u1 != v1 && u0 != v0)
    {
      unsigned long int r;
      if (u1 > v1)
  {
    sub_ddmmss (u1, u0, u1, u0, v1, v0);
    count_trailing_zeros (r, u0);
    u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
    u1 >>= r;
  }
      else  /* u1 < v1.  */
  {
    sub_ddmmss (v1, v0, v1, v0, u1, u0);
    count_trailing_zeros (r, v0);
    v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
    v1 >>= r;
  }
    }

  gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);

  /* If U == V == GCD, done.  Otherwise, compute GCD (V, |U - V|).  */
  if (u1 == v1 && u0 == v0)
    return gn;

  v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
  gp[0] = mpn_gcd_1 (gp, gn, v0);

  return 1;
}

CodePaste of the above.

Answer 1

Just roll your own code specifically for this problem, why not? (10^9)^2 fits into 64-bit int, so you can work with base-(10^9) digits, each held in a 64-bit int. To represent 2^(10^5)-bits values, 2^(10^5) ~= 10^30103, i.e. values with ~ 30103 decimal digits, you will only need 30103/9 ~= 3350 ints which is a ~ 27 kB array in memory, for each of the two numbers involved.

According to WP, for binary GCD algorithm you only need minus and /2 which is trivial to implement by halving each digit, with occasional carry of 5*10^8 into the lower digit (94 / 2 = 47 = {4,5+2}). The final multiply by 2^k can be done with a naive algorithm, since it needs to be done only once.

Printing in base-10 will be trivial. If you don't care for the printing of the final result, then you won't need the final multiply (or if you report your result as 2^k*x) and you could then work with base-10^18 digits, halving the number of digits to work with.

You will only need the usual C integer arithmetic to work with the digits.