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, n1 and n2,
we can express them as k1d and k2d,
where d is the greatest common divisor. From the definition of
greatest common divisor, k1 and k2
are relatively prime.
Without loss of generality we can assume
k1 > k2 > 0. We can then express k1
in terms of k2:
k1 = c1 k2 + k3
such that c1 > 0 and k2 > k3
> 0.
Not only that, k3 is relatively prime to k2
, because d2d (where d2 is
the greatest common divisor of k2 and k3)
would be a common divisor of n1 and n2 greater
than d, contrary to our definition.
We can repeat the above step over and over:
k2 = c2 k3 + k4
k3 = c3 k4 + k5
and so on, getting a cascading series of k's such that
ki>
ki+1 > 0, and where ki and ki+1
are relatively prime (if they weren't, diwould
divide ki-1 which would then not be relatively prime
to ki).
Eventually, we must reach an m such that km
= 1.
Then, km-1 = km-1 km + 0.
Working backwards, we can assign ni = kid
and get a traditional representation of Euclid's algorithm:
n1 = c1 n2 + n3
n2 = c2 n3 + n4
...
nm-2 = cm-2 nm-1 + nm
nm-1 = cm-1 nm + 0.
where d = nm. And of course, if nm
= 1, n1 and n2 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:
ni = dki+1 * ci + dki+2
we can also say
ni = dki+1*(ci+1) - d(ki+1-ki+2).
since ki+1 > ki+1-ki+2
> 0, continuing the algorithm from ki+1-ki+2 will
also provide a correct result, and will save a step whenever ki+1-ki+2
< ki+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