Euclid's algorithm

Given two integers, a and b, returns their greatest common divisor, gcd(a,b). It is perhaps the first known algorithm, dating back well over two millenia! Commonly attributed to Euclid.

if a>b
   swap(a,b)
while a>0 {
   x = b mod a
   b = a
   a = x
}
return b

(Strictly speaking, this is only for natural numbers; you should first take the absolute value of a and b...)

Obviously this algorithm terminates, since at every pass through the loop a is assigned the result of a "mod a" operation, which necessarily decreases it. a is always non-negative, so only a finite number of passes may occur.

The way it works is simple. At the beginning of every pass, the gcd of the values of a and b is the same as the gcd of the original inputs. To show this, we note that b mod a is b minus a multiple of a; therefore, any common divisor of a and b is a common divisor of a and b mod a, and vice versa, so the gcds must be identical. At the last step, it returns b, which is just right! (gcd(0,b)=b, since any natural number divides 0 -- it's not just a quibble: consider the last step, which starts with b a multiple of a)

Note that this algorithm will work in any ring, not just the integers.

Historically (and to some extent mathematically), Euclid's algorithm is really interesting because it is still the best algorithm we have today. In fact, it is used extensively in factorization (see factoring algorithms) and primality testing, two modern, fast paced fields in mathematics. Very impressive, considering its origin is placed around 200 BC.

Unbeknownst to some, a couple extra steps let you use to Euclid's Algorithm to find not only the greatest common denominator g of two integers b and c, but also the values of x and y for which bx + cy = g. We here utilize the fact that the gcd of two (or more) numbers is also their smallest positive linear combination.

A simple demonstration, which you can probably convert into C faster than I can:

qi  ri    xi   yi
    71    1    0
1   50    0    1
2   21    1   -1
2    8   -2    3
1    5    5   -7
1    3   -7   10
1    2   12  -17
     1  -19   27

gcd(71,50) = 1 = (71)(-19) + (50)(27)

Each row requires four swift calculations:

• q_i+1 = Floor(r_i/r_i+1)
• r_i+1 = r_i-1 - q_i*r_i
• x_i+1 = x_i-1 - q_i*x_i
• y_i+1 = y_i-1 - q_i*y_i

Of course, r₀ and r₁ are b and c, in order of decreasing magnitude

Euclid's Algorithm can be difficult to understand if you don't have an example. So I've prepared one. I hope it helps.

Problem: Find the greatest common divisor of the numbers 636 and 483.

Solution: In essence, this means "find the greatest number that fits evenly into both 636 and 483". You take the following steps:

636 = 483 * 1 + 153      /* 153 is the remainder from 636 ÷ 483. */
483 = 153 * 3 + 24
153 =  24 * 6 + 9
24  =   9 * 2 + 6
9   =   6 * 1 + 3
6   =   3 * 2 + 0        /* Stop when remainder = 0 */

So 3 is the GCD of 636 and 483.

Any questions, suggestions or errors, please msg me.

Wonder why Euclid's Algorithm works?

We aleady know that any pair of positive integers has a greatest common divisor; the trick is to find it.

Starting with two distinct numbers, n₁ and n₂, we can express them as k₁d and k₂d, where d is the greatest common divisor. From the definition of greatest common divisor, k₁ and k₂ are relatively prime.

Without loss of generality we can assume k₁ > k₂ > 0. We can then express k₁ in terms of k₂:

k₁ = c₁ k₂ + k₃ such that c₁ > 0 and k₂ > k₃ > 0.

Not only that, k₃ is relatively prime to k₂ , because d₂d (where  d₂ is the greatest common divisor of k₂ and k₃) would be a common divisor of n₁ and n₂ greater than d, contrary to our definition.

We can repeat the above step over and over:

k₂ = c₂ k₃ + k₄
k₃ = c₃ k₄ + k₅

and so on, getting a cascading series of k's  such that k_i> k_i+1 > 0, and where k_i and k_i+1 are relatively prime (if they weren't, d_iwould divide k_i-1 which would then not be relatively prime to k_i).

Eventually, we must reach an m such that  k_m = 1.

Then,  k_m-1 = k_m-1 k_m + 0.  Working backwards, we can assign n_i = k_id and get a traditional representation of Euclid's algorithm:

n₁ = c₁ n₂ + n₃
n₂ = c₂ n₃ + n_4
...
n_m-2 = c_m-2 n_m-1 + n_m
n_m-1 = c_m-1 n_m + 0.

where d = n_m.  And of course, if n_m = 1, n₁ and n₂ are relatively prime.

---

A slight modification of the original algorithm lets you save a couple of steps in some cases. Since the algorithm takes fewer steps for smaller numbers, you should use the "least positive remainder" for each step. That is, if step i is represented:

n_i = dk_i+1 * c_i + dk_i+2

we can also say

n_i = dk_i+1*(c_i+1) - d(k_i+1-k_i+2).

since k_i+1 > k_i+1-k_i+2 > 0, continuing the algorithm from k_i+1-k_i+2 will also provide a correct result, and will save a step whenever k_i+1-k_i+2 < k_i+2


Applying the modification to schmik's example:

636 = 483 * 1 + 153
483 = 153 * 3 +  24
153 =  24 * 6 +   9
24  =   9 * 3 -   3        /* 3 is less than 6 */
9   =   3 * 3 +   0        /* Stop */

So, again, 3 is the GCD of 636 and 483, but we've saved a step. Let's try a more difficult pair of numbers:

973 = 593 * 1 + 380
593 = 380 * 1 + 213
380 = 213 * 1 + 167
213 = 167 * 1 +  46
167 =  46 * 3 +  29
 46 =  29 * 1 +  17
 29 =  17 * 1 +  12
 17 =  12 * 1 +   5
 12 =   5 * 2 +   2
  5 =   2 * 2 +   1
  2 =   2 * 1 +   0  /* unnecessary last step */

showing that 973 and 593 are relatively prime. This becomes:
 

973 = 593 * 2 - 213
593 = 213 * 3 -  46
213 =  46 * 5 -  17
 46 =  17 * 3 -   5
 17 =   5 * 3 +   2
  5 =   2 * 2 +   1
  2 =   2 * 1 +   0

saving 4 steps. Notice that most of the numbers that appear in the second sequence also appeared in the first sequence, but a step earlier.

Finally, another little tidbit: The worst case for Euclid's Algorithm is always on two successive terms in a Fibonacci sequence. Why? The coefficients are always 1, and the remainders take longer to shrink than with other numbers. And they're always relatively prime to boot, damn them!

---

Euclid's Algorithm in C++ at compile time, using recursive templates. Yes, I know it's also done here but a lot of compilers will choke on the simple version. This version also uses the 'least positive remainder' shortcut to save on template specializations.

#include <ostream>

using std::ostream;
using std::endl;

#ifdef STATIC_CONST_WORKS_THE_WAY_IT_SHOULD

#  define MAGIC_NUM(n, value) static const unsigned long n = value

#else

//
// enum hack macro
//
#  define MAGIC_NUM(n, value) enum n ## _ { n = value }

#endif

template<unsigned long i, unsigned long j>
struct euclids
{
 MAGIC_NUM(  lesser, (i<j)?i:j) );
 MAGIC_NUM( greater, (i<j)?:j:i );
 MAGIC_NUM(       d, (lesser>0)?lesser:1) );
 MAGIC_NUM(       r, (lesser>0)?greater%d:0 );
 MAGIC_NUM(      r2, lesser - r );
 MAGIC_NUM(     lpr, (r2<r)?r2:r );

 typedef euclids<lesser,lpr> recur_type;

 MAGIC_NUM(     gcd, (lpr>0)?((lpr>1)?recur_type::gcd:1):lesser );
 MAGIC_NUM(     lcm, greater*(lesser/gcd) );

 static ostream announce (ostream& os)
 {
  os << "euclids <" << i << ", " << j << ">" << endl;
  // os << "   d: " << d   << endl;
  // os << "   m: " << m   << endl;
  // os << "   r: " << r   << endl;
  // os << " lpr: " << lpr << endl;
  // os << " gcd: " << gcd << endl;
  // os << " lcm: " << lcm << endl;

  if (lpr>0) recur_type::announce (os);
  return os;
 } // announce
}; // struct euclids

Euclid's Algorithm has a very elegant recursive definition:

• Base Case: GCD(A, 0) = A
• Recusive Case: GCD(A, B) = GCD(B, A%B)

This of course leads to a one-line algorithm in most any high level programming language. The function is tail recursive, so if you have a good compiler no call stack is built up.