Editing Greatest common divisor

The '''greatest common divisor (GCD)''' of a set of integers is the greatest integer that divides all the integers in the set, unless the set contains only zeroes, in which the GCD is defined to be zero. For example, the GCD of 12 and 20 is 4, because 4 divides both 12 and 20, and no integer larger than 4 divides both 12 and 20.

==Properties==
* <math>\gcd(x,y) = \gcd(y,x)</math> (commutativity)
* <math>\gcd(0,x) = |x|</math>
* <math>\gcd(1,x) = 1</math>
* <math>\gcd(x,x) = |x|</math>
* All common divisors of a set of integers divide the GCD of that set. For example, the common divisors of 24 and 42 are -6, -3, -2, -1, 1, 2, 3, and 6, and the GCD is 6. Notice that all the common divisors are also divisors of 6.
* The GCD is never negative, because if <math>-x</math> is a common divisor, then so is <math>+x</math>. Also, <math>\gcd(x,y) = \gcd(|x|,|y|)</math>.
* If <math>S</math> and <math>T</math> are sets, then <math>\gcd(S \cup T) = \gcd(\gcd(S), \gcd(T))</math>. For example, <math>\gcd(x,y,z) = \gcd(x,\gcd(y,z)) = \gcd(\gcd(x,y),z)</math> (associativity).
* <math>\gcd(x,y)\,\operatorname{lcm}(x,y) = |xy|</math>.
* <math>\gcd(x,\operatorname{lcm}(x,y)) = \operatorname{lcm}(x,\gcd(x,y))</math>. (This property is the absorption law of the lattice over the positive integers with GCD as meet and LCM as join.)
* If <math>x</math> has prime factorization <math>(-1)^a p_1^{k_1} p_2^{k_2} ... p_m^{k_m}</math> and <math>y</math> has prime factorization <math>(-1)^b p_1^{l_1} p_2^{l_2} ... p_m^{l_m}</math>, then <math>\gcd(x,y) = p_1^{\min(k_1,l_1)} p_2^{\min(k_2,l_2)} ... p_n^{\min(k_n,l_n)}</math>. Therefore, it is possible to find the GCD of two integers by factorizing them; but much more efficient methods exist.
* <math>\gcd(cx, cy) = |c|\gcd(x,y)</math> whenever <math>c \in \mathbb{Z}</math>.
* <math>\gcd(x, y) = \gcd(x+ky, y)</math> for any <math>k \in \mathbb{Z}</math>.
* There exist integers <math>a, b</math> such that <math>ax + by = 1</math> if and only if <math>\gcd(x,y) = 1</math>. (This is called ''Bézout's identity''. We will not prove it directly, but we will exhibit an algorithm for determining <math>a</math> and <math>b</math>.)

==Euclidean algorithm==
Means of computing the GCD of two integers efficiently have been known since antiquity. The well-known ''Euclidean algorithm'' uses the property that <math>\gcd(x,y) = \gcd(x-y,y)</math> for any <math>x, y</math>. When given two positive integers, we repeatedly subtract the smaller from the larger, thus reducing the size of the larger integer; and we repeat this until the two integers become identical. For example, <math>\gcd(24, 42) = \gcd(24, 18) = \gcd(6, 18) = \gcd(6, 12) = \gcd(6, 6) = 6</math>. The following is a naive implementation:
<syntaxhighlight lang="c">
int gcd(int x, int y)
{
    if (x < 0) x = -x;
    if (y < 0) y = -y;
    if (x == 0) return x;
    if (y == 0) return y;
    if (x == y) return x;
    if (x > y) return gcd(x - y, x);
    else return gcd(x, y - x);
}
</syntaxhighlight>
In practice, negative numbers are rarely encountered. Also, when we obtain two identical integers, if we perform one additional subtraction, we will get zero, which we have to handle in any case, so:
<syntaxhighlight lang="c">
int gcd(int x, int y)
{
    if (x == 0 || y == 0) return x + y;
    if (x >= y) return gcd(x - y, y);
    else return gcd(x, y - x);
}
</syntaxhighlight>
To speed this up, we note that if we start out with <math>x \gg y</math>, then the final result of subtracting <math>y</math> from <math>x</math> over and over again will be the least non-negative residue of <math>x</math> ''modulo'' <math>y</math>. Hence, instead of performing repeated subtraction, we use the mod operator. Finally, with some cleverness, we obtain this common form:
<syntaxhighlight lang="c">
int gcd(int x, int y)
{
    if (x == 0) return y;
    else return gcd(y%x, x);
}
</syntaxhighlight>
Also, although we have written this algorithm [[Recursive function|recursively]], we can also write it iteratively, since it is tail recursive:
<syntaxhighlight lang="c">
while (x > 0)
{
    y %= x;
    swap(x, y);
}
/* y now contains the GCD of original (x, y) */
</syntaxhighlight>
Now, regarding performance:

''Theorem'': The GCD of two positive integers <math>a</math> and <math>b</math> can be computed using the Euclidean algorithm in <math>O(\log\max(a,b))</math> mod operations.

''Sketch Proof'': The idea is to use induction on <math>\max(a,b)</math>. Without loss of generality, we assume <math>a \geq b</math>. Then, there are two possibilities. If <math>b \leq \frac{a}{2}</math>, then <math>\max(a,b)</math> will be cut at least in half after the first mod operation. But if <math>b > \frac{a}{2}</math>, then <math>b</math> goes into <math>a</math> only once, and we obtain <math>a-b < \frac{a}{2}</math> and <math>b</math> after the first mod operation; a second will then take <math>b</math> below <math>\frac{a}{2}</math>. So in either case, after at most two mod operations, <math>\max(a,b)</math> is reduced by at least a factor of 2.

''Note'': A more careful analysis shows that the worst-case number of mod operations required to find <math>\gcd(a,b)</math> using the Euclidean algorithm goes approximately as <math>\log_{\varphi}(\max(a,b))</math>, where <math>\varphi = \frac{1+\sqrt{5}}{2}</math>.

The bit complexity of the Euclidean algorithm is <math>O(uv)</math>, so that finding the GCD of two large numbers with <math>u</math> and <math>v</math> bits, respectively, requires <math>O(uv)</math> time.<ref>Jeffrey Shallit. ''The Greatest Common Divisor'', pp. 22-22. CMI Introductory Workshop, Algorithmic Number Theory. Retrieved from www.cs.uwaterloo.ca/~shallit/Talks/gcd.ps</ref> For two numbers of approximately equal length <math>n</math>, this is <math>O(n^2)</math>. By comparison, the straightforward algorithm based on factorizing the input integers does not run in polynomial time, and hence will perform far worse for large <math>n</math>.

===Extended Euclidean algorithm===
The ''extended Euclidean algorithm'' is the name given to the well-known procedure used to find integer solutions to the equation <math>ax + by = 1</math>, where <math>a</math> and <math>b</math> are given with <math>\gcd(a,b) = 1</math>. For example, we can use it to solve the equation <math>13x + 17y = 1</math> for <math>x</math> and <math>y</math> (one possible solution is <math>x = 4, y = -3</math>.

This equation only has solutions when <math>\gcd(a,b) = 1</math>, because if <math>\gcd(a,b) = d > 1</math>, then <math>ax + by</math> will always be divisible by <math>d</math>, and hence can never be 1. However, in this case <math>\gcd(a/d, b/d) = 1</math>, so we can find solutions to the equation <math>(a/d)x + (b/d)y = 1</math>, that is, <math>ax + by = d</math>. So in general we can always solve <math>ax + by = \gcd(a,b)</math> in integers.

The reason why this is known as the extended Euclidean algorithm is that it is based on reducing <math>a</math> and <math>b</math> as in the Euclidean algorithm, recursively, until they become small enough so that the equation is easy to solve; then reconstructing the solution to the original equation using the results.

In particular, when we take the remainder of <math>b</math> modulo <math>a</math> and obtain <math>r</math>, this tells us that <math>b = ka + r</math> for some integer <math>k</math>. So each step of the Euclidean algorithm expresses the fact that <math>\gcd(a, ka+r) = \gcd(a, r) = \gcd(r,a)</math>.

Now suppose we can find integers <math>x', y'</math> such that <math>x'r + y'a = d</math>. Then, recall that <math>r = b - ka</math>; therefore, <math>x'(b - ka) + y'a = d</math>. Expanding, we obtain <math>(y' - x'k)a + x'b = d</math>. So set <math>x = y' - x'k, y = x'</math>, and we now have a solution to <math>ax + by = d</math>.

This suggests that, in general, to solve <math>ax + by = d</math>, we first compute quotient <math>k</math> and remainder <math>r</math> when <math>b</math> is divided by <math>a</math>; and then (recursively) find some solution to <math>rx' + ay' = d</math>; and finally obtain <math>x = y' - kx', y = x'</math>.

This gives the following implementation, where we assume that <math>a</math> and <math>b</math> are positive:
<syntaxhighlight lang="c">
/* Returns the GCD. */
int gcd_extended(int a, int b, int* x, int* y)
{
    if (a == 0)
    {
        *x = 0, *y = 1; /* 0(0) + 1(b) = b = gcd(a,b) */
        return b;
    }
    else
    {
        int k = b/a;
        int r = b%a;
        int _x, _y;
        int d = gcd_extended(r, a, &_x, &_y);
        *x = _y - k*_x;
        *y = _x;
        return d;
    }
}
</syntaxhighlight>
One specific application of the extended Euclidean algorithm is that of computing modular inverses. That is, given some integer <math>p</math> and some modulus <math>m</math>, we wish to find another integer <math>q</math> such that <math>pq \equiv 1 \pmod{m}</math>. For example, if <math>p = 4</math> and <math>m = 21</math>, then <math>q = -5</math> works, since <math>pq = -20 \equiv 1 \pmod{21}</math>.

All we need to do to solve this is to note that the statement that <math>pq \equiv 1 \pmod{m}</math> is equivalent to <math>pq = km + 1</math> for some integer <math>k</math>, or <math>pq + (-m)k = 1</math>, to be solved for integers <math>q, k</math> (even though we might not care what <math>k</math> is). We then use the extended Euclidean algorithm. Note that this suggests that a solution exists only when <math>\gcd(p,m) = 1</math>, which is indeed the case.

==Binary GCD algorithm==
This is another popular GCD algorithm, also recursive. 
* If <math>x</math> and <math>y</math> are both even, then <math>\gcd(x,y) = 2\gcd\left(\frac{x}{2},\frac{y}{2}\right)</math>.
* If <math>x</math> is even, but <math>y</math> is odd, then <math>\gcd(x,y) = \gcd\left(\frac{x}{2},y\right)</math>. This is because every common divisor of <math>x</math> and <math>y</math> is odd, and hence will also divide <math>\frac{x}{2}</math>.
* If <math>x</math> and <math>y</math> are both odd, then <math>\gcd(x,y) = \gcd(x-y, y)</math>, as always.

This suggests the following naive implementation, assuming that <math>x</math> and <math>y</math> are both nonnegative:
<syntaxhighlight lang="c">
int gcd(int x, int y)
{
    if (x == 0) return y;
    if (x < y) return gcd(y, x);
    if (x%2 == 1)
        if (y%2 == 1)
            return gcd(x - y, y);
        else
            return gcd(x, y/2);
    else
        if (y%2 == 1)
            return gcd(x/2, y);
        else
            return gcd(x/2, y/2)*2;
}
</syntaxhighlight>
Computing the GCD of <math>a</math> and <math>b</math> using this algorithm will take <math>O(\log\max(a,b))</math> steps, just as the Euclidean algorithm. This becomes clear when we write out the binary representations of <math>a</math> and <math>b</math>:
* If either <math>a</math> or <math>b</math> is even, then it will be divided by two, which reduces its length by one.
* If both <math>a</math> and <math>b</math> are odd, then one of them will become even after a subtraction, and then it will be divided by two, which again reduces its length by one. In effect <math>\gcd(a,b) = \gcd((a-b)/2, b)</math> (for the case <math>a = b</math>).
Therefore, we obtain a bound of <math>2(\log_2 a + \log_2 b) = 2 \log_2 \max(a,b)</math> iterations, since after every two iterations it is guaranteed that the binary length of at least one of the two numbers is reduced.

Furthermore, each operation (either halving an even integer, or performing a subtraction) has bit complexity proportional to the length of the integers involved. If we halve <math>a</math>, then <math>O(\log a)</math> time is required for this step; if we subtract <math>b</math> from <math>a</math>, then <math>O(\log a)</math> time is required for this step, and in the next step <math>a</math> will be halved, also requiring <math>O(\log a)</math> time. So at most <math>O(\log a)</math> operations of <math>O(\log a)</math> bit complexity each can be performed on <math>a</math>, and at most <math>O(\log b)</math> operations of <math>O(\log b)</math> bit complexity each; so the overall bit complexity is <math>O(u^2 + v^2)</math>, where <math>u = \log a</math> and <math>v = \log b</math>. This will often be expressed as <math>O((u+v)^2)</math> (quadratic in total input size) or <math>O(\max(u,v)^2)</math>, matching the bound on the Euclidean algorithm.

The binary GCD algorithm can be considered to be motivated by the fact that reduction modulo 2, division by 2, and subtraction admit simpler bit complexities than the general modulo operation that the Euclidean algorithm uses, which suggests that it should be faster in practice. However, the naive implementation shown above is unlikely to outperform the optimized Euclidean algorithm in the previous section; optimization is required in order for the binary GCD algorithm to be competitive. 

==References==
<references/>