Extended Euclidean Algorithm, in depth: the invariant behind Bezout coefficients, a full trace, sign handling, overflow bounds, modular and batch inverses, and CRT

The idea: carry the combination along

Plain Euclid rests on one fact: gcd(a, b) = gcd(b, a mod b), because any number dividing a and b also divides a - q*b, and the reverse. Repeat until the second number is zero and the first is the gcd. For 240 and 46: 240 = 5*46 + 10, 46 = 4*10 + 6, 10 = 1*6 + 4, 6 = 1*4 + 2, 4 = 2*2 + 0, so the gcd is 2.

The extension adds bookkeeping. Every remainder in that sequence is built from a and b by subtraction, so each one can be written as s*a + t*b for some integers s and t. Track s and t alongside the remainder. Start with two rows whose combinations are obvious: a = 1*a + 0*b and b = 0*a + 1*b. Each new remainder is the row two above minus q times the row above, where q is the integer quotient of those two remainders. Apply that same operation to the s and t columns. Because subtraction of linear combinations gives a linear combination, the invariant r = s*a + t*b holds on every row. When r reaches zero, the previous row holds the gcd and its coefficients.

That is the whole proof of correctness: the remainder column is exactly plain Euclid, so it ends at the gcd, and the invariant says the s and t on that row satisfy Bezout's identity.

A full trace by hand

Check the gcd row: -9*240 + 47*46 = -2160 + 2162 = 2. The row after it, where r = 0, is also informative: 23*240 - 120*46 = 0, and 23 = 46/2 and 120 = 240/2. Those are b/g and a/g, which is the step between consecutive solutions, the starting point for the family of all solutions covered on the Diophantine page.

Notice the signs: s and t alternate in sign from row to row and grow in magnitude. That pattern is what keeps fixed-width implementations safe, as shown below.

Iterative code and the sign of the gcd

The iterative form keeps only the last two rows, so it uses constant extra space.

def ext_gcd(a, b):
    """Return (g, x, y) with a*x + b*y == g == gcd(a, b) >= 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
    if old_r < 0:                      # floor division can leave -gcd
        old_r, old_s, old_t = -old_r, -old_s, -old_t
    return old_r, old_s, old_t

print(ext_gcd(240, 46))   # (2, -9, 47)
print(ext_gcd(46, 240))   # (2, 47, -9)
print(ext_gcd(5, 0))      # (5, 1, 0)
print(ext_gcd(0, 0))      # (0, 1, 0): gcd(0, 0) is 0 by convention

The sign fix at the end is not cosmetic. With Python's floor division the remainders take the sign of the divisor, so without it the loop returns (-2, 9, 47) for inputs 240 and -46: a valid identity, since 240*9 + (-46)*47 = -2, but a negative gcd that breaks every caller testing g == 1. Languages with truncating division, such as C and Java, have the same issue with different signs. The contract you want is a non-negative gcd, so normalise once at the end by flipping all three values.

A recursive version reads more like the math and is fine for small numbers: if b is zero return (a, 1, 0); otherwise take (g, x, y) from (b, a mod b) and return (g, y, x - (a // b) * y). Its recursion depth equals the number of division steps, which is small, as the next section shows, so stack depth is not a concern; the iterative version is preferred mostly because it is easy to translate to fixed-width code.

How many steps, and how big the numbers get

Two questions decide whether this is safe in production code: how many steps it takes, and how large the numbers get.

Steps. The worst case is consecutive Fibonacci numbers, where every quotient is 1 and remainders shrink as slowly as possible. Lame's theorem turns that into a bound: the number of division steps is at most about five times the number of decimal digits of the smaller input, so O(log min(a, b)). Measured: 832040 and 514229, the 30th and 29th Fibonacci numbers, take 28 steps. Typical inputs take far fewer; 10^18 and 10^18 - 1 take 2. Each step is one division, so on machine words the algorithm is effectively instant; for multi-thousand-bit numbers, as in RSA, the division cost dominates and libraries use faster variants.

Sizes. For positive a and b, not equal and neither dividing the other, the returned coefficients satisfy |x| <= b/g and |y| <= a/g. A check on 100,000 random pairs below 10^9 with the code above confirms it, and the trace shows the final zero row reaching exactly b/g and a/g. Better still, because consecutive coefficients alternate in sign, each update s_new = s_old - q*s is really |s_new| = |s_old| + q*|s|: magnitudes only add, and the largest one is the final b/g. No intermediate value exceeds the inputs, so if a and b fit in a signed 64-bit integer, every s and t does too. That is the reason the extended algorithm is safe in C or Java with long, as long as you do not compute the product a*x itself in that width.

Modular inverses, one at a time and in batches

The main use is the modular inverse: the number x with a*x = 1 (mod m). Reduce Bezout's identity mod m: if a*x + m*y = 1 then a*x leaves remainder 1. So an inverse exists exactly when gcd(a, m) = 1, and the extended algorithm finds it.

def mod_inverse(a, m):
    g, x, _ = ext_gcd(a % m, m)
    if g != 1:
        raise ValueError(f"{a} has no inverse mod {m}")
    return x % m                      # bring into [0, m)

print(mod_inverse(17, 3120), pow(17, -1, 3120))   # 2753 2753

def batch_inverse(xs, m):
    """Inverses of every x in xs with ONE extended-Euclid call."""
    prefix = [1] * (len(xs) + 1)
    for i, v in enumerate(xs):
        prefix[i + 1] = prefix[i] * v % m
    acc = mod_inverse(prefix[-1], m)  # inverse of the whole product
    out = [0] * len(xs)
    for i in range(len(xs) - 1, -1, -1):
        out[i] = acc * prefix[i] % m  # strip everything except xs[i]
        acc = acc * xs[i] % m
    return out

M = 10**9 + 7
print(batch_inverse([2, 3, 5, 7, 11], M))
# [500000004, 333333336, 400000003, 142857144, 818181824]

The worked number is the classic textbook RSA example: with p = 61 and q = 53, phi = 3120, public exponent e = 17, and the private exponent d is the inverse of 17 mod 3120, which is 2753. Check: 17*2753 = 46801 = 15*3120 + 1. Python 3.8 and later computes the same thing with pow(a, -1, m), and Java has BigInteger.modInverse; use them in real code and keep your own version for fixed-width languages and for understanding.

Batch inversion (often called Montgomery's trick) matters when you need many inverses mod the same m, for example normalising thousands of points in elliptic-curve code or computing inverse factorials. Multiply everything together, invert the product once, and peel inverses off with two passes of multiplications. One extended-Euclid call plus about 3n multiplications replaces n calls. If any element is not invertible the single inversion fails, so handle zeros before batching.

Two cautions. When m is prime you can also invert with Fermat's little theorem, a^(m-2) mod m, using modular exponentiation, which is simpler but slower than Euclid on machine words. And in cryptography the textbook algorithm leaks timing, because its number of steps depends on the secret input; production crypto libraries use constant-time inversion methods instead. Never hand-roll inversion for secret keys.

Merging two congruences

The second classic use is merging two congruences, the building block of the Chinese remainder theorem. Write x = r1 + m1*k and require r1 + m1*k = r2 (mod m2), that is m1*k = r2 - r1 (mod m2). That is solvable only when g = gcd(m1, m2) divides r2 - r1; then dividing through by g makes m1/g invertible mod m2/g, and the Bezout coefficient p from ext_gcd(m1, m2) is exactly that inverse.

def crt(r1, m1, r2, m2):
    """Solve x = r1 (mod m1), x = r2 (mod m2). Return (x, lcm) or None."""
    g, p, _ = ext_gcd(m1, m2)          # m1*p + m2*q = g
    if (r2 - r1) % g:
        return None                     # congruences contradict each other
    lcm = m1 // g * m2
    k = (r2 - r1) // g * p % (m2 // g)
    return (r1 + m1 * k) % lcm, lcm

print(crt(2, 3, 3, 5))   # (8, 15)
print(crt(2, 4, 3, 6))   # None: x even and x odd at once
print(crt(3, 4, 5, 6))   # (11, 12)

The outputs: x = 2 (mod 3) and x = 3 (mod 5) give 8 mod 15. The moduli 4 and 6 share a factor 2, so the remainders 2 and 3 demand that x be even and odd at once: no solution, which the divisibility check catches instead of returning garbage. Remainders 3 and 5 are compatible, giving 11 mod 12, where 12 is the lcm, not the product. Folding this function over a list merges any number of congruences. Watch the size of the lcm: with many moduli it grows past 64 bits quickly, so use arbitrary-precision integers or keep it under a bound you have checked.

Failure modes

• Negative gcd from signed inputs and floor or truncating division. Normalise the sign once at the end.
• Negative inverse: x from the identity can be negative; reduce with x % m, and in C or Java add m when the remainder is negative.
• Overflow in the check, not the algorithm: a*x may not fit in 64 bits even though a and x do. Verify with 128-bit or arbitrary-precision arithmetic.
• Assuming an inverse exists: if gcd(a, m) is not 1, there is none. Raise an error instead of returning a wrong number.
• Reducing a mod m first and forgetting the zero case: if a is a multiple of m, the inverse does not exist and ext_gcd returns g = m.
• Timing leaks when used on secrets. Use a vetted constant-time library.

What to do next

1. Type in the iterative ext_gcd, run the 240 and 46 trace, and check each row against r = s*a + t*b.
2. Run it on inputs with every sign combination and confirm the gcd is never negative.
3. Write mod_inverse and test it against pow(a, -1, m) on random coprime pairs, and on a non-coprime pair to see the error.
4. Port ext_gcd to C, Java or Rust with 64-bit integers and verify on Fibonacci pairs near the top of the range.
5. Implement batch_inverse and use it to compute inverse factorials mod 10^9 + 7.
6. Write crt, fold it over three or more congruences, and include a contradictory pair in your tests.
7. Continue with Linear Diophantine Equations to turn one Bezout solution into all of them.