The Euclidean Algorithm
The Euclidean algorithm is a simple but useful procedure to find the greatest common divisor (GCD) of two positive integers. (There are variations for polynomials and other more general numbers.)
This package uses the Euclidean algorithm on numbers of an arbitrary unsigned integer type called UITYPE. You should add an appropriate "typedef" statement to define this type. If your application uses two or more different unsigned integer types, you might consider using templates instead.
The package accesses numbers of type UITYPE by reference because in some applications (such as cryptography) they may be quite large.
The package defines two functions. The first one merely computes the GCD of two positive integers:
GreatestCommonDivisor(x, y, g); const UITYPE &x; positive integer const UITYPE &y; positive integer UITYPE &g; filled in with the GCD of x and y
The second one is a little more sophisticated:
EuclideanAlgorithm(x, y, a, b, g); const UITYPE &x; positive integer const UITYPE &y; positive integer UITYPE &a; see below UITYPE &b; see below UITYPE &g; filled in with the GCD of x and y
This function computes the GCD g of x and y and also two integers a and b such that:
ax - by = g, 1 <= a <= y, 0 <= b < x.
These functions will fail in undefined ways if either x or y is zero, or if there is a numeric overflow. However, if the UITYPE type can hold x and y, it can hold all intermediate values and results without overflow.
The Euclidean algorithm is recusive. If x <= y, we divide y by x to obtain a quotient and remainder with the properties
y = qx + r, 0 <= r < x.
If r = 0, the GCD of x and y is x, so a = 1, b = 0, and g = x.
If r is not 0, then we notice first that every divisor of x and y is also a divisor of r, and that every divisor of x and r is also a divisor of y. Hence the GCD of x and r is the same as the GCD of x and y. We use the algorithm on x and r to obtain:
a'x - b'r = g, 1 <= a' <= r, 0 <= b' < x.
Then it is easy to show that the following values of a and b meet all the requirements:
a = a' + b'q, b = b'.
If x > y, we can use the algorithm with x and y reversed to obtain b'y - a'x = g, 1 <= b' <= x, 0 <= a' < y.
(-a')x - (-b')y = g,
but this isn't quite what we want, because -a' and -b' are out of range in nearly all cases. We can fix this by adding the identity yx - xy = 0 to obtain
(y-a')x - (x-b')y = g,
in which the coefficients are in the desired ranges.
The coefficients are unique within their specified ranges if g = 1. If g > 1, they are not unique, but application of the same algorithm to x/g and y/g yields a and b such that
a(x/g) - b(y/g) = 1, 1 <= a <= (y/g), 0 <= b <= (x/g)-1,
ax - by = g, 1 <= a <= (y/g), 0 <= b <= (x/g)-1,
and these coefficients are unique within these narrower ranges.