Euclidean Algorithm, in depth: why gcd(a, b) = gcd(b, a mod b), a traced example, safe code for negatives, binary GCD and real uses

The idea and why it is correct

The obvious method is to try every candidate from the smaller number down to 1. For numbers near a billion that is up to a billion divisions, and for cryptographic sizes it never finishes. Factoring both numbers and intersecting the prime factors is no better, because factoring large integers is itself hard. Euclid's insight avoids factoring entirely.

The key fact: for integers a and b with b not zero, gcd(a, b) = gcd(b, a mod b).

Proof. Write a = q b + r, where r = a mod b. If d divides a and b, then d divides a - q b, which is r, so every common divisor of a and b is a common divisor of b and r. Conversely, if d divides b and r, it divides q b + r, which is a. The two pairs therefore have exactly the same set of common divisors, so they have the same greatest one.

Termination. The remainder satisfies 0 <= r < b, so the second argument strictly decreases while staying non-negative. It must reach 0. When it does, gcd(g, 0) = g, because every integer divides 0 and the largest divisor of g is g itself. That is the whole algorithm: replace (a, b) by (b, a mod b) until b is 0, then return a.

Euclid's original form subtracted the smaller number from the larger instead of dividing. That is the same idea (subtracting b q times is what the division does in one step) but it can be catastrophically slow: gcd(1,000,000,000, 1) takes a billion subtractions and one division.

A trace by hand

Take a = 1071 and b = 462.

The answer is 21. You can check it: 1071 = 3 x 3 x 7 x 17 and 462 = 2 x 3 x 7 x 11, so the common part is 3 x 7 = 21. The algorithm reached it in three divisions without knowing any factor. If the first argument is the smaller one, the first step simply swaps them: gcd(462, 1071) begins with 462 mod 1071 = 462, giving the pair (1071, 462).

Code that handles zero and negative numbers

Python's integers are unbounded and its % operator returns a result with the sign of the divisor, so a short loop works, provided the sign of the answer is fixed at the end:

def gcd(a: int, b: int) -> int:
    """Greatest common divisor; gcd(0, 0) == 0 by convention, result is never negative."""
    while b:
        a, b = b, a % b
    return abs(a)

assert gcd(1071, 462) == 21
assert gcd(-12, 18) == 6 and gcd(12, -18) == 6
assert gcd(0, 5) == 5 and gcd(0, 0) == 0

Without the abs, gcd(12, -18) returns -6, because the remainders take the sign of the negative divisor. The convention that gcd(0, 0) = 0 is not arbitrary: it keeps identities such as gcd(a, b) x lcm(a, b) = |a b| true and is what Python, C++ and most libraries return.

In C and C++ the % operator truncates toward zero, so the remainder has the sign of the dividend. The safest approach is to convert to unsigned magnitudes first. That also handles the most negative value, whose absolute value does not fit in the signed type, so calling abs(INT64_MIN) is undefined behaviour:

#include <stdint.h>

static uint64_t mag(int64_t x) {           /* |x| without signed overflow */
    return x < 0 ? (uint64_t)0 - (uint64_t)x : (uint64_t)x;
}

uint64_t gcd_i64(int64_t a, int64_t b) {
    uint64_t x = mag(a), y = mag(b);
    while (y != 0) {
        uint64_t r = x % y;
        x = y;
        y = r;
    }
    return x;                               /* gcd(INT64_MIN, 0) = 2^63 fits */
}

The result is returned as unsigned on purpose: gcd(INT64_MIN, 0) is 2 to the 63rd power, which no signed 64-bit integer can hold. A recursive version is equally correct and the depth is small (under a hundred for 64-bit inputs), but the loop avoids any question of stack use.

How many steps it takes

Each step at least halves something. If r = a mod b, then r < a / 2: either b <= a / 2 and r < b, or b > a / 2 and r = a - b < a / 2. So every two steps the larger number at least halves, giving at most about 2 log2(min(a, b)) steps. The exact worst case is sharper. Gabriel Lame proved in 1844 that the number of division steps never exceeds five times the number of decimal digits of the smaller number, and the inputs that come closest are consecutive Fibonacci numbers, where every quotient is 1 and the remainders walk down the Fibonacci sequence. gcd(89, 55) takes nine divisions: 89, 55, 34, 21, 13, 8, 5, 3, 2, 1. Random inputs are much kinder; the average number of steps is close to 0.84 times the natural log of the smaller number.

For fixed-size machine integers this means a 64-bit gcd takes at most about 92 divisions and usually far fewer. For big integers the count of steps stays logarithmic, but each step costs a division of multi-word numbers. The bounds on intermediate values and step counts are analysed in more detail on the extended algorithm page.

Binary GCD and big integers

Integer division is one of the slowest instructions on a CPU: a 64-bit divide costs tens of cycles on many cores, while shifts and subtractions cost one. Stein's binary GCD algorithm (1967) avoids division by using three facts: gcd(2a, 2b) = 2 gcd(a, b); gcd(2a, b) = gcd(a, b) when b is odd; and gcd(a, b) = gcd(a - b, b), where a - b is even when both are odd.

#include <stdint.h>

uint64_t gcd_binary(uint64_t a, uint64_t b) {
    if (a == 0) return b;
    if (b == 0) return a;
    int shift = __builtin_ctzll(a | b);     /* common factors of two */
    a >>= __builtin_ctzll(a);               /* make a odd */
    do {
        b >>= __builtin_ctzll(b);           /* make b odd */
        if (a > b) { uint64_t t = a; a = b; b = t; }
        b -= a;                             /* odd - odd = even, or zero */
    } while (b != 0);
    return a << shift;
}

__builtin_ctzll counts trailing zero bits in one instruction on GCC and Clang (MSVC has _BitScanForward64 and C++20 has std::countr_zero). Whether this beats the division loop depends on the CPU: on cores with fast hardware dividers the gap is small, so measure on your target before switching. The Python equivalent of the shift count is (x & -x).bit_length() - 1, but in Python the built-in is faster than either hand-written version.

For numbers of thousands of digits, libraries use other methods. Lehmer's algorithm runs the Euclidean steps on the leading machine word of each number to predict several quotients at once, then applies them to the full numbers in one combined update, replacing many multi-word divisions with a few multiplies. GMP goes further and uses a subquadratic half-gcd method for very large operands.

Library functions

Use the library function when it exists. It handles signs and zero correctly and, for big integers, uses faster methods than any short hand-written loop.

Where gcd is used

Reducing fractions. Divide numerator and denominator by their gcd, and make the denominator positive. Python's fractions.Fraction does exactly this after every operation, which is why fraction arithmetic stays exact.

Least common multiple without overflow. lcm(a, b) = |a b| / gcd(a, b), but computing a b first can overflow even when the lcm fits. Divide first:

from math import gcd

def lcm(a: int, b: int) -> int:
    if a == 0 or b == 0:
        return 0
    return abs(a // gcd(a, b) * b)

def gcd_all(xs):
    g = 0
    for x in xs:
        g = gcd(g, x)
        if g == 1:          # cannot get smaller; stop early
            break
    return g

# gcd of differences: the largest m such that all values are equal mod m
vals = [17, 41, 89, 113]
m = gcd_all([v - vals[0] for v in vals[1:]])   # 24, 72, 96 -> 24
assert all((v - vals[0]) % m == 0 for v in vals)

Starting a fold with 0 works because gcd(0, x) = |x|, and the early exit at 1 makes the gcd of a long array cheap in practice. The gcd-of-differences trick answers questions such as which period explains a set of timestamps, or what modulus makes all values congruent.

Array rotation. Rotating an array of length n by k positions in place splits the indices into exactly gcd(n, k) cycles, each of length n / gcd(n, k). The juggling rotation algorithm walks each cycle once, moving every element exactly once.

Number theory and cryptography. Two numbers are coprime when their gcd is 1, which is the condition for a modular inverse to exist. Key generation in RSA checks that the public exponent is coprime to the totient, solving linear Diophantine equations starts by testing whether the gcd divides the constant, and combining congruences with the Chinese remainder theorem needs gcd of the moduli.

Common bugs

What to do next

1. Trace gcd(252, 105) by hand, then confirm the answer by factoring both numbers.
2. Implement the division loop in your main language and test it on zero, negative inputs and the most negative integer.
3. Replace any a * b / gcd(a, b) in your code with a / gcd(a, b) * b.
4. Time the division loop against the binary version on your target CPU before choosing one for hot code.
5. Count steps for consecutive Fibonacci pairs and check them against Lame's bound.
6. Move on to the extended algorithm and modular inverses, then to modular exponentiation to see how they combine in RSA.