Euclid's algorithm: GCD, LCM and Bézout

The greatest common divisor (GCD) of two whole numbers, also called the highest common factor (HCF), is the largest number that divides both of them: gcd(1920, 1080) = 120. Their least common multiple (LCM) is the smallest positive number that both of them divide: lcm(4, 6) = 12. Euclid’s algorithm repeats a single step, turning (a, b) into (b, a mod b), and stops at a remainder of 0; it needs only O(log min(a, b)) divisions, and consecutive Fibonacci numbers are its slowest inputs. The extended version also finds integers x and y with ax + by = gcd(a, b), the Bézout coefficients, which solve equations such as 7x + 11y = 100 and give the modular inverses used in the next lessons. The LCM then follows from the GCD as a / gcd(a, b) × b, dividing first so that no intermediate value overflows.

Euclid’s rule: gcd(a, b) = gcd(b, a mod b)

Divide a by b: a = q × b + r, with a remainder r from 0 to b − 1. Any number that divides a and b also divides r = a − q × b, and any number that divides b and r also divides a = q × b + r. The pairs (a, b) and (b, r) therefore have exactly the same common divisors, and in particular the same largest one. Since r < b, every step lowers the second number of the pair, which has to reach 0 sooner or later; the first number is then the answer, because every number divides 0 and gcd(g, 0) = g. Book VII of Euclid’s Elements gives the method as taking the smaller number away from the larger one again and again; one division does many of those subtractions at once.

A screen’s aspect ratio is its width and its height divided by their GCD. This program prints every division for three screen sizes:

import math


def gcd_steps(a, b):
    """Euclid's algorithm: returns the GCD and the division steps a = q × b + r."""
    steps = []
    while b != 0:
        q, r = divmod(a, b)
        steps.append((a, q, b, r))
        a, b = b, r
    return a, steps


for width, height in [(1920, 1080), (2560, 1080), (3440, 1440)]:
    g, steps = gcd_steps(width, height)
    print(f"{width} × {height}")
    for a, q, b, r in steps:
        print(f"  {a:>5} = {q} × {b:>4} + {r:>4}")
    print(f"  gcd = {g}, aspect ratio {width // g}:{height // g}, math.gcd agrees: {g == math.gcd(width, height)}")

1920 × 1080
   1920 = 1 × 1080 +  840
   1080 = 1 ×  840 +  240
    840 = 3 ×  240 +  120
    240 = 2 ×  120 +    0
  gcd = 120, aspect ratio 16:9, math.gcd agrees: True
2560 × 1080
   2560 = 2 × 1080 +  400
   1080 = 2 ×  400 +  280
    400 = 1 ×  280 +  120
    280 = 2 ×  120 +   40
    120 = 3 ×   40 +    0
  gcd = 40, aspect ratio 64:27, math.gcd agrees: True
3440 × 1440
   3440 = 2 × 1440 +  560
   1440 = 2 ×  560 +  320
    560 = 1 ×  320 +  240
    320 = 1 ×  240 +   80
    240 = 3 ×   80 +    0
  gcd = 80, aspect ratio 43:18, math.gcd agrees: True

1920 × 1080 reduces to 16:9 after four divisions, and the two wider screens reduce to 64:27 and 43:18, about 2.37 and 2.39 to 1. The order of the two numbers does not matter: when a is smaller than b, the first step gives a mod b = a and only swaps them.

The picture shows the same four steps as geometry. Cutting the largest squares that fit from a 1920 × 1080 rectangle leaves a strip of 840 × 1080, then 840 × 240, then 120 × 240, which two squares of 120 fill exactly. The side of the last square is the GCD: squares of 120 tile the whole rectangle, and no larger square does.

Zero and negative numbers have fixed answers. Python’s math.gcd never returns a negative number, returns 0 when all its arguments are 0, and takes any number of arguments since Python 3.9; C++17’s std::gcd and Java’s BigInteger.gcd also work on the absolute values of their arguments.

How many steps: Fibonacci pairs are the worst case

The numbers shrink quickly. Take a ≥ b > 0. If b is at most a / 2, the remainder is smaller than b and so smaller than a / 2. If b is larger than a / 2, the quotient is 1 and the remainder is a − b, again smaller than a / 2. Either way, the first number of the pair has at least halved two steps later, so a pair whose smaller number is b needs about 2 log₂ b divisions at most: O(log min(a, b)).

The slowest pairs are consecutive Fibonacci numbers, such as 55 and 34. Every quotient is 1, so each remainder is only the difference of the two numbers before it, the smallest drop there is, and the remainders walk down the Fibonacci sequence one number at a time. Gabriel Lamé proved a bound from this: the number of divisions never exceeds five times the number of decimal digits of the smaller number. The program counts the divisions for the largest Fibonacci pair of each size, and for 10,000 random pairs whose smaller number has the same number of digits:

import random


def steps(a, b):
    """The number of divisions Euclid's algorithm makes on (a, b)."""
    count = 0
    while b != 0:
        a, b = b, a % b
        count += 1
    return count


fib = [0, 1]
while fib[-1] < 10**19:
    fib.append(fib[-1] + fib[-2])

rng = random.Random(2026)
print("digits  Fibonacci pair: steps   random pairs: average  max   5 × digits")
for d in (3, 6, 9, 12, 15, 18):
    # the largest Fibonacci number with d digits, and the one after it
    n = max(i for i, f in enumerate(fib) if len(str(f)) == d)
    fib_steps = steps(fib[n + 1], fib[n])
    counts = []
    for _ in range(10_000):
        b = rng.randrange(10 ** (d - 1), 10**d)  # the smaller number has d digits
        a = rng.randrange(b + 1, 10 ** (d + 1))
        counts.append(steps(a, b))
    print(f"{d:>6}  F({n + 1}), F({n}): {fib_steps:>5}   {sum(counts) / len(counts):>21.1f}  {max(counts):>3}   {5 * d:>10}")

digits  Fibonacci pair: steps   random pairs: average  max   5 × digits
     3  F(17), F(16):    15                     6.2   13           15
     6  F(31), F(30):    29                    12.0   21           30
     9  F(45), F(44):    43                    17.8   29           45
    12  F(60), F(59):    58                    23.6   37           60
    15  F(74), F(73):    72                    29.5   48           75
    18  F(88), F(87):    86                    35.4   52           90

The Fibonacci pairs sit at Lamé’s bound: 15 divisions for 3 digits, where the bound allows 15, and 86 for 18 digits, where it allows 90. Random pairs of the same size need well under half as many, 35.4 on average for 18 digits, and the slowest of the 10,000 took 52. Even the worst case is small: two numbers below 10¹⁸ never need more than 90 divisions. For numbers with thousands of digits each division costs more, and NIST’s dictionary puts the whole algorithm at O((log a)(log b)) bit operations.

Extended Euclid: the x and y with ax + by = gcd(a, b)

The same divisions also show how to build the GCD out of the two numbers. Write every remainder r as a combination r = a·s + b·t of the original pair. The first two rows need no work: a = a·1 + b·0 and b = a·0 + b·1. Each new remainder is the row before minus q times the current row, and doing the same subtraction to the s and t columns keeps every row true, because the combination is linear. When the remainder reaches 0, the row above holds the GCD with its coefficients:

import random


def extended_gcd(a, b):
    """Returns (g, x, y) with a·x + b·y = g = gcd(a, b), for a, b >= 0 not both 0."""
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    while r != 0:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    return old_r, old_s, old_t


def trace(a, b):
    """Prints every remainder r with the s and t that keep r = a·s + b·t."""
    print(f"{'r':>4} {'q':>3} {'s':>4} {'t':>4}   {a}·s + {b}·t")
    rows = [(a, None, 1, 0), (b, None, 0, 1)]
    while rows[-1][0] != 0:
        (r0, _, s0, t0), (r1, _, s1, t1) = rows[-2], rows[-1]
        q = r0 // r1
        rows[-1] = (r1, q, s1, t1)
        rows.append((r0 - q * r1, None, s0 - q * s1, t0 - q * t1))
    for r, q, s, t in rows:
        print(f"{r:>4} {'' if q is None else q:>3} {s:>4} {t:>4}   {a * s + b * t:>8}")


trace(11, 7)
g, x, y = extended_gcd(11, 7)
print(f"gcd = {g} = 11·({x}) + 7·({y})")
print(f"so 7 has the inverse {y % 11} modulo 11: 7·{y % 11} = {7 * (y % 11)} = {7 * (y % 11) // 11}·11 + 1;",
      f"pow(7, -1, 11) = {pow(7, -1, 11)}")

rng = random.Random(7)
worst = 0.0
for _ in range(10_000):
    a, b = rng.randrange(1, 10**12), rng.randrange(1, 10**12)
    g, x, y = extended_gcd(a, b)
    assert a * x + b * y == g and a % g == 0 and b % g == 0
    worst = max(worst, abs(x) * g / b, abs(y) * g / a)
print(f"10,000 random pairs below 10^12: a·x + b·y = g every time; largest |x|/(b/g) or |y|/(a/g): {worst:.5f}")

   r   q    s    t   11·s + 7·t
  11        1    0         11
   7   1    0    1          7
   4   1    1   -1          4
   3   1   -1    2          3
   1   3    2   -3          1
   0       -7   11          0
gcd = 1 = 11·(2) + 7·(-3)
so 7 has the inverse 8 modulo 11: 7·8 = 56 = 5·11 + 1; pow(7, -1, 11) = 8
10,000 random pairs below 10^12: a·x + b·y = g every time; largest |x|/(b/g) or |y|/(a/g): 0.49997

Read the table from the top. The q of a row is how many times that row fits into the row above it, and every row matches its last column: 4 = 11·1 + 7·(−1), 3 = 11·(−1) + 7·2 and finally 1 = 11·2 + 7·(−3). These x and y are the Bézout coefficients. No positive combination ax + by can be smaller than gcd(a, b), since the GCD divides every combination, and the algorithm reaches the GCD itself, so gcd(a, b) is the smallest positive value that ax + by can take: that is Bézout’s identity. The coefficients also stay small: over 10,000 random pairs, x was never more than half of b / g in size, nor y more than half of a / g.

The coefficients matter most for inverses. When gcd(a, m) = 1, the equation a·x + m·y = 1 says that a·x leaves the remainder 1 when divided by m, so x is the inverse of a modulo m, the number that undoes a multiplication by a. Here 7·(−3) leaves remainder 1 modulo 11, and −3 mod 11 = 8: 7 × 8 = 56 = 5 × 11 + 1. Python’s pow(7, -1, 11) returns the same 8, a feature added in Python 3.8, and the lesson on modular arithmetic builds on it.

Solving ax + by = c in integers

A school buys pens at 7 each and notebooks at 11 each and wants to spend exactly 100. In whole numbers x ≥ 0 and y ≥ 0, that is 7x + 11y = 100. Three facts solve every equation of this kind:

1. An integer solution exists exactly when g = gcd(a, b) divides c. Every value of ax + by is a multiple of g, so g must divide c; when it does, the Bézout coefficients multiplied by c / g give a solution.
2. One solution gives all the others. From (x₀, y₀), every solution is x = x₀ + k·(b / g), y = y₀ − k·(a / g) for an integer k: raising x by b / g adds a·b / g to the left side, and lowering y by a / g takes the same amount away.
3. Bounds become a range of k. Conditions such as x ≥ 0 and y ≥ 0 give a lowest and a highest k.

def extended_gcd(a, b):
    """Returns (g, x, y) with a·x + b·y = g = gcd(a, b)."""
    old_r, r, old_s, s, old_t, t = a, b, 1, 0, 0, 1
    while r != 0:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    return old_r, old_s, old_t


def solve(a, b, c):
    """Prints every way to write c = a·x + b·y with whole numbers x, y >= 0 (a, b > 0)."""
    g, x, y = extended_gcd(a, b)
    print(f"{a}x + {b}y = {c}: gcd({a}, {b}) = {g} = {a}·({x}) + {b}·({y})")
    if c % g != 0:
        print(f"  {g} does not divide {c}: no whole-number solution at all")
        return
    x0, y0 = x * (c // g), y * (c // g)
    dx, dy = b // g, a // g
    print(f"  one solution: x = {x0}, y = {y0}; all of them: x = {x0} + {dx}k, y = {y0} - {dy}k")
    k_low = -(x0 // dx)  # the smallest k with x >= 0: ceil(-x0 / dx)
    k_high = y0 // dy  # the largest k with y >= 0: floor(y0 / dy)
    found = [(x0 + dx * k, y0 - dy * k) for k in range(k_low, k_high + 1)]
    print(f"  x >= 0 and y >= 0 need {k_low} <= k <= {k_high}: {found if found else 'none'}")


solve(7, 11, 100)  # pens at 7 and notebooks at 11, exactly 100 spent
solve(7, 11, 200)
solve(6, 10, 8)
solve(6, 10, 7)

7x + 11y = 100: gcd(7, 11) = 1 = 7·(-3) + 11·(2)
  one solution: x = -300, y = 200; all of them: x = -300 + 11k, y = 200 - 7k
  x >= 0 and y >= 0 need 28 <= k <= 28: [(8, 4)]
7x + 11y = 200: gcd(7, 11) = 1 = 7·(-3) + 11·(2)
  one solution: x = -600, y = 400; all of them: x = -600 + 11k, y = 400 - 7k
  x >= 0 and y >= 0 need 55 <= k <= 57: [(5, 15), (16, 8), (27, 1)]
6x + 10y = 8: gcd(6, 10) = 2 = 6·(2) + 10·(-1)
  one solution: x = 8, y = -4; all of them: x = 8 + 5k, y = -4 - 3k
  x >= 0 and y >= 0 need -1 <= k <= -2: none
6x + 10y = 7: gcd(6, 10) = 2 = 6·(2) + 10·(-1)
  2 does not divide 7: no whole-number solution at all

The extended algorithm gives 7·(−3) + 11·2 = 1, so x = −300 and y = 200 spend exactly 100, with a negative number of pens. Stepping k moves to x = −300 + 11k and y = 200 − 7k, and both are at least 0 only for k = 28: 8 pens and 4 notebooks are the one way to spend exactly 100. With 200 to spend there are three ways. The last two equations fail in two different ways: 6x + 10y = 8 has integer solutions, but every one of them has a negative part, while 6x + 10y = 7 has none at all, because the left side is always even.

LCM: divide before you multiply

For each prime, the GCD of two positive numbers takes the smaller of its two powers and the LCM the larger one, so together they take both: gcd(a, b) × lcm(a, b) = a × b. The formula a × b / gcd(a, b) is right on paper and risky in code, because the product can be far larger than the LCM. Dividing first, a / gcd(a, b) × b, is exact, since the GCD divides a, and nothing it computes is larger than the result. The same two pairs in five languages:

import math

CASES = [(6_000_000_006, 4_000_000_004), (99_999_989, 99_999_971)]

for a, b in CASES:
    g = math.gcd(a, b)
    print(f"a = {a}, b = {b}, gcd = {g}")
    print(f"  a * b / gcd  = {a * b // g}")
    print(f"  a / gcd * b  = {a // g * b}")
    print(f"  math.lcm     = {math.lcm(a, b)}")
print(f"math.lcm of 1..30 = {math.lcm(*range(1, 31))}")

a = 6000000006, b = 4000000004, gcd = 2000000002
  a * b / gcd  = 12000000012
  a / gcd * b  = 12000000012
  math.lcm     = 12000000012
a = 99999989, b = 99999971, gcd = 1
  a * b / gcd  = 9999996000000319
  a / gcd * b  = 9999996000000319
  math.lcm     = 9999996000000319
math.lcm of 1..30 = 2329089562800

// JavaScript has no gcd built in; Euclid's loop works on Numbers and on BigInts alike.
function gcd(a, b) {
  while (b) [a, b] = [b, a % b];
  return a;
}

const CASES = [[6000000006, 4000000004], [99999989, 99999971]];

for (const [a, b] of CASES) {
  const g = gcd(a, b);
  console.log(`a = ${a}, b = ${b}, gcd = ${g}`);
  console.log(`  a * b / gcd  = ${(a * b) / g}`);
  console.log(`  a / gcd * b  = ${(a / g) * b}`);
  const [A, B] = [BigInt(a), BigInt(b)];
  console.log(`  with BigInt  = ${(A / gcd(A, B)) * B}`);
}
console.log(`largest exact integer of a Number: ${Number.MAX_SAFE_INTEGER}`);

a = 6000000006, b = 4000000004, gcd = 2000000002
  a * b / gcd  = 12000000012
  a / gcd * b  = 12000000012
  with BigInt  = 12000000012
a = 99999989, b = 99999971, gcd = 1
  a * b / gcd  = 9999996000000320
  a / gcd * b  = 9999996000000320
  with BigInt  = 9999996000000319
largest exact integer of a Number: 9007199254740991

import java.math.BigInteger;

public class Lcm {
    static long gcd(long a, long b) {
        while (b != 0) {
            long r = a % b;
            a = b;
            b = r;
        }
        return a;
    }

    public static void main(String[] args) {
        long[][] cases = {{6_000_000_006L, 4_000_000_004L}, {99_999_989L, 99_999_971L}};
        for (long[] c : cases) {
            long a = c[0], b = c[1], g = gcd(a, b);
            System.out.println("a = " + a + ", b = " + b + ", gcd = " + g);
            System.out.println("  a * b / gcd  = " + a * b / g);
            System.out.println("  a / gcd * b  = " + a / g * b);
            try {
                System.out.println("  multiplyExact(a, b) = " + Math.multiplyExact(a, b));
            } catch (ArithmeticException e) {
                System.out.println("  multiplyExact(a, b) throws ArithmeticException: " + e.getMessage());
            }
            BigInteger A = BigInteger.valueOf(a), B = BigInteger.valueOf(b);
            System.out.println("  BigInteger   = " + A.divide(A.gcd(B)).multiply(B));
        }
    }
}

#include <cstdint>
#include <iostream>
#include <numeric>

int main() {
    const std::int64_t cases[][2] = {{6'000'000'006, 4'000'000'004}, {99'999'989, 99'999'971}};
    for (const auto& c : cases) {
        const std::int64_t a = c[0], b = c[1], g = std::gcd(a, b);
        std::cout << "a = " << a << ", b = " << b << ", gcd = " << g << "\n";
        std::int64_t product;
        const bool overflow = __builtin_mul_overflow(a, b, &product);
        std::cout << "  a * b fits in 64 bits: " << (overflow ? "no" : "yes") << "\n";
        std::cout << "  a / gcd * b  = " << a / g * b << "\n";
        std::cout << "  std::lcm     = " << std::lcm(a, b) << "\n";
    }
}

package main

import (
	"fmt"
	"math/big"
)

func gcd(a, b int64) int64 {
	for b != 0 {
		a, b = b, a%b
	}
	return a
}

func main() {
	cases := [][2]int64{{6_000_000_006, 4_000_000_004}, {99_999_989, 99_999_971}}
	for _, c := range cases {
		a, b := c[0], c[1]
		g := gcd(a, b)
		fmt.Printf("a = %d, b = %d, gcd = %d\n", a, b, g)
		fmt.Printf("  a * b / gcd  = %d\n", a*b/g)
		fmt.Printf("  a / gcd * b  = %d\n", a/g*b)
		x, y := new(big.Int), new(big.Int)
		A, B := big.NewInt(a), big.NewInt(b)
		G := new(big.Int).GCD(x, y, A, B)
		fmt.Printf("  big.Int: gcd = %v = a*(%v) + b*(%v), lcm = %v\n", G, x, y, new(big.Int).Mul(new(big.Int).Quo(A, G), B))
	}
}

from session import run, versions

# Lcm.java runs as a single source file, lcm.cpp is built with the sanitizers on, go run builds lcm.go.
versions("java -version", "clang++ --version", "go version")
run("""
java -XX:-UsePerfData Lcm.java
clang++ -std=c++23 -O2 -Wall -Wextra -fsanitize=address,undefined lcm.cpp -o lcm && ./lcm
go run lcm.go
""")

# openjdk version "25.0.4.1" LTS
# Apple clang version 21.0.0 (clang-2100.3.34.2)
# go version go1.27.1 darwin/arm64
$ java -XX:-UsePerfData Lcm.java
a = 6000000006, b = 4000000004, gcd = 2000000002
  a * b / gcd  = 2776627984
  a / gcd * b  = 12000000012
  multiplyExact(a, b) throws ArithmeticException: long overflow
  BigInteger   = 12000000012
a = 99999989, b = 99999971, gcd = 1
  a * b / gcd  = 9999996000000319
  a / gcd * b  = 9999996000000319
  multiplyExact(a, b) = 9999996000000319
  BigInteger   = 9999996000000319
$ clang++ -std=c++23 -O2 -Wall -Wextra -fsanitize=address,undefined lcm.cpp -o lcm && ./lcm
a = 6000000006, b = 4000000004, gcd = 2000000002
  a * b fits in 64 bits: no
  a / gcd * b  = 12000000012
  std::lcm     = 12000000012
a = 99999989, b = 99999971, gcd = 1
  a * b fits in 64 bits: yes
  a / gcd * b  = 9999996000000319
  std::lcm     = 9999996000000319
$ go run lcm.go
a = 6000000006, b = 4000000004, gcd = 2000000002
  a * b / gcd  = 2776627984
  a / gcd * b  = 12000000012
  big.Int: gcd = 2000000002 = a*(1) + b*(-1), lcm = 12000000012
a = 99999989, b = 99999971, gcd = 1
  a * b / gcd  = 9999996000000319
  a / gcd * b  = 9999996000000319
  big.Int: gcd = 1 = a*(5555554) + b*(-5555555), lcm = 9999996000000319

In the first pair, a × b is about 2.4 × 10¹⁹, beyond 2⁶³ − 1 = 9,223,372,036,854,775,807, the top of a signed 64-bit integer. Java and Go wrap around without a word, divide the wrong product and print 2,776,627,984. Java’s Math.multiplyExact throws an ArithmeticException instead, and C++’s __builtin_mul_overflow reports that the product does not fit. Dividing first gives 12,000,000,012 in every language, and so do C++17’s std::lcm and Python’s math.lcm. The second pair breaks JavaScript instead. A Number is exact for whole numbers no larger than 9,007,199,254,740,991, which is 2⁵³ − 1, and the LCM 9,999,996,000,000,319 lies beyond it: it comes out as 9,999,996,000,000,320, one too large, with no error at all, while a BigInt keeps it exact. Go’s big.Int method GCD returns the Bézout coefficients as well.

Each language offers a different part of this toolkit:

The LCM of many numbers grows quickly: the LCM of 1 to 30 is already 2,329,089,562,800. Fold a list one number at a time, lcm(lcm(a, b), c) and so on, and in a 64-bit language check every step for overflow or move to big integers.

Exercise

def bezout(a, b):
    """(g, x, y) with a * x + b * y == g == gcd(a, b), for whole numbers a, b >= 0, not both 0."""
    raise NotImplementedError("write bezout")


def solutions(a, b, c, max_x, max_y):
    """Every (x, y) with a * x + b * y == c, 0 <= x <= max_x and 0 <= y <= max_y, sorted by x."""
    raise NotImplementedError("write solutions")

import math
import random

from bezout import bezout, solutions


def brute(a, b, c, max_x, max_y):
    """Tries every x: right, and far too slow for large numbers."""
    out = []
    for x in range(max_x + 1):
        rest = c - a * x
        if rest >= 0 and rest % b == 0 and rest // b <= max_y:
            out.append((x, rest // b))
    return out


def test_bezout_small():
    """small pairs in both orders, equal numbers and a zero"""
    for a, b in [(240, 46), (46, 240), (7, 11), (11, 7), (5, 5), (1, 1), (0, 7), (7, 0), (12, 18)]:
        g, x, y = bezout(a, b)
        assert g == math.gcd(a, b), (a, b)
        assert a * x + b * y == g, (a, b, x, y)


def test_bezout_large():
    """consecutive Fibonacci numbers near 10^18, the slowest case of Euclid's algorithm"""
    a, b = 1100087778366101931, 679891637638612258
    g, x, y = bezout(a, b)
    assert g == 1 and a * x + b * y == 1


def test_bezout_random():
    """500 random pairs below 10^12"""
    rng = random.Random(11)
    for _ in range(500):
        a, b = rng.randrange(1, 10**12), rng.randrange(1, 10**12)
        g, x, y = bezout(a, b)
        assert g == math.gcd(a, b) and a * x + b * y == g, (a, b)


def test_pens_and_notebooks():
    """the lesson's examples, then the same equations with tighter bounds"""
    assert solutions(7, 11, 100, 100, 100) == [(8, 4)]
    assert solutions(7, 11, 200, 100, 100) == [(5, 15), (16, 8), (27, 1)]
    assert solutions(7, 11, 200, 20, 100) == [(5, 15), (16, 8)]
    assert solutions(7, 11, 200, 100, 10) == [(16, 8), (27, 1)]


def test_no_solution():
    """the gcd does not divide c, or every whole-number solution has a negative part"""
    assert solutions(6, 10, 7, 100, 100) == []
    assert solutions(6, 10, 8, 100, 100) == []


def test_zero_and_tight_bounds():
    """c = 0, and bounds of 0"""
    assert solutions(4, 6, 0, 5, 5) == [(0, 0)]
    assert solutions(4, 6, 12, 0, 5) == [(0, 2)]
    assert solutions(4, 6, 12, 5, 0) == [(3, 0)]
    assert solutions(4, 6, 12, 0, 0) == []


def test_random_against_brute_force():
    """400 random small cases, compared with trying every x"""
    rng = random.Random(3)
    for _ in range(400):
        a, b = rng.randint(1, 30), rng.randint(1, 30)
        c, max_x, max_y = rng.randint(0, 300), rng.randint(0, 40), rng.randint(0, 40)
        assert solutions(a, b, c, max_x, max_y) == brute(a, b, c, max_x, max_y), (a, b, c, max_x, max_y)


def test_near_ten_to_the_eighteen():
    """two primes near 10^9, a c near 10^18 with exactly one solution, and a case with 10,001 solutions"""
    a, b = 999_999_937, 999_999_929
    c = a * 123_456_789 + b * 987_654_321
    assert solutions(a, b, c, 10**18, 10**18) == [(123_456_789, 987_654_321)]
    many = solutions(1, 1, 10_000, 10**18, 10**18)
    assert len(many) == 10_001 and many[0] == (0, 10_000) and many[-1] == (10_000, 0)

Interview questions

Fresher · Why does Euclid’s algorithm always stop, and which inputs make it slowest?

Each step replaces (a, b) with (b, a mod b), and the remainder is smaller than b, so the second number falls at every step; whole numbers cannot fall for ever, so it reaches 0. The first number is then the GCD, because the pair keeps the same common divisors at every step. It also stops quickly: two steps at least halve the first number, so it needs O(log min(a, b)) divisions. The slowest inputs are consecutive Fibonacci numbers, where every quotient is 1 and the remainders fall as slowly as they can. Lamé’s bound says that the number of divisions never exceeds five times the number of decimal digits of the smaller number, and the pair 1,597 and 987 needs 15 divisions, exactly that bound for three digits.

Key takeaways

• gcd(a, b) = gcd(b, a mod b) and gcd(g, 0) = g: the remainders shrink, and the last one before 0 is the GCD.
• Two steps at least halve the numbers, so Euclid needs O(log min(a, b)) divisions. Consecutive Fibonacci numbers are the slowest pairs, and no pair needs more than five divisions per digit of its smaller number.
• The extended algorithm keeps r = a·s + b·t on every row and returns x and y with ax + by = gcd(a, b); when gcd(a, m) = 1, x is the inverse of a modulo m.
• For ax + by = c to have integer solutions, gcd(a, b) must divide c; from one solution the others step by b / g and a / g, and bounds on x and y become a range of k.
• Compute lcm(a, b) as a / gcd(a, b) × b. Multiplying first overflowed in Java and Go and lost the last digit in JavaScript.