Extended Euclidean Algorithm: Bezout's Identity & Modular Inverses

Have you ever needed to find a number that, when multiplied by another, gives a remainder of 1? That's a modular inverse, and it's crucial in cryptography, hash functions, and more. The Extended Euclidean Algorithm is the tool that makes this possible. It's a classic number theory algorithm that not only computes the greatest common divisor (GCD) of two integers but also finds the coefficients of Bezout's identity: integers x and y such that ax + by = gcd(a, b). This might sound abstract, but it has direct applications: from solving linear Diophantine equations to computing modular inverses for RSA encryption. In this tutorial, you'll learn how the algorithm works, see it in action with Python code, and understand how to use it in production systems. We'll also cover common pitfalls and debugging strategies. By the end, you'll be able to implement the Extended Euclidean Algorithm confidently and apply it to real-world problems.

Understanding the Euclidean Algorithm

The Euclidean algorithm is a method for finding the greatest common divisor (GCD) of two integers. It is based on the principle that gcd(a, b) = gcd(b, a mod b). The algorithm repeatedly replaces (a, b) with (b, a % b) until b becomes 0, at which point a is the GCD. For example, to find gcd(30, 12): 30 % 12 = 6, so next pair is (12, 6) 12 % 6 = 0, so next pair is (6, 0) => GCD = 6. This is simple and efficient with O(log min(a, b)) time complexity.

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

print(gcd(30, 12))  # Output: 6
print(gcd(17, 5))   # Output: 1

Bezout's Identity: The Core Concept

Bezout's Identity states that for any integers a and b, there exist integers x and y such that ax + by = gcd(a, b). These coefficients are not unique; the Extended Euclidean Algorithm finds one such pair. For example, for a=30 and b=12, gcd=6, and we have 301 + 12(-2) = 6. This identity is fundamental in solving linear Diophantine equations and computing modular inverses. The coefficients can be found by backtracking through the Euclidean algorithm steps.

# Manual backtracking for gcd(30,12):
# Step1: 30 = 2*12 + 6  => 6 = 30 - 2*12
# Step2: 12 = 2*6 + 0   => stop
# So 6 = 30*1 + 12*(-2) => x=1, y=-2
print(f"30*1 + 12*(-2) = {30*1 + 12*(-2)}")  # Output: 6

Implementing the Extended Euclidean Algorithm

The Extended Euclidean Algorithm computes the GCD along with the coefficients x and y. It can be implemented recursively or iteratively. The recursive version follows the recurrence: if b == 0, return (a, 1, 0). Otherwise, recursively compute (g, x1, y1) for (b, a % b), then return (g, y1, x1 - (a // b) * y1). The iterative version uses a loop to update coefficients. Both have O(log min(a, b)) time complexity. Here's a Python implementation:

def extended_gcd(a, b):
    if b == 0:
        return (a, 1, 0)
    else:
        g, x1, y1 = extended_gcd(b, a % b)
        x = y1
        y = x1 - (a // b) * y1
        return (g, x, y)

print(extended_gcd(30, 12))  # Output: (6, 1, -2)
print(extended_gcd(17, 5))   # Output: (1, -2, 7)

Computing Modular Inverses

A modular inverse of a modulo m is an integer x such that ax ≡ 1 (mod m). This exists only if gcd(a, m) = 1. Using the Extended Euclidean Algorithm, we find x and y such that ax + my = 1. Then x mod m is the modular inverse. For example, to find the inverse of 17 modulo 43: extended_gcd(17, 43) returns (1, -5, 2) because 17(-5) + 43*2 = 1. So the inverse is -5 mod 43 = 38. In Python, we can compute it as:

def mod_inverse(a, m):
    g, x, y = extended_gcd(a, m)
    if g != 1:
        raise ValueError('Inverse does not exist')
    else:
        return x % m

print(mod_inverse(17, 43))  # Output: 38
print(mod_inverse(30, 12))  # Raises ValueError

Solving Linear Diophantine Equations

A linear Diophantine equation is of the form ax + by = c. It has integer solutions if and only if gcd(a, b) divides c. The Extended Euclidean Algorithm provides a particular solution (x0, y0) for ax + by = g. Then all solutions are given by x = x0(c/g) + (b/g)t, y = y0(c/g) - (a/g)t for integer t. For example, solve 30x + 12y = 18. Since gcd=6 divides 18, solutions exist. From extended_gcd(30,12) we have (6,1,-2). Multiply by 3: x0=3, y0=-6. General solution: x=3 + 2t, y=-6 -5t.

def solve_diophantine(a, b, c):
    g, x0, y0 = extended_gcd(a, b)
    if c % g != 0:
        return None  # No solution
    x0 *= c // g
    y0 *= c // g
    return (x0, y0, b // g, a // g)  # returns particular solution and step sizes

sol = solve_diophantine(30, 12, 18)
if sol:
    x0, y0, step_x, step_y = sol
    print(f"Particular: x={x0}, y={y0}")
    print(f"General: x = {x0} + {step_x}*t, y = {y0} - {step_y}*t")
else:
    print("No solution")

Iterative Implementation for Production

The recursive implementation is elegant but may cause stack overflow for large inputs. An iterative version is more robust for production. It uses a loop to update coefficients similarly to the recursive version. Here's an iterative implementation:

def extended_gcd_iterative(a, b):
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    while r != 0:
        quotient = old_r // r
        old_r, r = r, old_r - quotient * r
        old_s, s = s, old_s - quotient * s
        old_t, t = t, old_t - quotient * t
    return (old_r, old_s, old_t)

print(extended_gcd_iterative(30, 12))  # Output: (6, 1, -2)
print(extended_gcd_iterative(17, 5))   # Output: (1, -2, 7)