Extended Euclidean algorithm
The Extended Euclidean algorithm is an efficient computational procedure in number theory that extends the classical Euclidean algorithm to not only determine the greatest common divisor (GCD) of two integers a and b, denoted d = gcd(a, b), but also to find integers x and y such that a x + b y = d.[1][2] This linear combination expresses d as a Bézout combination of a and b, proving Bézout's identity, which states that such coefficients always exist for any pair of integers.[1][2]
The algorithm operates by iteratively applying the division algorithm, tracking quotients and remainders while maintaining expressions for each remainder as a linear combination of the original inputs.[1][2] It begins with r0 = a = 1·a + 0·b and r1 = b = 0·a + 1·b, then computes subsequent remainders ri = ri-2 - qi · ri-1, updating the coefficients accordingly until a remainder of zero is reached, at which point the penultimate remainder is the GCD and its coefficients provide the solution.[1][2] The process runs in O(log min(a, b)) time complexity, making it highly practical for large integers.[3][4]
Originating as an enhancement to Euclid's algorithm from Elements (circa 300 BCE), the extended version explicitly constructs the coefficients implied by the original method, with the theorem formalized as Bézout's identity by Étienne Bézout in 1779, though the core idea predates him.[5][6] Key applications include solving linear Diophantine equations a x + b y = c (feasible when d divides c), computing modular multiplicative inverses when gcd(a, m) = 1 (vital for systems like RSA encryption), and facilitating primality testing and factorization in computational number theory.[1][2][7] These uses underscore its foundational role in modern cryptography, computer science, and algebraic structures like rings.[8][9]
Fundamentals
Definition and Description
The extended Euclidean algorithm builds upon the classical Euclidean algorithm, originally described by the ancient Greek mathematician Euclid in the 3rd century BCE for computing the greatest common divisor (GCD) of two integers. Early extensions of this method to find integer coefficients expressing the GCD as a linear combination trace back to the Indian mathematician Aryabhata around 500 CE, who developed a procedure known as the "pulverizer" (kuttaka) for solving linear Diophantine equations.[10] In the 18th century, French mathematician Étienne Bézout formalized the underlying identity in his 1779 work Théorie générale des équations algébriques, establishing it for both integers and polynomials, though the algorithmic back-substitution process predates this attribution. The algorithm's significance in modern computational number theory emerged in the 20th century, where it became essential for efficient implementations in areas like cryptography and algebraic computation.[11]
Formally, given two integers a and b with a > b \geq 0 and not both zero, the extended Euclidean algorithm computes d = \gcd(a, b) along with integers x and y satisfying Bézout's identity:
ax + by = d.
This equation guarantees that the GCD can always be expressed as an integer linear combination of a and b, reflecting the ideal generated by a and b in the ring of integers.[10] The algorithm proceeds at a high level by iteratively applying the division step of the Euclidean algorithm—replacing (a, b) with (b, r) where r = a - q b for quotient q, until the remainder is zero—while tracking coefficients through back-substitution to reconstruct the linear combination for the final non-zero remainder, which is d.[12]
The Bézout coefficients x and y are not unique; if (x', y') is another pair satisfying the equation, then x' = x + k (b/d) and y' = y - k (a/d) for some integer k, making the solutions unique up to multiples of (b/d, -a/d).[13] This modulo structure arises directly from the equation's homogeneity and the properties of the GCD.[12]
Illustrative Example
To illustrate the extended Euclidean algorithm, consider the problem of computing \gcd(240, 46) and finding integers x and y such that $240x + 46y = \gcd(240, 46). This process follows Bézout's identity, which guarantees the existence of such integers for any pair of integers.[14]
The algorithm proceeds in two phases: first, applying the Euclidean algorithm to find the gcd via successive divisions, tracking quotients and remainders; second, back-substituting to express the gcd as a linear combination. The following table summarizes the forward phase, where each remainder r_i is computed as r_{i} = r_{i-2} - q_i r_{i-1}, with initial values r_0 = 240 and r_1 = 46.
| Step | r_i (remainder) | Quotient q_i | Expression for r_i |
|---|
| 0 | 240 | - | $240 = 1 \cdot 240 + 0 \cdot 46 |
| 1 | 46 | - | $46 = 0 \cdot 240 + 1 \cdot 46 |
| 2 | 10 | 5 | $10 = 240 - 5 \cdot 46 |
| 3 | 6 | 4 | $6 = 46 - 4 \cdot 10 |
| 4 | 4 | 1 | $4 = 10 - 1 \cdot 6 |
| 5 | 2 | 1 | $2 = 6 - 1 \cdot 4 |
| 6 | 0 | 2 | $0 = 4 - 2 \cdot 2 |
Mathematical Proof
The correctness of the extended Euclidean algorithm, which computes integers x and y such that ax + by = \gcd(a, b) for nonnegative integers a \geq b \geq 0, can be established by mathematical induction on b, leveraging Bézout's identity that \gcd(a, b) is the smallest positive linear combination of a and b.[15][16]
Base case: If b = 0, then \gcd(a, 0) = a, and the algorithm returns x = 1, y = 0, satisfying a \cdot 1 + 0 \cdot 0 = a. This holds since any integer divides zero, and the greatest divisor of a and zero is a itself.[15][17]
Inductive step: Assume the algorithm is correct for all pairs with second argument less than b > 0; that is, for the recursive call on \gcd(b, r) where r = a \mod b = a - qb and q = \lfloor a/b \rfloor, it returns integers x' and y' such that b x' + r y' = \gcd(b, r) = d = \gcd(a, b). Substituting r = a - qb yields b x' + (a - qb) y' = d, or a y' + b (x' - q y') = d. Thus, setting x = y' and y = x' - q y' satisfies a x + b y = d, confirming the algorithm's output for the original pair.[15][16]
The algorithm terminates because the underlying Euclidean algorithm does: each remainder strictly decreases (r < b) and remains nonnegative, eventually reaching zero after finitely many steps, at most O(\log b) iterations by properties of the Fibonacci sequence bounding the worst case.[18][15]
The solution (x, y) is not unique; all integer solutions to a x + b y = d are given by x = x_0 + k (b/d) and y = y_0 - k (a/d) for integer k, where (x_0, y_0) is any particular solution and d = \gcd(a, b), since adding multiples of b/d to x and subtracting multiples of a/d from y preserves the equation due to a (b/d) = b (a/d).[19][20]
Algorithmic Implementation
The recursive formulation of the extended Euclidean algorithm computes the greatest common divisor d = \gcd(a, b) of two integers a and b (assuming a \geq b \geq 0), along with the Bézout coefficients x and y satisfying the equation ax + by = d. This version naturally extends the recursive structure of the standard Euclidean algorithm by tracking the coefficients through back-substitution in the recursion.[21]
The function, typically denoted as extended_gcd(a, b), returns a tuple (d, x, y). In the base case, if b = 0, it returns (a, 1, 0), since \gcd(a, 0) = a and a \cdot 1 + 0 \cdot 0 = a.[21]
For the recursive case, the algorithm calls itself on the pair (b, a \mod b) to obtain (d, x_1, y_1) where d = \gcd(b, a \mod b) and b x_1 + (a \mod b) y_1 = d. It then computes the coefficients for the original pair as x = y_1 and y = x_1 - \left\lfloor \frac{a}{b} \right\rfloor y_1, before returning (d, x, y). This adjustment ensures the linear combination holds for a and b, mirroring the division step in the Euclidean algorithm. The following pseudocode illustrates this structure:
function extended_gcd(a, b):
if b == 0:
return (a, 1, 0)
else:
(d, x1, y1) = extended_gcd(b, a % b)
x = y1
y = x1 - (a // b) * y1
return (d, x, y)
function extended_gcd(a, b):
if b == 0:
return (a, 1, 0)
else:
(d, x1, y1) = extended_gcd(b, a % b)
x = y1
y = x1 - (a // b) * y1
return (d, x, y)
[21]
The time complexity of this recursive formulation is O(\log \min(a, b)) arithmetic operations, matching that of the standard Euclidean algorithm, as the number of recursive calls equals the number of division steps until the base case is reached, and each call performs constant work beyond the recursion. The recursion depth is also O(\log \min(a, b)), since the arguments decrease rapidly.[22]
This recursive approach offers advantages in conceptual clarity, as its structure directly parallels the inductive proof of the Euclidean algorithm's correctness, making it particularly suitable for understanding and verifying the computation of Bézout coefficients through mathematical induction on the input size.[23]
Iterative Pseudocode
The iterative formulation of the extended Euclidean algorithm computes the greatest common divisor d = \gcd(a, b) along with Bézout coefficients x and y such that ax + by = d, using a loop to simulate the divisions and back-substitutions without recursion. This approach is particularly useful in programming environments where recursion depth limits could be exceeded for large inputs.[24]
The following language-agnostic pseudocode illustrates the core structure, initializing coefficients for the initial remainders and updating them in each iteration while shifting the values of a and b:
function extended_gcd(a, b):
if b == 0:
return a, 1, 0
prevx, x = 1, 0
prevy, y = 0, 1
while b != 0:
q = a // b
a, b = b, a % b
x, prevx = prevx - q * x, x
y, prevy = prevy - q * y, y
return a, prevx, prevy // returns d, x, y where original_a * x + original_b * y = d
function extended_gcd(a, b):
if b == 0:
return a, 1, 0
prevx, x = 1, 0
prevy, y = 0, 1
while b != 0:
q = a // b
a, b = b, a % b
x, prevx = prevx - q * x, x
y, prevy = prevy - q * y, y
return a, prevx, prevy // returns d, x, y where original_a * x + original_b * y = d
In this implementation, the variables track the coefficients for the current and previous remainders, ensuring the invariant holds after each update.[24]
An array-based variant builds sequences for the remainders r, coefficients s, and coefficients t iteratively, akin to constructing a table of values. Initialize r_0 = a, r_1 = b, s_0 = 1, s_1 = 0, t_0 = 0, t_1 = 1. Then, for each subsequent index i \geq 2, compute the quotient q_i = \lfloor r_{i-2} / r_{i-1} \rfloor, and update r_i = r_{i-2} - q_i r_{i-1}, s_i = s_{i-2} - q_i s_{i-1}, t_i = t_{i-2} - q_i t_{i-1}. Continue until r_i = 0; the GCD is the last non-zero remainder r_{i-1}, with corresponding s_{i-1} and t_{i-1}. This method explicitly stores all intermediate coefficients, facilitating verification or further computations.
To handle edge cases, if a = 0 and b \neq 0, return d = |b|, x = 0, y = \operatorname{sign}(b); if b = 0, return d = |a|, x = \operatorname{sign}(a), y = 0. For negative inputs, apply the algorithm to |a| and |b| and adjust the signs of the coefficients based on the original signs to ensure a x + b y = d > 0.[25]
The iterative versions achieve the same time complexity as the underlying Euclidean algorithm, O(\log \min(|a|, |b|)), due to the logarithmic number of divisions, while avoiding recursion overhead and potential stack issues.[26]
Generalizations
Polynomial Version
The extended Euclidean algorithm generalizes to univariate polynomials over a field K, such as the rational numbers or real numbers, in the polynomial ring K. Given two polynomials f, g \in K, not both zero, the algorithm computes a monic greatest common divisor d \in K and polynomials s, t \in K satisfying Bézout's identity f s + g t = d. This adaptation relies on the fact that K forms a Euclidean domain, where the degree function serves as the Euclidean norm, enabling division with remainder analogous to the integer case.[27][28]
The algorithm proceeds recursively, mirroring the structure for integers but using polynomial division. Without loss of generality, assume \deg(f) \geq \deg(g). Divide f by g to obtain the quotient q and remainder r such that f = q g + r with \deg(r) < \deg(g). Recursively apply the algorithm to g and r, yielding polynomials s' and t' where g s' + r t' = d. Back-substituting the expression for r gives f s' + g (t' - q s') = d, so set s = s' and t = t' - q s'. The base cases occur when g = 0, in which d = f (normalized to be monic) with s = 1 and t = 0, or when g is a nonzero constant, in which d = 1 (the monic unit) and the coefficients are scaled accordingly by the leading coefficient of g. Normalization ensures d is monic by dividing by its leading coefficient throughout.[29][28]
For illustration, consider f(x) = x^2 + 1 and g(x) = x + 1 over the real numbers \mathbb{R}. Dividing f by g yields quotient x - 1 and constant remainder $2. The recursion on g and $2 produces the monic gcd d = 1, and back-substitution provides explicit s and t satisfying the identity, though the full computation involves scaling by the leading coefficient of the remainder.[29]
The computational complexity of this algorithm is O(\deg(f) \cdot \deg(g)) arithmetic operations in the field K, assuming standard polynomial division algorithms. With faster multiplication techniques, such as those based on fast Fourier transforms, the complexity can improve, but the basic version aligns with the quadratic scaling in degrees.[30]
Extension to Multiple Arguments
The extended Euclidean algorithm can be generalized to compute the greatest common divisor d of n > 2 integers a_1, a_2, \dots, a_n and express d as an integer linear combination d = \sum_{i=1}^n x_i a_i, where the coefficients x_i are integers (not necessarily unique).[31] This generalization relies on Bézout's identity extended to multiple integers, which holds by induction: the base case for two integers is given by the standard extended Euclidean algorithm, and for additional integers, the result follows from expressing the GCD of the first n-1 as a linear combination and then incorporating the nth via another application of the algorithm.[31]
The method proceeds iteratively by pairwise application. Begin by applying the extended Euclidean algorithm to the first two integers a_1 and a_2 to obtain d_2 = \gcd(a_1, a_2) = x_1 a_1 + x_2 a_2. Then, compute d_3 = \gcd(d_2, a_3) = y_1 d_2 + y_3 a_3 using the extended Euclidean algorithm again, and substitute the expression for d_2 to yield d_3 = (y_1 x_1) a_1 + (y_1 x_2) a_2 + y_3 a_3. Continue this process sequentially for each subsequent a_i, updating the coefficients for prior terms by multiplication with the new intermediate coefficient and adding the new term's coefficient. This yields the full linear combination after n-1 applications.[32]
For example, consider the integers 48, 18, and 30. First, apply the extended Euclidean algorithm to 48 and 18:
\begin{align*}
48 &= 2 \cdot 18 + 12, \\
18 &= 1 \cdot 12 + 6, \\
12 &= 2 \cdot 6 + 0.
\end{align*}
Back-substituting gives \gcd(48, 18) = 6 = -1 \cdot 48 + 3 \cdot 18. Now incorporate 30 by computing \gcd(6, 30):
\begin{align*}
30 &= 5 \cdot 6 + 0.
\end{align*}
Thus, $6 = 1 \cdot 6 + 0 \cdot 30. Substituting the prior expression yields $6 = -1 \cdot 48 + 3 \cdot 18 + 0 \cdot 30. The coefficients are updated accordingly, confirming d = 6.[32]
In this iterative approach, the coefficients x_i can grow rapidly—potentially exponentially in n—due to successive multiplications during substitution, especially when intermediate quotients are large, although minimal coefficients satisfying the equation are bounded by the sum of the |a_i|.[33] This growth makes the method less efficient for large n, as the bit length of coefficients may increase significantly. If only the GCD is needed without coefficients, an alternative is to use Stein's algorithm (also known as the binary GCD algorithm) iteratively in a pairwise fashion, which avoids divisions and uses bit shifts and subtractions for faster computation on large integers, though it does not directly provide Bézout coefficients.[34]
Applications
Rational Number Simplification
To simplify a rational number expressed as a fraction \frac{p}{q} where p and q are integers and q \neq 0, the extended Euclidean algorithm (EEA) is applied to compute the greatest common divisor d = \gcd(|p|, |q|).[35] The simplified form is then \frac{p/d}{q/d}, with signs adjusted to ensure the denominator is positive (e.g., if q < 0, multiply numerator and denominator by -1).[25] This process reduces the fraction to its lowest terms by dividing both numerator and denominator by their common divisor d.[36]
The EEA is particularly efficient for this task because it computes the GCD in O(\log \min(|p|, |q|)) steps using repeated division, outperforming naive factorization methods for large integers.[35] Although the EEA also outputs Bézout coefficients expressing d as an integer linear combination of |p| and |q|, these are incidental for simplification and not required.[25]
For example, consider the fraction \frac{240}{46}. Applying the EEA yields \gcd(240, 46) = 2, so the simplified form is \frac{240/2}{46/2} = \frac{120}{23}.[1]
Edge cases must be handled carefully: if q = 0, the fraction is undefined; if p = 0, it simplifies to \frac{0}{1}.[25]
The steps of the EEA on p and q directly mirror the continued fraction expansion of \frac{p}{q}, where the quotients from successive divisions form the continued fraction coefficients, terminating finitely for rationals.[37]
Historically, the EEA's precursor—the Euclidean algorithm for GCD—has been used since Euclid's Elements (Book VII, c. 300 BCE) in early arithmetic for handling exact fractions in rational computations.[35]
Modular Multiplicative Inverses
In modular arithmetic, an integer a has a multiplicative inverse modulo m if and only if \gcd(a, m) = 1, meaning there exists an integer x such that a x \equiv 1 \pmod{m}.[38] The extended Euclidean algorithm computes this inverse by finding integers x and y satisfying Bézout's identity a x + m y = 1, from which x \mod m (adjusted to a positive representative between 0 and m-1) serves as the inverse.[38]
To compute the inverse, apply the extended Euclidean algorithm to a and m:
- Perform the standard Euclidean algorithm to find \gcd(a, m); if it is not 1, no inverse exists.
- Use the extended version to back-substitute and express 1 as a linear combination a x + m y = 1.
- The coefficient x is the inverse, taken modulo m.
For example, to find the inverse of 5 modulo 17:
\begin{align*}
17 &= 3 \cdot 5 + 2, \\
5 &= 2 \cdot 2 + 1, \\
2 &= 2 \cdot 1 + 0.
\end{align*}
Back-substituting:
$1 = 5 - 2 \cdot 2, \quad 2 = 17 - 3 \cdot 5,
$1 = 5 - 2 \cdot (17 - 3 \cdot 5) = 7 \cdot 5 - 2 \cdot 17.
Thus, $5 \cdot 7 \equiv 1 \pmod{17}, so the inverse is 7.[38]
If \gcd(a, m) = d > 1, no multiplicative inverse exists for a modulo m when solving a x \equiv 1 \pmod{m}, as d does not divide 1; more generally, the linear congruence a x \equiv b \pmod{m} is solvable if and only if d divides b.[38]
The extended Euclidean algorithm's efficiency in computing these inverses is crucial in cryptography, such as RSA decryption, where the private exponent d is the modular inverse of the public exponent e modulo \phi(n), computed via the algorithm during key generation.[39] Similarly, in Diffie-Hellman key exchange variants like implicitly-certified public keys, inverses modulo p-1 are required to derive private keys from shared parameters.[40]
Wilson's theorem connects to modular inverses by stating that for a prime p, (p-1)! \equiv -1 \pmod{p}, reflecting the structure of the multiplicative group modulo p where every non-zero element has a unique inverse.[41]
Inverses in Algebraic Extensions
In algebraic extensions of a field K, such as simple extensions K[\alpha] where \alpha is a root of an irreducible polynomial m(x) \in K, elements can be represented as polynomials in \alpha of degree less than the degree of m(x). The extended Euclidean algorithm applied to polynomials over K enables the computation of multiplicative inverses for elements that are coprime to m(x).[42]
To invert an element \beta \in K[\alpha], represent \beta as a polynomial \beta(x) \in K with \deg(\beta(x)) < \deg(m(x)). Apply the polynomial extended Euclidean algorithm to \beta(x) and m(x), which computes the greatest common divisor and Bézout coefficients s(x), t(x) \in K such that s(x) \beta(x) + t(x) m(x) = \gcd(\beta(x), m(x)). Since m(x) is irreducible and \beta is invertible if \gcd(\beta(x), m(x)) = 1 (a constant in K), the coefficients yield s(x) \beta(x) \equiv 1 \pmod{m(x)}, so the inverse is s(\alpha) reduced modulo m(x). This process leverages the division algorithm in K to generate the remainder sequence until the gcd is reached, then back-substitutes to express the gcd.[43][42]
For example, consider the quadratic extension \mathbb{Q}(\sqrt{2}), where \alpha = \sqrt{2} satisfies the minimal polynomial m(x) = x^2 - 2. To invert \beta = 1 + \sqrt{2}, represent it as \beta(x) = x + 1. Applying the extended Euclidean algorithm over \mathbb{Q}:
\begin{align*}
x^2 - 2 &= (x - 1)(x + 1) + (-1).
\end{align*}
The next step divides x + 1 by -1, yielding quotient -(x + 1) and remainder 0, so the GCD is -1 (a unit). Back-substituting:
-1 = x^2 - 2 - (x - 1)(x + 1),
$1 = -(x^2 - 2) + (x - 1)(x + 1).
Thus, s(x) = x - 1, so the inverse is \sqrt{2} - 1. Verification: (1 + \sqrt{2})(\sqrt{2} - 1) = 2 - 1 = 1.
This method generalizes to finite-degree extensions K(\alpha_1, \dots, \alpha_n), where inverses reduce to applications of the polynomial extended Euclidean algorithm in towers of simple extensions, provided the base field K supports exact arithmetic.[43]
In computer algebra systems like SageMath, this technique is implemented for efficient inversion in quadratic and higher-degree number fields, supporting symbolic computations in algebraic geometry and number theory. Analogs appear in elliptic curve cryptography, where inverses in finite field extensions (e.g., \mathbb{F}_{p^k}) are computed via the extended Euclidean algorithm for point arithmetic over composite fields.[44][45]
The approach requires coefficients in a field K (e.g., \mathbb{Q} or \mathbb{F}_p) for exact division; it does not directly apply to composite extensions without primitive element decomposition or handling non-simple structures.[42]