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

References

[^1] Handbook of Applied Cryptography