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
sandtsuch thatgcd(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) = aand(b, 0) = b.(2·a, 2·b) = 2·gcd(a, b)(2·a, b) = (a, b), ifbis odd (2 is not a common divisor).(a, b) = (|a − b|, min(a, b)), ifaandbare 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