# Euclidean Algorithm

Created:
2017-04-16

Updated:
2017-04-16

A fundamental and efficient algorithm to compute the greatest common divisor
(`gcd`

from now on) of two integer numbers which forms the basis of several
cryptographic algorithms.

Given two integers `a`

and `b`

, the algorithm comes in two variants:

**basic**: yields the`gcd(a, b)`

;**extended**: yields two integers`s`

and`t`

such that`gcd(a,b) = a·s + b·t`

As a notational shortcut we’ll write `(a, b)`

instead of `gcd(a, b)`

.

## Basic Euclidean Algorithm

Assuming `a > b`

:

```
d | a → a = d·x
d | b → b = d·y
→ d | (a - b) = d·(x - y)
```

The procedure is iterated while `a ≥ b·q`

, with `q > 0`

an integer such that
`a < b·(q + 1)`

```
b = d·y → b·q = d·y·q = d·k → d | b·q
→ d | (a - b·q) = d·(x - k)
```

If `q`

is the bigger integer such that `a ≥ b·q`

then the above is saying that
`d`

divides the remainder of `a`

divided by `b`

. In other words `d | (a mod b)`

The described procedure is equivalent to the following lemma:

**Lemma (Euclid’s algorithm)**. Given `b ≥ 1`

, `(a,b) = (b, a mod b)`

```
d|a ∧ d|b ↔ d|b ∧ d|(a mod b)
```

*Proof*

(⇒) `d|a ∧ d|b`

```
a = d·x and b = d·y
a mod b = a - b·q , for max{q: a - b·q ≥ 0}
a - b·q = d·(x - y·q) → d | (a mod b)
```

(⇐) `d|b ∧ d|(a mod b)`

```
b = d·y
d|(a mod b) = a - b·q → a - b·q = d·x
→ a = d·x + q·b = d·x + q·d·y = d·(x + q·y) → d|a
```

∎

Note that if `b > a`

then `(a, b) = (b, a mod b) = (b, a)`

.
That is, in this case the Euclidean algorithm step has the effect of a swap.

By repeating the procedure, since the second number steadily decreases, we’ll
finally reach the point where it is zero. And `(k, 0) = k`

.

Example:

```
a = 3^3 · 5 · 11^2 = 16335
b = 2 · 3 · 5^2 · 7 = 1050
(16335, 1050) = (1050, 16335 mod 1050) =
(1050, 585) = (585, 1050 mod 585) =
(585, 465) = (465, 585 mod 465) =
(465, 120) = (120, 465 mod 120) =
(120, 105) = (105, 120 mod 105) =
(105, 15) = (15, 105 mod 15) =
(15, 0) = 15
```

## Extended Euclidean Algorithm (EEA)

**Theorem. Bezout’s Identity**

Given two integers `a`

and `b`

with `b ≥ 1`

and `d = (a, b)`

there exists two
integers `x`

and `y`

such that `d = a·x + b·y`

.

*Proof*.

By construction using the basic Euclid’s algorithm.

We initially set `r0 = a`

and `r1 = b`

```
r0 = q1·r1 + r2 (r0, r1) = (r1, r0 mod r1) = (r1, r2)
r1 = q2·r2 + r3 (r1, r2) = (r2, r1 mod r2) = (r2, r3)
...
r(n-2) = q(n-1)·r(n-1) + rn
r(n-1) = qn·rn + r(n+1) (r(n-1), rn) = (rn, r(n+1))
With r(n+1) = 0 → (rn, 0) = rn = d
```

Extended:

```
r0 = a
r1 = b
r2 = r0 - q1·r1
r3 = r1 - q2·r2
...
rn = r(n-2) - q(n-1)·r(n-1)
```

We can thus write `rn`

as a linear combination of `a`

and `b`

.

```
r0 = x0·a + y0·b , with x0=1 and y0=0
r1 = x1·a + y1·b, with x1=0 and y1=1
r2 = r0 - q1·r1 = x2·a + y2·b
r3 = r1 - q2·r2 = x3·a + x3·b
...
ri = r(i-2) - q(i-1)·r(i-1) =
= [x(i-2)·a + y(i-2)·b] - q(i-1)·[x(i-1)·a + y(i-1)·b] =
= [x(i-2) - q(i-1)·x(i-1)]·a + [y(i-2) - q(i-1)·y(i-1)]·b =
= xi·a + yi·b
...
rn = ... = a·x + b·y = d
```

∎

Note that:

```
xi = x(i-2) - q(i-1)·x(i-1)
yi = y(i-2) - q(i-1)·y(i-1)
```

**Corollary**. EEA to compute multiplicative inverse.

*Proof*

If `1 = a·x + b·y`

then `1 ≡ a·x (mod b)`

, thus `x`

is the multiplicative
inverse of `a`

modulo `b`

.

∎

Note that the quotients of the Euclidean algorithm should be computed anyway.

Example. Compute the inverse of 60 modulo 17:

```
r0 = 60, r1 = 17
r2 = 60 - 3·17 = 9 (q1 = 3)
r3 = 17 - 1·9 = 8 (q2 = 1)
r4 = 9 - 1·8 = 1 (q3 = 1)
r5 = 8 - 8·1 = 0 (q4 = 8)
→ r4 = (60, 17) = 1
x0 = 1
x1 = 0
x2 = x0 - q1·x1 = 1 - 3·0 = 1
x3 = x1 - q2·x2 = 0 - 1·1 = -1
x4 = x2 - q3·x3 = 1 - 1·(-1) = 2
The inverse is thus 2: 60·2 mod 17 = 1
```

## Efficient Alternative GCD Implementations (HAC[^1] 14.4)

### Binary GCD Algorithm (Stein’s Algorithm)

Stein’s algorithm uses simpler arithmetic operations than the conventional Euclidean algorithm; it replaces division with arithmetic shifts, comparisons, and subtraction.

The algorithm reduces the problem of finding the `gcd`

of two **non-negative**
numbers `a`

and `b`

by repeatedly applying these identities:

`(0, a) = a`

and`(b, 0) = b`

.`(2·a, 2·b) = 2·gcd(a, b)`

`(2·a, b) = (a, b)`

, if`b`

is odd (2 is not a common divisor).`(a, b) = (|a − b|, min(a, b))`

, if`a`

and`b`

are both odd (this last one is a single step of the classic Euclidean algorithm).

Pseudo Code

```
g = 1
while a and b are both even do: a = a/2, b = b/2, g = 2·g
while a ≠ 0 do:
while a is even do: a = a/2
while b is even do: b = b/2
t = |a - b|/2
if a ≥ b then then a = t else b = t
return g·b
```

### Binary Extended Euclid’s Algorithm

The classic EEA algorithm has the drawback of requiring relatively costly multiple-precision divisions. The following algorithm eliminates the divisions at the expense of more iterations.

Pseudo Code

```
1. g = 1
2. while a and b are bot even do:
a = a/2, b = b/2, g = 2·g
3. u = a, v = b
A = 1, B = 0, C = 0, D = 1
4. while u is even do:
u = u/2
if A ≡ B ≡ 0 (mod 2) then:
A = A/2, B = B/2
else:
A = (A + b)/2, B = (B - a)/2
5. while v is even do:
v = v/2
if C ≡ D ≡ 0 (mod 2) then:
C = C/2, D = D/2
else:
C = (C + b)/2, D = (D - a)/2
6. if u ≥ v then:
u = u - v, A = A - C, B = B - D
else
v = v - u, C = C - A, D = D - B
7. if u ≠ 0:
goto step 4
a = C, b = D
return(a, b, g·v)
```

## Python Code

### Basic Euclid’s Algorithm

```
def euclid_gcd(a, b):
while b != 0:
a, b = b, a % b
return a
```

### Recursive EEA

Given:

```
gcd, x1, y1 = EEA(b, a mod b)
```

Such that `gcd = b·x1 + (a mod b)·y1`

then:

```
gcd = b·x1 + (a - b·q)·y1 = b·x1 + a·y1 - b·q·y1 = a·y1 + b·(x1 - q·y1) = a·x2 + b·y2
```

```
def extended_euclidean(a, b):
if b == 0:
return a, 1, 0
gcd, x1, y1 = extended_euclidean(b, a % b)
x2 = y1
y2 = x1 - (a // b) * y1
return gcd, x2, y2
```

### Non-recursive EEA

```
def extended_euclidean(a, b):
x0, y0, x1, y1 = 1, 0, 0, 1
while b != 0:
q = a // b
a, b = b, a - q * b
x0, x1 = x1, x0 - q * x1
y0, y1 = y1, y0 - q * y1
return a, x0, y0
```

### Modular Inverse

Computes the inverse of `a`

modulo `b`

assuming that `(a, b) = 1`

.

```
def inverse(a, b):
x0, x1 = 1, 0
while b != 0:
q = a // b
a, b = b, a - q * b
x0, x1 = x1, x0 - q * x1
return x0
```

## References

- Simple C implementation of modular inverse