title: GCD algorithms for arbitrary-precision arithmetic
body: <p>I am completely stuck with this question, so I am looking for any help.</p>
  
  <p>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.</p>
  
  <p>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".</p>
  
  <p>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.</p>
  
  <p><strong>Question:</strong> 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? </p>
  
  <p>I'd like to consider any advice or any other solution to solve that problem.</p>
  
  <p>Here is the part of GNU MP binary GCD method for (a = b = 2 limbs) if anyone would take a look:</p>
  
  <pre><code>/* Use binary algorithm to compute G &lt;-- 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 &amp; 1);
    ASSERT (v0 &amp; 1);
  
    /* Check for u0 != v0 needed to ensure that argument to
     * count_trailing_zeros is non-zero. */
    while (u1 != v1 &amp;&amp; u0 != v0)
      {
        unsigned long int r;
        if (u1 &gt; v1)
    {
      sub_ddmmss (u1, u0, u1, u0, v1, v0);
      count_trailing_zeros (r, u0);
      u0 = ((u1 &lt;&lt; (GMP_NUMB_BITS - r)) &amp; GMP_NUMB_MASK) | (u0 &gt;&gt; r);
      u1 &gt;&gt;= r;
    }
        else  /* u1 &lt; v1.  */
    {
      sub_ddmmss (v1, v0, v1, v0, u1, u0);
      count_trailing_zeros (r, v0);
      v0 = ((v1 &lt;&lt; (GMP_NUMB_BITS - r)) &amp; GMP_NUMB_MASK) | (v0 &gt;&gt; r);
      v1 &gt;&gt;= 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 &amp;&amp; u0 == v0)
      return gn;
  
    v0 = (u0 == v0) ? ((u1 &gt; v1) ? u1-v1 : v1-u1) : ((u0 &gt; v0) ? u0-v0 : v0-u0);
    gp[0] = mpn_gcd_1 (gp, gn, v0);
  
    return 1;
  }
  </code></pre>
  
  <p><a href="http://codepaste.ru/13601/" rel="nofollow">CodePaste of the above</a>.</p>