Extended Euclidean algorithm

The extended Euclidean algorithm is an extension to the Euclidean algorithm. Besides finding the greatest common divisor of integers a and b, as the Euclidean algorithm does, it also finds integers x and y (one of which is typically negative) that satisfy Bézout's identity

$ax + by = \gcd(a, b).$

The extended Euclidean algorithm is particularly useful when a and b are coprime, since x is the modular multiplicative inverse of a modulo b, and y is the modular multiplicative inverse of b modulo a. This has value in a key calculation of the RSA public-key encryption algorithm; finding an integer decryption exponent d that is the modular multiplicative inverse of the chosen encryption exponent e modulo φ(n), where e and φ(n) are coprime.

Informal formulation of the algorithm

Dividend	Divisor	Quotient	Remainder
120	23	5	5
23	5	4	3
5	3	1	2
3	2	1	1
2	1	2	0

To illustrate the extension of Euclid's algorithm, consider the computation of gcd(120, 23), which is shown on the table on the left. Notice that the quotient in each division is recorded as well alongside the remainder. In this example, the divisor in the last line (which is equal to 1) indicates that the gcd is 1; that is, 120 and 23 are coprime (also called relatively prime).

Dividend	Divisor	Quotient	Remainder
120	25	4	20
25	20	1	5
20	5	4	0

The table on the right shows the computation for gcd(120, 25). In this case, the gcd is 5, the divisor in the last line.

There are two methods to proceed, both using an integer division algorithm as a subprocedure, which will be discussed separately.

Iterative method

This method computes expressions of the form $r_i = a x_i + b y_i$ for the remainder in each step $i$ of the Euclidean algorithm. Each successive number $r_i$ can be written as the remainder of the division of the previous two such numbers, which remainder can be expressed using the whole quotient q_i of that division as follows:

$r_i = r_{i-2} - q_i r_{i-1}.\,$

By substitution, this gives:

$r_i = (ax_{i-2} + by_{i-2}) - q_i (ax_{i-1} + by_{i-1}),\,$ which can be written

$r_i = a(x_{i-2} - q_i x_{i-1}) + b(y_{i-2} - q_i y_{i-1}).\,$

The first two values are the initial arguments to the algorithm:

$r_1 = a = a\times1 + b\times0\,$

$r_2 = b = a\times0 + b\times1.\,$

So the coefficients start out as x₁ = 1, y₁ = 0, x₂ = 0, and y₂ = 1, and the others are given by

$x_i = x_{i-2} - q_i x_{i-1},\,$

$y_i = y_{i-2} - q_i y_{i-1}.\,$

The expression for the last non-zero remainder gives the desired results since this method computes every remainder in terms of a and b, as desired.

Example: Compute the GCD of 120 and 23.

The computation proceeds as follows:

Step	Quotient	Remainder	Substitute	Combine terms
1		120		120 = 120 × 1 + 23 × 0
2		23		23 = 120 × 0 + 23 × 1
3	[120 ÷ 23] = 5	120 − 23 × 5 = 5	5 = (120 × 1 + 23 × 0) − (120 × 0 + 23 × 1) × 5	5 = 120 × 1 + 23 × −5
4	[23 ÷ 5] = 4	23 − 5 × 4 = 3	3 = (120 × 0 + 23 × 1) − (120 × 1 + 23 × −5) × 4	3 = 120 × −4 + 23 × 21
5	[5 ÷ 3] = 1	5 − 3 × 1 = 2	2 = (120 × 1 + 23 × −5) − (120 × −4 + 23 × 21) × 1	2 = 120 × 5 + 23 × −26
6	[3 ÷ 2] = 1	3 − 2 × 1 = 1	1 = (120 × −4 + 23 × 21) − (120 × 5 + 23 × −26) × 1	1 = 120 × −9 + 23 × 47
7	[2 ÷ 1] = 2	2 − 1 × 2 = 0	End of algorithm

The last line reads 1 = 120 × −9 + 23 × 47, which is the required solution: x = −9 and y = 47.

Since the GCD turns out to be 1, this also means that 47 gives the multiplicative inverse of 23 in modulo 120, and also that in modulo 23, −9 (or 14, which represents the same element) gives the multiplicative inverse of 120 (or equivalently of 5).

47 × 23 ≡ 1 (mod 120) and also −9 × 120 ≡14 × 5 ≡ 1 (mod 23).

For an integer a to be invertible modulo b, it is necessary that a and b be coprime, so if a GCD greater than 1 is found, one may conclude that such a modular inverse does not exist.

Note that the equations defining the values x_i are independent of those defining the values y_i, and that the Bézout identity found at the end

$\mathit{GCD} = r_k = a x_k + b y_k\,$

allows deducing y_k when x_k is known. One can therefore omit the values y_i from the computation (although they can be useful to check for computational errors). This is particularly relevant for applications, such as computing modular inverses, that require the value of only one Bézout coefficient.

Recursive method

This method attempts to solve the original equation directly, by reducing the dividend and divisor gradually, from the first line to the last line, which can then be substituted with trivial value and work backward to obtain the solution.

Consider the original equation:

120	x	+	23	y	=	1
(5×23 + 5)	x	+	23	y	=	1
23	(5x + y)	+	5	x	=	1
...
1	xk	+	0	yk	=	1

Notice that the equation remains unchanged after decomposing the original dividend in terms of the divisor plus a remainder, and then regrouping terms. If we have a solution to the equation in the second line, then we can work backward to find x and y as required. Although we don't have the solution yet to the second line, notice how the magnitude of the terms decreased (120 and 23 to 23 and 5). Hence, if we keep applying this, eventually we'll reach the last line, which obviously has (1, 0) as a trivial solution. Then we can work backward and gradually find out x and y.

Dividend	=	Quotient	x	Divisor	+	Remainder
120	=	5	x	23	+	5
23	=	4	x	5	+	3
...

For the purpose of explaining this method, the full working will not be shown. Instead some of the repeating steps will be described to demonstrate the principle behind this method.

Start by rewriting each line from the first table using the division algorithm, focusing on the dividend this time (because we'll be substituting the dividend).

120	x0	+	23	y0	=	1
(5×23+5)	x0	+	23	y0	=	1
23	(5x0+y0)	+	5	x0	=	1
23	x1	+	5	y1	=	1
(4×5+3)	x1	+	5	y1	=	1
5	(4x1+y1)	+	3	x1	=	1
5	x2	+	3	y2	=	1

Assume that we were given x2 = 2 and y2 = −3 already, which is indeed a valid solution. x1 = y2 = −3 Solve 4x1+y1=x2 by substituting x1 = −3, which gives y1 = 2 − 4(−3) = 14 x0 = y1 = 14 Solve 5x0 + y0 = x1 by substituting x0=14, so y0 = −3 − 5(14)= −73

Table method

The table method, also known as "The Magic Box", is probably the simplest method to carry out with a pencil and paper. It is similar to the iterative method, although it does not directly require algebra to use. Let $\operatorname{rem}(a, b)$ denote the remainder $r$ in the division $a = bq + r$. The main idea is to think of the equation chain

$\gcd(a, b) = \gcd(b, \operatorname{rem}(a, b)) = \dots = \gcd(c, 1)$

as a sequence of divisors $a, b, \operatorname{rem}(a, b), \dots, 1$. In the running example we have the sequence 120, 23, 5, 3, 2, 1. Any element in this chain can be written as a linear combination of the original $a$ and $b$, most notably, the last element, $\gcd(a, b)$, can be written in this way. The table method involves keeping a table of each divisor, written as a linear combination. The algorithm starts with the table as follows:

x	y	d
1	0	120
0	1	23

The elements in the $d$ column of the table will be the divisors in the sequence. Each $d_i$ can be represented as the linear combination

$d_i = x_i \cdot a + y_i \cdot b.$

The $x$ and $y$ values are obvious for the first two rows of the table, which represent $a$ and $b$ themselves. To compute $d_i$ for any $i > 2$, notice that

$d_i = \operatorname{rem}(d_{i-2}, d_{i-1}).$

Suppose $d_i = d_{i-2} - k_{i-1} \cdot d_{i-1}$. Then it must be that

$x_i = x_{i-2} - k_{i-1} \cdot x_{i-1}$ and

$y_i = y_{i-2} - k_{i-1} \cdot y_{i-1}.$

This is easy to verify algebraically with a simple substitution.

Actually carrying out the table method though is simpler than the above equations would indicate. To find the third row of the table in the example, just notice that 120 divided by 23 goes 5 times plus a remainder. This gives us k, the multiplying factor for this row. Now, each value in the table is the value two rows above it, minus k times the value immediately above it. This correctly leads to

$x_3 = 1 - 5 \cdot 0 = 1$,

$y_3 = 0 - 5 \cdot 1 = -5$, and

$d_3 = 120 - 5 \cdot 23 = 5.$

After repeating this method to find each line of the table (note that the remainder written in the table and the multiplying factor are two different numbers!), the final values for $x$ and $y$ will solve $ax + by = \gcd(a, b)\,$:

x	y	d	k
1	0	120
0	1	23	5
1	-5	5	4
-4	21	3	1
5	-26	2	1

This method is simple, requiring only the repeated application of one rule, and leaves the answer in the final row of the table with no backtracking. Note also that if you are intending to find a modular inverse for a modulo b and you end up with x being negative, you then need to add the modulus b to x in order to bring it into the range $0< x < b$. This does not affect the validity of the solution as we have $a(x+b)+b(y-a)=ax+by$.

Formal description of the algorithm

Iterative method

By routine algebra of expanding and grouping like terms (refer to last section), the following algorithm for iterative method is obtained:

1. Apply Euclidean algorithm, and let q_n (n starts from 1) be a finite list of quotients in the division.
2. Initialize x₀, x₁ as 1, 0, and y₀, y₁ as 0,1 respectively.
   1. Then for each i so long as q_i is defined,
   2. Compute x_i+1 = x_i−1 − q_ix_i
   3. Compute y_i+1 = y_i−1 − q_iy_i
   4. Repeat the above after incrementing i by 1.
3. The answers are the second-to-last of x_n and y_n.

Pseudocode for this method is shown below:

function extended_gcd(a, b)
    x := 0    lastx := 1
    y := 1    lasty := 0
    while b ≠ 0
        quotient := a div b
        (a, b) := (b, a mod b)
        (x, lastx) := (lastx - quotient*x, x)
        (y, lasty) := (lasty - quotient*y, y)       
    return (lastx, lasty)

Extra steps are needed when working with negative a and/or b.

Recursive method

Solving the general case of the equation in the last corresponding section, the following algorithm results. Given any pair of nonnegative integers a, b, it finds integers values x, y such that ax + by is nonnegative and divides both a and b, which implies that it is the greatest common divisor of a and b.

1. If b = 0, the algorithm ends, returning the solution x = 1, y = 0.
2. Otherwise:
   • Determine the quotient q and remainder r of dividing a by b using an integer division algorithm.
   • Then recursively find coefficients s, t such that bs + rt divides both b and r.
   • Finally the algorithm returns the solution x = t, and y = s − qt.

This can be written in pseudocode as follows:

function extended_gcd(a, b)
    if b == 0
        return (1, 0)
    else
        (q, r) := divide (a, b)
        (s, t) := extended_gcd(b, r)
        return (t, s - q * t)

This assumes a "divide" procedure exists that returns a (quotient,remainder) pair (one could alternatively put q := a div b, and then r = a − b * q).

Proof of correctness

Let a and b be nonnegative integers. We wish to prove that the algorithm terminates, and returns (x, y) such that ax + by divides both a and b. We proceed by induction on b.

• If b = 0, then the algorithm terminates immediately with x = 1 and y = 0, so ax + by = a; this is clearly nonnegative and divides both a and 0 (as 0a = 0).
• Otherwise let q, r be obtained from the division of a by b as in the algorithm description, so a = bq + r and r < b by the properties of Euclidean division.
  • The inequality r < b means we can apply the induction hypothesis for (b,r) in the place of (a,b), and this guarantees that the recursive application terminates.
  • Let (s,t) be the values it returns; by induction we know that bs + rt is nonnegative and divides both b and r.
  • Now the algorithm returns x = t and y = s − qt. We have ax + by = at + b(s − qt) = bs + (a − bq)t = bs + rt, which is (still) nonnegative and divides both b and r. The same number therefore also divides r + bq = a, which terminates the proof.

The proof shows that the recursive algorithm can be adapted to deal with arbitrary integers a, b with just a slight modification: in the terminating case b = 0 it should test the sign of a, and if it is negative return x = −1 instead of x = 1, to ensure that ax + by is always nonnegative. This does assume that the selected division algorithm works if a and/or b are negative, and ensures that |r| < |b| in all cases (the usual specification of Euclidean division in fact requires r to be nonnegative, but this is not necessary in the proof, which is now by recurrence on |b|).

Computing a multiplicative inverse in a finite field

The extended Euclidean algorithm can also be used to calculate the modular multiplicative inverse in a finite field.

Pseudocode

Given the irreducible polynomial f(x) used to define the finite field, and the element a(x) whose inverse is desired, then a form of the algorithm suitable for determining the inverse is given by the following. NOTE: remainder() and quotient() are functions different from the arrays remainder[ ] and quotient[ ]. remainder() refers to the remainder when two numbers are divided, and quotient() refers to the integer quotient when two numbers are divided. Division (with remainder) must be calculated in terms of polynomial arithmetic and not as with "normal" numbers.

pseudocode:

remainder[1] := f(x)
remainder[2] := a(x)
auxiliary[1] := 0
auxiliary[2] := 1
i := 2
while remainder[i] > 1
    i := i + 1
    remainder[i] := remainder(remainder[i-2] / remainder[i-1])
    quotient[i] := quotient(remainder[i-2] / remainder[i-1])
    auxiliary[i] := -quotient[i] * auxiliary[i-1] + auxiliary[i-2]
inverse := auxiliary[i]

Note: In some finite fields, for example GF(2ⁿ), the operations of addition (+) and subtraction (−) are the same. Consequently, some algorithms intended for those fields will not show a minus sign before -quotient[i].

Example

For example, if the polynomial used to define the finite field GF(2⁸) is f(x) = x⁸ + x⁴ + x³ + x + 1, and x⁶ + x⁴ + x + 1 = {53} in big-endian hexadecimal notation is the element whose inverse is desired, then performing the algorithm results in the following:

i	remainder[i]	quotient[i]	auxiliary[i]
1	x8 + x4 + x3 + x + 1		0
2	x6 + x4 + x + 1		1
3	x2	x2 + 1	x2 + 1
4	x + 1	x4 + x2	x6 + x4 + x4 + x2 + 1
5	1	x + 1	x7 + x6 + x3 + x2 + x2 + x + 1 + 1
6	0	x + 1

Note: Addition in a binary finite field is XOR.

Thus, the inverse is x⁷ + x⁶ + x³ + x = {CA}, as can be confirmed by multiplying the two elements together.

The case of more than two numbers

One can handle the case of more than two numbers iteratively. First we show that $\gcd(a,b,c) = \gcd(\gcd(a,b),c)$. To prove this let $d=\gcd(a,b,c)$. By definition of gcd $d$ is a divisor of $a$ and $b$. Thus $\gcd(a,b)=k d$ for some $k$. Similarly $d$ is a divisor of $c$ so $c=jd$ for some $j$. Let $u=\gcd(k,j)$. By our construction of $u$, $ud | a,b,c$ but since $d$ is the greatest divisor $u$ is a unit. And since $ud=\gcd(gcd(a,b),c)$ the result is proven.

So if $na + mb = \gcd(a,b)$ then there are $x$ and $y$ such that $x\gcd(a,b) + yc = \gcd(a,b,c)$ so the final equation will be

$x(na + mb) + yc = (xn)a + (xm)b + yc = \gcd(a,b,c).\,$

So then to apply to n numbers we use induction

$\gcd(a_1,a_2,\dots,a_n) =\gcd(a_1,\, \gcd(a_2,\, \gcd(a_3,\dots, \gcd(a_{n-1}\,,a_n)))\dots),$

with the equations following directly.

See also

• Euclidean domain
• Linear congruence theorem

References

• Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. Introduction to Algorithms, Second Edition. MIT Press and McGraw-Hill, 2001. ISBN 0-262-03293-7. Pages 859–861 of section 31.2: Greatest common divisor.

External links

• Source for the form of the algorithm used to determine the multiplicative inverse in GF(2^8)