The Euclidean Algorithm in Python: The Oldest Algorithm Still in Daily Use

Euclid wrote this one down around 300 BC, in Book VII of the Elements. It is still the code that runs when Python reduces a fraction, and still the code that runs when a machine generates an RSA key. No other algorithm in this series is that old and that current at the same time.

The problem it solves is small: given two whole numbers, find the largest number that divides both of them exactly. The greatest common divisor of 48 and 18 is 6. The obvious approach is to try every candidate from the smaller number downwards until one divides both, which for two eighteen-digit numbers means up to a quintillion tests. Euclid's method answers the same question in at most 90 divisions, and usually far fewer.

That gap — from a quintillion to under a hundred — comes from a single observation about divisors that takes two lines to prove. Once you have it, least common multiples, modular inverses and RSA private keys all fall out of the same loop.

The idea

A common divisor of a and b is a number that divides both with no remainder. The common divisors of 48 and 18 are 1, 2, 3 and 6, so the greatest common divisor — the gcd — is 6.

Euclid's insight is that you can shrink the problem without changing the answer.

Take two positive numbers with a larger than b, and suppose some number d divides both. Then d divides their difference a - b too — a whole number of ds minus a whole number of ds is still a whole number of ds. The argument runs backwards just as easily: anything dividing b and a - b divides their sum, which is a.

So the pair (a, b) and the pair (a - b, b) have exactly the same set of common divisors. Not a similar set — the same one. Identical sets have identical largest members:

gcd(a, b) = gcd(a - b, b)

That is the whole algorithm. Keep replacing the larger number with the difference. The numbers shrink, the answer never moves, and when the two finally become equal you are looking at the gcd, because the greatest divisor a number shares with itself is the number.

From subtracting to dividing

Repeated subtraction is correct but wasteful: the gcd of 1,000,000 and 2 costs half a million identical subtractions.

Look at what those subtractions achieve. They strip out copies of b until what is left is smaller than b — which is division with remainder, spelled slowly. Doing it in one step gives Euclid's real rule:

gcd(a, b) = gcd(b, a mod b)

where a mod b is the remainder after dividing a by b. When the remainder hits zero, b divides a exactly, and the gcd of the pair is b.

Two properties make this a real algorithm rather than a nice identity.

It is correct at every step. Write a = q * b + r. The same divisor argument still applies, with q copies of b removed instead of one: r = a - q * b is a difference of multiples of any common divisor d, and q * b + r = a is a sum of them. Same set of common divisors, same greatest member.

It always terminates. The remainder is never negative and is always strictly smaller than b, so the second number of the pair strictly decreases every step and is bounded below by zero. A strictly decreasing sequence of non-negative integers has to reach 0. It cannot take longer than b steps; it will turn out to take dramatically fewer.

Watching it work

Start with the subtraction form on 48 and 18, since it is the version you can follow without doing any arithmetic in your head.

• 48 and 18 — 48 is larger, so replace it: 48 - 18 = 30. Now (30, 18).
• 30 and 18 — replace again: 30 - 18 = 12. Now (12, 18).
• 12 and 18 — now 18 is the larger one: 18 - 12 = 6. Now (12, 6).
• 12 and 6 — 12 - 6 = 6. Now (6, 6).
• The two are equal. The answer is 6.

Four subtractions. Every intermediate pair — (30, 18), (12, 18), (12, 6), (6, 6) — has the same common divisors as the original: 1, 2, 3 and 6.

Now the division form on a harder pair, 1071 and 462. By subtraction this takes 11 steps. By division it takes three:

• 1071 = 2 * 462 + 147, so the remainder is 147. The pair becomes (462, 147).
• 462 = 3 * 147 + 21, remainder 21. The pair becomes (147, 21).
• 147 = 7 * 21 + 0, remainder 0. Done — the gcd is 21.

Notice the shift that happens between rows: the old b becomes the new a, and the remainder becomes the new b. That single move is the entire loop body. And notice the last non-zero remainder, 21 — that is always the answer. Check it: 1071 = 21 × 51 and 462 = 21 × 22, and 51 and 22 share nothing.

The code

Here is the subtraction version first, because it maps one-to-one onto the argument above.

def gcd_subtraction(a: int, b: int) -> int:
    """Greatest common divisor by repeated subtraction.

    Both inputs must be positive integers. Replacing the larger number with
    the difference leaves the set of common divisors untouched, so the answer
    never changes; the two numbers just get smaller until they meet.
    """
    while a != b:
        if a > b:
            a -= b
        else:
            b -= a
    return a


print(gcd_subtraction(48, 18))
print(gcd_subtraction(1071, 462))

6
21

Correct, but it is only ever a teaching version — feed it 1,000,000 and 2 and it performs 499,999 subtractions, and feed it a zero and it never returns at all.

The division form fixes both problems and is shorter:

def gcd(a: int, b: int) -> int:
    """Greatest common divisor by Euclid's division rule.

    Repeatedly replaces the pair (a, b) with (b, a mod b). The second value
    strictly shrinks every step, so it reaches 0, and the first value is the
    answer at that point. gcd(0, 0) is 0 by convention, matching math.gcd.
    """
    a, b = abs(a), abs(b)
    while b:
        a, b = b, a % b
    return a


def gcd_recursive(a: int, b: int) -> int:
    """The same rule written as its own definition: gcd(a, b) = gcd(b, a mod b)."""
    if b == 0:
        return abs(a)
    return gcd_recursive(b, a % b)


print(gcd(1071, 462), gcd_recursive(1071, 462))
print(gcd(48, 18), gcd(18, 48))
print(gcd(0, 5), gcd(5, 0), gcd(0, 0))
print(gcd(-48, 18), gcd(270, 192))

21 21
6 6
5 5 0
6 6

Both versions are five lines of real work. The recursive one is the identity gcd(a, b) = gcd(b, a mod b) typed out verbatim, which makes it the better one to read; the iterative one is the better one to ship, since it uses no stack.

How the code maps to the idea

while b: is the termination test. In Python an integer is falsy only when it is 0, so this loop runs until the remainder reaches zero — exactly the stopping condition from the proof. When it exits, b is 0 and the answer sits in a.

a, b = b, a % b is the whole algorithm in one line. The right-hand side is fully evaluated before either assignment happens, so a % b still uses the old a. Split it into two statements and you destroy a before computing the remainder — that is the most common way this function gets broken.

The order of the arguments does not matter. If you call gcd(18, 48), the first iteration computes 18 mod 48, which is 18, so the pair becomes (48, 18) and the algorithm carries on normally. Passing them the wrong way round costs exactly one extra division, never a wrong answer.

abs on the way in handles negative inputs. Divisors do not care about sign — 6 divides both 48 and −48 — so the gcd is conventionally non-negative.

Zero needs no special case. gcd(5, 0) skips the loop and returns 5, which is right: every number divides 0, so the greatest divisor shared by 5 and 0 is 5. gcd(0, 5) takes one iteration to swap into the same state. gcd(0, 0) returns 0, the convention math.gcd also uses.

Now count the work instead of trusting the claim:

def count_subtraction_steps(a: int, b: int) -> int:
    """How many subtractions the naive form performs."""
    steps = 0
    while a != b:
        if a > b:
            a -= b
        else:
            b -= a
        steps += 1
    return steps


def count_division_steps(a: int, b: int) -> int:
    """How many divisions the modulo form performs."""
    steps = 0
    while b:
        a, b = b, a % b
        steps += 1
    return steps


print(f"{'pair':<20}{'subtractions':>14}{'divisions':>11}")
for first, second in [(48, 18), (1071, 462), (1_000_000, 2)]:
    pair = f"gcd({first}, {second})"
    print(f"{pair:<20}{count_subtraction_steps(first, second):>14}"
          f"{count_division_steps(first, second):>11}")

pair                  subtractions  divisions
gcd(48, 18)                      4          3
gcd(1071, 462)                  11          3
gcd(1000000, 2)             499999          1

That last row is the argument for division in one line: 499,999 steps against one.

Complexity

Time: O(log min(a, b)). Here is where the logarithm comes from.

Claim: when a is at least b, the remainder a mod b is smaller than a / 2. Two cases, and they cover everything. If b is at most a / 2, then the remainder is smaller than b, which is at most a / 2. If b is bigger than a / 2, then b fits into a exactly once, so the remainder is a - b, which is again smaller than a / 2. Either way the remainder loses at least half. (The loop guarantees the precondition: from the first division onwards, the new first number is the old second one, which is the larger of the pair.)

Now follow the pair through two steps. Starting from (a, b) you get (b, a mod b), then (a mod b, ...). So after two steps the first slot holds a value smaller than half of what it held before. A quantity that halves every two steps reaches 1 after about 2 * log2(a) steps, and the loop stops there. After the very first division both numbers are at most the smaller input, so the bound is O(log min(a, b)) divisions.

Concretely: an eighteen-digit number is about 60 bits, so this argument caps the run at roughly 120 divisions — against 10¹⁸ trial divisions for the naive search. The true bound is sharper still, and it has a name.

The worst case is consecutive Fibonacci numbers, and this is not a curiosity — it is exactly the input that makes every quotient as small as possible. Each Fibonacci number is the sum of the two before it, so dividing one by its predecessor gives a quotient of 1 and a remainder equal to the one before that. Quotient 1 is the least progress a division step can make.

print(f"{'pair':<16}{'divisions':>10}")
previous, current = 1, 1
for _ in range(8):
    previous, current = current, previous + current
    pair = f"gcd({current}, {previous})"
    print(f"{pair:<16}{count_division_steps(current, previous):>10}")

worst_steps, worst_pair = 0, (0, 0)
for larger in range(2, 1000):
    for smaller in range(1, larger):
        steps = count_division_steps(larger, smaller)
        if steps > worst_steps:
            worst_steps, worst_pair = steps, (larger, smaller)
print(f"worst pair below 1000: gcd{worst_pair} takes {worst_steps} divisions")

pair             divisions
gcd(2, 1)                1
gcd(3, 2)                2
gcd(5, 3)                3
gcd(8, 5)                4
gcd(13, 8)               5
gcd(21, 13)              6
gcd(34, 21)              7
gcd(55, 34)              8
worst pair below 1000: gcd(987, 610) takes 14 divisions

Each rung of the Fibonacci ladder costs exactly one more division than the last, and the brute-force search over every pair below 1,000 lands on 987 and 610 — Fibonacci numbers again.

This is what Lamé's theorem (1844) formalises: if the algorithm needs n division steps on a pair whose smaller member is b, then b is at least the (n+1)-th Fibonacci number. Turned around, n is never more than five times the number of decimal digits of b — so an eighteen-digit pair takes at most 90 divisions, tightening the 120 the halving argument gave. The 14-step run above therefore needed a b of at least 610, which is exactly what the search found. Lamé's proof is generally counted as the first complexity analysis of an algorithm in history, written a century before there was a computer to run one.

Space: O(1) for the iterative version — two integers and a temporary. The recursive version uses one stack frame per division, so O(log min(a, b)) frames; for any pair of 64-bit integers that is under 100 frames, comfortably inside Python's default recursion limit of 1,000.

One honest caveat: counting divisions assumes a division is one operation, which holds only while the numbers fit in a machine word. RSA keys are 2048 bits and up. Counting bit operations instead, the whole run costs about O(k²) for k-bit inputs with schoolbook division — practical, but the reason production libraries use the refinements below.

Two things you get for free

Least common multiple

The lowest common multiple of 4 and 6 is 12. There is no separate algorithm for it, because gcd and lcm are two halves of one fact:

a * b = gcd(a, b) * lcm(a, b)

The reason is visible in the prime factorisations. For each prime, the gcd takes the smaller exponent and the lcm takes the larger. The smaller plus the larger is the sum, and the sum is what the product a * b has. So the lcm is the product divided by the gcd.

def lcm(a: int, b: int) -> int:
    """Lowest common multiple, derived from the gcd.

    Divides before multiplying so the intermediate value never exceeds the
    answer, which matters in every language with fixed-width integers.
    """
    if a == 0 or b == 0:
        return 0
    return abs(a // gcd(a, b) * b)


print(lcm(4, 6), lcm(21, 6), lcm(1071, 462))
print(lcm(0, 7))

12 42 23562
0

Write it as a // gcd(a, b) * b, not a * b // gcd(a, b). The division is exact either way, but dividing first keeps the intermediate value no larger than the answer. In Python that is a habit; in C or Java it is the difference between a correct result and a wrapped one. The zero guard matters too — without it, lcm(0, 0) divides by a gcd of zero.

The extended Euclidean algorithm

Bézout's identity says something stronger than "the gcd exists": for any a and b there are integers x and y with

a * x + b * y = gcd(a, b)

Those coefficients are already hiding in the division trace — you just have to read it backwards. From the run on 1071 and 462:

• 21 = 462 - 3 * 147 (the second division, rearranged)
• 147 = 1071 - 2 * 462 (the first division, rearranged)

Substitute the second into the first: 21 = 462 - 3 * (1071 - 2 * 462) = 7 * 462 - 3 * 1071. So x = -3 and y = 7. Multiply it out if you like: −3,213 + 3,234 = 21.

The recursive implementation does that substitution automatically. Each call gets its child's coefficients and rewrites them in terms of its own two numbers.

def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return (g, x, y) such that a * x + b * y == g == gcd(a, b).

    The base case is free: gcd(a, 0) is a, and a * 1 + 0 * 0 == a. Every
    other level rewrites its child's answer in terms of its own two numbers.
    """
    if b == 0:
        return a, 1, 0
    g, x, y = extended_gcd(b, a % b)
    # The child solved b * x + (a mod b) * y == g. Substituting
    # a mod b == a - (a // b) * b and collecting terms gives the line below.
    return g, y, x - (a // b) * y


g, x, y = extended_gcd(1071, 462)
print(f"gcd = {g}, x = {x}, y = {y}")
print(f"1071 * ({x}) + 462 * {y} = {1071 * x + 462 * y}")

gcd = 21, x = -3, y = 7
1071 * (-3) + 462 * 7 = 21

The recursion descends exactly as far as the plain version does — same divisions, same bound, three extra multiplications per level. Extended Euclid is O(log min(a, b)) too.

Modular inverses, and why RSA needs them

This is the payoff. In ordinary arithmetic the inverse of 17 is 1/17. Modular arithmetic has no fractions, so the inverse of a modulo m is the whole number x with a * x leaving a remainder of 1 when divided by m.

Such an x exists if and only if gcd(a, m) = 1, and the reason is Bézout. If the gcd is 1, then a * x + m * y = 1 for some x and y, and m * y vanishes modulo m — so a * x leaves remainder 1. If the gcd is bigger than 1, no multiple of a can ever land on 1, because everything in sight is divisible by that common factor.

So the extended algorithm is the modular inverse algorithm. Take x, reduce it into range, done.

def modular_inverse(value: int, modulus: int) -> int:
    """The number that multiplies `value` back to 1, working modulo `modulus`."""
    g, x, _ = extended_gcd(value % modulus, modulus)
    if g != 1:
        raise ValueError(f"{value} has no inverse modulo {modulus}")
    return x % modulus


public_exponent = 17
totient = 3120  # (61 - 1) * (53 - 1) for the textbook RSA key with n = 3233
private_exponent = modular_inverse(public_exponent, totient)
print(f"private exponent d = {private_exponent}")
print(f"17 * {private_exponent} mod 3120 = {(public_exponent * private_exponent) % totient}")

try:
    modular_inverse(6, 9)
except ValueError as error:
    print(error)

private exponent d = 2753
17 * 2753 mod 3120 = 1
6 has no inverse modulo 9

Those numbers are the textbook RSA key: primes 61 and 53, modulus 3,233, public exponent 17. The private exponent 2,753 is not chosen or searched for — it is the modular inverse of 17, produced by extended Euclid in four steps. Generating an RSA key pair means picking a public exponent coprime to the totient (a gcd test) and then inverting it (extended Euclid). Both halves are this post.

When to use it, and when not to

Use it whenever you need a gcd, an lcm, or a modular inverse — in practice, fractions, ratios, cycle lengths and any modular arithmetic. It is one of the few algorithms with no real competition in its weight class.

But in Python, call the standard library. math.gcd is written in C, accepts any number of arguments since 3.9, and for large integers CPython switches to Lehmer's algorithm — a refinement that handles several quotients at a time from the leading digits alone, cutting the number of full-width divisions sharply. math.lcm arrived in 3.9, and since 3.8 pow(a, -1, m) gives you a modular inverse directly.

import math
from fractions import Fraction

print(math.gcd(1071, 462), math.lcm(4, 6, 10), math.gcd(12, 18, 24))
print(Fraction(1071, 462))
print(pow(17, -1, 3120))

for width, height in [(1920, 1080), (2560, 1440), (1366, 768)]:
    divisor = math.gcd(width, height)
    print(f"{width}x{height} -> {width // divisor}:{height // divisor}")

21 60 6
51/22
2753
1920x1080 -> 16:9
2560x1440 -> 16:9
1366x768 -> 683:384

Write it yourself when you need the Bézout coefficients — pow(a, -1, m) hands back the inverse but throws away x and y, which you need for solving equations of the form a * x + b * y = c. Or in a language with no gcd in its standard library. Or when you are learning, which is the best reason of all.

Do not use the subtraction form. Its only defence is hardware where division is unavailable or ruinously expensive, and even there the right answer is binary GCD (Stein's algorithm), which replaces division with subtraction, halving and parity tests.

Do not expect it to factorise anything. The gcd tells you what two numbers share, not what either one is made of. Finding the prime factors of a single number is a vastly harder problem — the one RSA's security rests on.

Where it shows up in the real world

Python's fractions.Fraction normalises every fraction on construction by dividing numerator and denominator by their gcd. That is why Fraction(1071, 462) printed 51/22 above. CPython's statistics module adds its partial sums as Fraction objects to stay exact, so gcds run on the way to a plain mean too.

RSA key generation. Choosing a public exponent means checking that it is coprime to (p - 1) * (q - 1), which is a gcd test, and the private exponent is that public exponent's modular inverse, which is extended Euclid. Both halves of generating a key pair are this algorithm.

pow(a, -1, m) in CPython is an extended-Euclid implementation, and the same routine underpins modular division across cryptographic libraries — including the Chinese remainder theorem shortcut that makes RSA decryption several times faster.

Aspect ratios. Dividing a screen's width and height by their gcd is how 1920×1080 becomes 16:9. It also exposes the odd ones: 1366×768 reduces no further than 683:384, because 683 is prime. A true 16:9 panel 768 pixels tall would be 1365.33 pixels wide, so that laptop resolution is famously almost 16:9 and never exactly.

Anything with repeating cycles. Two signals with periods 12 and 18 realign after lcm(12, 18) = 36 units. Gear ratios, polyrhythms in music, and traffic-light cycle planning are the same computation with different units.

Common mistakes

Passing zero to the subtraction version. With a = 0 and b = 5, the two are never equal, and b -= a subtracts nothing forever. The loop hangs. The division form has no such failure mode.

Returning the wrong variable. After while b: ends, b is 0 by definition. The answer is a. Returning b gives you a very confident zero.

Splitting the swap. Writing a = b followed by b = a % b computes the remainder of b divided by b, which is 0, so the function returns after one step with the wrong answer. Keep it as a single tuple assignment.

Porting the modulo naively. Python's % returns a non-negative result when the divisor is positive; C, Java, Go and JavaScript take the sign of the dividend instead, so -48 % 18 is −12 in those languages and 6 in Python. A direct port can return a negative gcd, and the subtraction form can loop forever. Take absolute values on entry.

Forgetting to reduce the Bézout coefficient. The x that comes back is often negative — for 17 modulo 3120 it is −367. A modular inverse must be reported in range, so return x % modulus.

Assuming an inverse always exists. It exists only when the gcd is 1. Check the returned g before using x, as modular_inverse does above; skipping the check gives you a number that quietly fails to invert anything.

Practice

1. Write gcd_of_list that returns the gcd of any number of integers, using functools.reduce over the two-argument version.
2. Reduce a fraction to lowest terms without using fractions, making sure a negative denominator moves its sign to the numerator.
3. Count the divisions for every pair below 10,000 and confirm that the worst pair is again two consecutive Fibonacci numbers.
4. Use extended_gcd to solve a * x + b * y = c for given a, b and c, reporting that no solution exists when the gcd does not divide c.
5. Implement binary GCD (Stein's algorithm) using only subtraction, halving and even/odd tests, and compare its step count against the division form on a few large pairs.

Summary

The Euclidean algorithm is the best return on five lines of code in all of computing. One observation — that a and b share their divisors with b and a mod b — turns a search over a quintillion candidates into fewer than a hundred divisions, and the extended version throws in Bézout coefficients and modular inverses at no extra asymptotic cost. Twenty-three centuries later, nothing has replaced it.

Keep reading

• Fast Exponentiation — the other half of modular arithmetic, and how RSA actually encrypts once it has these keys.
• Sieve of Eratosthenes — the other ancient algorithm still in daily use, for finding the primes this one needs.
• Big O Notation — the full treatment of the counting arguments used above.