Binary GCD algorithm

From Wikipedia, the free encyclopedia

The binary GCD algorithm is an algorithm which computes the greatest common divisor of two nonnegative integers. It gains a measure of efficiency over the ancient Euclidean algorithm by replacing divisions and multiplications with shifts, which are cheaper when operating on the binary representation used by modern computers. This is particularly critical on embedded platforms that have no direct processor support for division. While the algorithm was first published in modern times by Josef Stein in 1961, it may have been known in first century China (Knuth, 1998).

1 The algorithm
2 Implementation in C
3 Implementation in assembly
4 See also
5 References
6 External links

[edit] The algorithm

The algorithm reduces the problem of finding the GCD by repeatedly applying these identities:

gcd(0, v) = v, because everything divides zero, and v is the largest number that divides v. Similarly, gcd(u, 0) = u. gcd(0, 0) is not defined.
If u and v are both even, then gcd(u, v) = 2·gcd(u/2, v/2), because 2 is a common divisor.
If u is even and v is odd, then gcd(u, v) = gcd(u/2, v), because 2 is not a common divisor. Similarly, if u is odd and v is even, then gcd(u, v) = gcd(u, v/2).
If u and v are both odd, and u ≥ v, then gcd(u, v) = gcd((u-v)/2, v). If both are odd and u < v, then gcd(u, v) = gcd(u, (v-u)/2). These are combinations of one step of the simple Euclidean algorithm, which uses subtraction at each step, and an application of step 3 above. The division by 2 results in an integer because the difference of two odd numbers is even.
Repeat steps 3-4 until u = v, or (one more step) until u = 0. In either case, the result is 2^kv, where k is the number of common factors of 2 found in step 2.

Since this definition is tail-recursive, a loop can be used to replace the recursion.

The algorithm requires O((log₂ uv)²) worst-case time, or in other words time proportional to the square of the number of bits in u and v together. Although each step reduces at least one of the operands by at least a factor of 2, the subtract and shift operations do not take constant time for very large integers (although they're still quite fast in practice, requiring about one operation per word of the representation).

Although this algorithm can outperform the traditional Euclidean algorithm on many platforms, its asymptotic performance is the same, it is considerably more complex, and it does not generalize to negative integers or other rings such as polynomials where division with remainder is defined. An extended version, analogous to the extended Euclidean algorithm, is given in (Knuth 1998, answer to exercise 39 of section 4.5.2, p. 646) along with pointers to other versions.

[edit] Implementation in C

Following is an implementation of the algorithm in the C programming language, taking two non-negative arguments u and v. It first removes all common factors of 2 using identity 2, computes the gcd of the remaining numbers using identities 3 and 4, then combines these to form the final answer.

The version on the left is more readable, the version on the right is heavily optimized. A compiler should be able to transform one in the other, but in practice the results will be considerably different depending on the optimization algorithms employed by the compiler.

int gcd (int a, int b)
{
  int shift = 0, diff;
  while ((a & 1) == 0 && (b & 1) == 0)
    a >>= 1, b >>= 1, shift++;

  while ((a & 1) == 0)
    a >>= 1;

  do
    {
      while ((b & 1) == 0)
        b >>= 1;

      diff = a - b;
      if (diff < 0)
        b = -diff;
      else
        a = b, b = diff;

      b >>= 1;
    }
  while (diff);
  return a << shift;
}

int gcd (int a, int b)
{
  int shift = 0, diff;
  while ((a & 1) == 0)
    {
      a >>= 1;
      if ((b & 1) == 0)
        b >>= 1, shift++;
    }

  while ((b & 1) == 0)
    {
     loop:
      b >>= 1;
    }
  diff = a - b;
  if (diff < 0)
    b = -diff;
  else
    a = b, b = diff;
  if (diff)
    goto loop;

  return a << shift;
}

Note that division-by-2 and modulo-2 operations are implemented as bitwise operations. Any modern C compiler can do that for unsigned operands.

Efficiency can be enhanced on platforms with special instructions to find the least significant 1 bit in a word, allowing trailing zero bits to be removed in much larger groups.

[edit] Implementation in assembly

This version of binary GCD in ARM assembly (using GNU Assembler syntax) highlights the benefit of branch predication. The loop to implement step 2 consists of 6 instructions, four of which are predicated. Step 3 consists of two loops, each 3 instructions long (one of the instructions being predicated); however, after the first iteration r0 is kept odd and need not be tested, and only one of the loops is executed. Finally, step 4 takes four instructions of which two are predicated.

gcd:
   @ Arguments arrive in registers r0 and r1
   mov     r3, #0            @ Initialize r3, the power of 2 in the result

remove_twos_loop:
   tst     r0, #1            @ Check least significant bit of r0.  If zero...
   tsteq   r1, #1            @ ... check least significant bit of r1 too.  If zero:
   addeq   r3, r3, #1        @ Increment r3, marking one more common 2 factor;
   moveq   r0, r0, lsr #1    @ Shift r0 right, removing the common 2 factor;
   moveq   r1, r1, lsr #1    @ Shift r1 right, removing the common 2 factor; ...
   beq     remove_twos_loop  @ ... And Check if we can remove one more

check_two_r0:
   tst     r0, #1            @ Check least significant bit of r0
   moveq   r0, r0, lsr #1    @ If r0 is even, shift r0 right, dividing it by 2, and...
   beq     check_two_r0      @ ... Check if we can remove one more

check_two_r1:
   tst     r1, #1            @ Check least significant bit of r1
   moveq   r1, r1, lsr #1    @ If r1 is even, shift it right by 1 and...
   beq     check_two_r1      @ ... Check if we can remove one more

   subs    r2, r0, r1        @ Set r2 = r0 - r1 and set flags
   rsblt   r2, r2, #0        @ Set r2 = |r0 - r1| (actually negate r2 if r0 < r1)
   movgt   r0, r1            @ Set r0 = min (r0, r1) (actually move r1 to r0 if r0 > r1)
   mov     r1, r2, lsr #1    @ Set r1 = |r0 - r1| / 2 (useless on last iteration)
   bne     check_two_r1      @ Uses flags from subs; was r0 = r1? If not, iterate
                             @ r0 is one of the old r0 and r1, so it's odd: check only r1

   mov     r0, r0, asl r3    @ Put r0 << r3 in the result register
   bx      lr                @ Return to caller