Greatest Common Divisor Algorithms, Ancient and Modern

by Michael Gilleland

Introduction

The easiest way to compute the greatest common divisor (or highest common factor, as it is sometimes called) of two non-negative integers in Java is to use the built-in gcd method of class java.math.BigInteger. But this is akin to looking up the answers to the daily crossword puzzle. It smacks of cheating.

There is a Latin expression "homines dum docent discunt," which means "people learn while they teach" (Seneca, Epistulae Morales, 7.8). But if you have no students to teach, another way to learn is to program. The mathematician Doron Zeilberger says, "Whenever I want to learn a new subject, I announce a graduate course in it, since the best way to learn is by teaching. But even better than teaching humans is teaching computers, i.e. program!"

I wrote the Java class below to teach myself more about greatest common divisor algorithms. It shows two ways to implement Euclid's algorithm, one iterative and the other recursive. It also contains an implementation of a modern binary algorithm, invented by Josef Stein in 1961, which finds the greatest common divisor without dividing. Both algorithms are discussed in Donald E. Knuth, The Art of Computer Programming, 3rd edition, volume 2 = Seminumerical Algorithms (Addison-Wesley, 1998), section 4.5.2 (The Greatest Common Divisor), pp. 333-353. The "main" test routine plays with a theorem by Dirichlet: If u and v are integers chosen at random, the probability that gcd (u, v) = 1 (i.e. that u and v are relatively prime) is 6 / (pi squared), or approximately .60793.

Source Code

//------------------------------------------------------
// This class uses different algorithms to compute the
// greatest common divisor of two non-negative integers.
//------------------------------------------------------

import java.math.BigInteger;
import java.util.Random;

public abstract class GreatestCommonDivisor {

  //------------------------------------
  // Euclid algorithm, Elements, Book 7,
  // Propositions 1-2 (circa 300 BC)
  //------------------------------------

  public static long euclidAlgorithmIterative (long u, long v) {

    trace ("Euclid algorithm,  iterative");
    long gcd = 0;
    long r = 0;

    u = Math.abs (u);
    v = Math.abs (v);

    while (true) {
      if (v == 0) {
        gcd = u;
        break;
      }
      else {
        r = u % v;
        trace ("u " + u + ", v " + v + ", u % v " + r);
        u = v;
        v = r;
      }
    }

    trace ("GCD " + gcd);
    return gcd;

  }

  public static long euclidAlgorithmRecursive (long u, long v) {

    trace ("Euclid algorithm,  recursive: u  " + u + ", v " + v);
    long gcd = 0;
    long r = 0;

    u = Math.abs (u);
    v = Math.abs (v);

    if (v == 0) {
      gcd = u;
      trace ("GCD " + gcd);
    }
    else {
      gcd = euclidAlgorithmRecursive (v, u % v);
    }

    return gcd;

  }

  //--------------------------------------------------------
  // Is number even or odd? We have to do this without
  // resorting to division, no matter how inefficient it is.
  //--------------------------------------------------------

  private static boolean isEven (long n) {
    String s = Long.toBinaryString (n);
    int length = s.length ();
    String last = s.substring (length - 1, length);
    return last.equals ("0");
  }

  private static boolean isOdd (long n) {
    return !isEven (n);
  }

  //-------------------------------------
  // Binary algorithm (uses no division),
  // invented by Josef Stein in 1961.
  // See J. Comp. Phys. 1 (1967) 397-405
  //-------------------------------------

  private static long binaryAlgorithm (long u, long v) {

    trace ("Binary algorithm, u " + u + ", v " + v);
    long gcd = 0;
    long k = 0;
    long t = 0;

    u = Math.abs (u);
    v = Math.abs (v);

    // Find power of 2

    while (true) {

      if (isOdd (u)) {
        break;
      }

      if (isOdd (v)) {
        break;
      }

      k++;
      u = u >> 1;
      v = v >> 1;
      
    }

    trace ("k " + k);

    // Initialize

    if (isOdd (u)) {
      t = -v;
    }
    else {
      t = u >> 1;
    }

    while (true) {

      trace ("u " + u + ", v " + v + ", t " + t);

      // Make t even

      if (isEven (t)) {
        t = t >> 1;
      }

      else {

        // Reset max (u, v)

        if (t > 0) {
          u = t;
        }
        else {
          v = -t;
        }

        // Subtract and check for completion

        t = u - v;
        if (t == 0) {
          gcd = u * (long) Math.pow (2, k);
          break;
        }
      }
    }

    trace ("GCD " + gcd);
    return gcd;

  }

  //---------------------------
  // Use Java's built in method
  //---------------------------

  public static long builtInMethod (long u, long v, boolean traceFlag) {
    BigInteger bigU = new BigInteger (Long.toString (u));
    BigInteger bigV = new BigInteger (Long.toString (v));
    long gcd = bigU.gcd (bigV).longValue ();
    if (traceFlag) {
      trace ("GCD (" + u + ", " + v + ") = " + gcd);
    }
    return gcd;
  }

  //------
  // Trace
  //------

  private static void trace (String s) {
    System.out.println (s);
  }

  //-----
  // Test
  //-----

  public static void main (String[] args) {

    // These numbers are from Knuth

    long u = 40902;
    long v = 24140;
    long gcd = 0;

    gcd = GreatestCommonDivisor.euclidAlgorithmIterative (u, v);
    gcd = GreatestCommonDivisor.euclidAlgorithmRecursive (u, v);
    gcd = GreatestCommonDivisor.binaryAlgorithm (u, v);
    gcd = GreatestCommonDivisor.builtInMethod (u, v, true);

    // Dirichlet's theorem

    double total = 1000;
    double relativelyPrime = 0;
    Random generator = new Random (System.currentTimeMillis ());
    for (int i = 0; i < (int) total; i++) {
      u = generator.nextLong ();
      v = generator.nextLong ();
      gcd = GreatestCommonDivisor.builtInMethod (u, v, false);
      if (gcd == 1) {
        relativelyPrime++;
      }
    }
    System.out.println ("Dirichlet's theorem, actual " +
                        relativelyPrime / total +
                        ", expected 0.608");

  }

}