The extended Euclidian algorithm

The extended Euclidian algorithm is based on the Euclidian algorithm, which finds the greatest common divisor of two numbers $a, b$ . The extended Euclidian algorithm also does some bookkeeping, which allows it to express the greatest common divisor as a linear combination of $a$ and $b$ . For example, if we run the Euclidian algorithm for $a = 112, b = 42$ , we find that $g cd (112, 42) = 14$ . Now, running the Euclidian algorithm tells us that $g cd (112, 42) = 14$ and additionally tells us that $14 = - 1 \cdot 112 + 3 \cdot 42$ .

Theory

If we suppose that the Euclidian algorithm is used to find the greatest common divisor of $a_{0}, b_{0}$ where $a_{0} \geq b_{0}$ , we can formulate the algorithm as a recurrent relation

a_{n + 1} = b_{n}

b_{n + 1} = mod_{b_{n}} (a_{n})

which is then applied until $b_{n} = 0$ , in which case $a_{n}$ is the greatest common divisor $g cd (a_{0}, b_{0})$ of $a_{0}$ and $b_{0}$ .

The extended Euclidian algorithm now keeps numbers $p_{n}, q_{n}, r_{n}, s_{n}$ such that

a_{n} = p_{n} a_{0} + q_{n} b_{0}

b_{n} = r_{n} a_{0} + s_{n} b_{0}

Obviously, we have $p_{0} = 1, q_{0} = 0$ and $r_{0} = 0, s_{0} = 1$ . Using that $a_{n + 1} = b_{n}$ , we see that

p_{n + 1} = r_{n}

q_{n + 1} = s_{n}

To find $r_{n + 1}$ and $s_{n + 1}$ , use that $b_{n + 1} = mod_{b_{n}} (a_{n}) = a_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ b_{n}$ . Using $a_{n} = p_{n} a_{0} + q_{n} b_{0}$ and $b_{n} = r_{n} a_{0} + s_{n} b_{0}$ we see that

b_{n + 1} = a_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ b_{n} = (p_{n} a_{0} + q_{n} b_{0}) - ⌊ \frac{a _{n}}{b _{n}} ⌋ (r_{n} a_{0} + s_{n} b_{0})

factoring out $a_{0}$ and $b_{0}$ , we see that

b_{n} = (p_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ r_{n}) a_{0} + (q_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ s_{n}) b_{0}

so we see that

r_{n + 1} = p_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ r_{n}

s_{n + 1} = q_{n} - ⌊ \frac{a _{n}}{b _{n}} ⌋ s_{n}

Implementation

Normally I try to give an implementation that works in C-style languages. However, the C-style languages are rather primitive in the sense that they provide no uniform and comfortable way to use tuples, so I will use Python for now.

Like in the implementation for Euclids algorithm, we have a wrapper that ensures that the first argument is always the biggest. The wrapper returns a triple (gcd, p, q), such that gcd is the greatest common divisor, and we have $g cd (a, b) = p a + q b$ .

def gcd_extended(a, b):
	if a < b:
		(gcd, p, q) = gcd_extended_helper(b, 1, 0, a, 0, 1)
		return (gcd, q, p)
	else: return gcd_extended_helper(a, 1, 0, b, 0, 1)

Like for Euclids algorithm, we can use a recursive implementation

def gcd_extended_helper(a, p, q, b, r, s):
	if b == 0: return (a, p, q)
	k = a // b;
	return gcd_extended_helper(b, r, s, a - k * b, p - k * r, q - k * s)

or an iterative one

def gcd_extended_helper(a, p, q, b, r, s):
	while b != 0:
		(a, p, q, b, r, s) = (b, r, s, a - k * b, p - k * r, q - k * s)
	return (a, p, q)

Applications

When modular arithmetic modulo some number $n$ , we may want to find a multiplicative inverse $a^{- 1}$ of the element $a$ . That is, we want to find $a^{- 1}$ such that

a a^{- 1} \equiv 1 (mod n)

Now, such an element only exists if $g cd (a, n) = 1$ , which is why it is often convenient to take $n$ prime (we can take any $a = 1, 2, ..., n - 1$ ). Now, note that if $p a + r n = 1$ , then we have $a p \equiv 1 (mod n)$ . So $a^{- 1} \equiv p (mod n)$ is the element we’re looking for.

Now we can define a function that computes the multiplicative inverse modulo $n$ as follows:

def inverse_modulo(a, n):
	(gcd, p, q) = gcd_extended(a, n)
	return p