 # 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:

1. `(0, a) = a` and `(b, 0) = b`.
2. `(2·a, 2·b) = 2·gcd(a, b)`
3. `(2·a, b) = (a, b)`, if `b` is odd (2 is not a common divisor).
4. `(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
``````