Extended Euclidean Algorithm (a.k.a. the Pulverizer)

Sam

Oct 05 2020 at 15:14 GMT

With Euclid's algorithm, we can find the greatest common divisor (GCD) of two integers $a$ and $b$ . It can be proven that the GCD is the smallest positive integer linear combination of $a$ and $b$ . Thus, we can write the GCD as a linear combination of $a$ and $b$ :

\gcd(a, b) = sa + tb \ \ \ \text { (for some } s,t \text{)}

The extended Euclidean algorithm allows to find not only the GCD, but also the values of the coefficients $s$ and $t$ .

I would like some explanations of how this algorithm works.

number-theory gcd linear-combination pulverizer euclid-algorithm

1 Answer

Carl

Oct 05 2020 at 17:47 GMT

The idea of the extended Euclidean algorithm is to keep track of how each encountered remainder can be written as a linear combination of $a$ and $b$ . This way, once we reach the last non-zero remainder, which corresponds to the GCD, we know how to express it as a linear combination of $a$ and $b$ .

Example: gcd(221, 91)

Let's find the GCD of 221 and 91 using the extended Euclidean algorithm.

\gcd(221, 91) = \gcd(91, \text{rem}(221, 91))

According to the Division theorem, we can write 221 as 91 times the quotient plus the remainder:

221 = 2 \cdot91 + 39

Thus, $\text{rem}(221, 91) = 39$ .

We can express this first remainder (39) as a linear combination of 221 and 91 by rearranging the terms:

39 = 221 + (-2)\cdot91

Let's continue with the algorithm:

\gcd(91, 39) = \gcd(39, \text{rem}(91, 39))

By the Division theorem:

91 = 2 \cdot 39 + 13

So, $\text{rem}(91, 39) = 13$ .

Let's express 13 as a linear combination of 221 and 91:

\begin{aligned} 13 &= 91 - 2\cdot 39 \\ &= 91 - 2 \cdot (221 - 2 \cdot 91) \\ &= 91 - 2 \cdot 221 + 4 \cdot 91 \\ &= -2 \cdot 221 + 5 \cdot 91 \end{aligned}

Let's continue:

\gcd(39, 13) = \gcd(13, \text{rem}(39, 13))

By the Division theorem:

39 = 3 \cdot 13 + 0

Thus, we reached a zero remainder: $\text{rem}(39, 13) = 0$ .

\gcd(13, 0) = 13

Thus, the GCD of 221 and 91 is 13, and we know how to express 13 as a linear combination of 221 and 91:

13 = -2 \cdot 221 + 5 \cdot 91

That is, $s=-2$ and $t=5$ are two possible coefficients in a linear combination of 221 and 91 that makes up the GCD.

C++ Implementation

Let's now see how to implement a function extendedGCD in C++ that computes both the GCD and the coefficients $s$ and $t$ using the extended Euclidean algorithm.

First we need a way to express the return value of such a function. It consists of 3 items:

The GCD
The coefficient $s$ in $sa+tb$
The coefficient $t$ in $sa + tb$

One way to express this is using a struct that has 3 fields. Since this struct represents the value of a linear combination along with the coefficients $s$ and $t$ , let's call it LinearCombination and define it as follows:

struct LinearCombination {
  int value;
  int s;
  int t;
};

So, the GCD of the previous example would correspond to:

LinearCombination gcd = {13, -2, 5};

We now have a way to express the GCD and the coefficients, so the signature of the extendedGCD function will be:

LinearCombination extendedGCD(int a, int b);

Let's think about how the extendedGCD function would work.

We start with $x = a$ and $y = b$ , which trivially are linear combinations of $a$ and $b$ :

\begin{aligned} a &= a + 0 \cdot b \\ b &= 0 \cdot a + b \end{aligned}

We compute the remainder using the Division theorem:

\begin{aligned} x &= qy + r \\ r &= x - qy \end{aligned}

Since $r$ is equal to the difference between linear combinations of $a$ and $b$ , $r$ is also a linear combination of $a$ and $b$ .

Next, we recursively compute $\gcd(y,r)$ . Thus, $y$ becomes the new $x$ and $r$ becomes the new $y$ .

Since $x$ and $y$ are initially linear combinations of $a$ and $b$ , and the remainder $r$ that we compute is also a linear combination of $a$ and $b$ , it follows that $x$ and $y$ are always linear combinations of $a$ and $b$ :

\begin{aligned} x &= s_x a + t_x b \\ y &= s_y a + t_y b \end{aligned}

Knowing this, we can determine how to compute the coefficients $s$ and $t$ for the remainder $r$ :

\begin{aligned} r &= x - qy \\ &= (s_x a + t_x b) - q(s_y a + t_y b) \\ &= s_x a + t_x b - qs_y a - qt_y b \\ &= (s_x - qs_y)a + (t_x - qt_y)b \end{aligned}

Thus, the coefficients for the remainder $r$ are $s = s_x - qs_y$ and $t = t_x - qt_y$ .

We keep doing this until we reach a zero-remainder, that is, until $y$ becomes 0. When that happens, we simply return $x$ , which is the last non-zero remainder, and thus it corresponds to the GCD.

Let's look at the implementation:

struct LinearCombination {
  int value;  // the value of sa + tb
  int s;  // the coefficient s in sa + tb
  int t;  // the coefficient t in sa + tb
};


LinearCombination extendedGCDRec(LinearCombination x, LinearCombination y) {
  if (y.value == 0) {
    return x;
  }

  int q = x.value / y.value;
  int s = x.s - q * y.s;
  int t = x.t - q * y.t;

  LinearCombination r = {x.value - q * y.value, s, t};

  return extendedGCDRec(y, r);
}


LinearCombination extendedGCD(int a, int b) {
  if (a == 0 && b == 0) {
    return LinearCombination({-1, 0, 0});
  }

  if (a == 0) {
    return LinearCombination({b, 0, 1});
  }

  LinearCombination x = {a, 1, 0};
  LinearCombination y = {b, 0, 1};

  return extendedGCDRec(x, y);
}

In extendedGCD, we first check whether both a and b are zero. If that's the case, the GCD is undefined, so we return a linear combination with the sentinel value -1 to indicate that.

Then, we check whether a is zero. If that's the case, the GCD is simply b with s = 0 and t = 1.

Otherwise, we express the initial values of x and y as a linear combinations of a and b, and call extendedGCDRec, which is the recursive function that does all the work.

In extendedGCDRec, we check if y is zero. If that's the case, we found the GCD, which is x.

Otherwise, we compute the remainder and its coefficients s and t.

Then, we call the function recursively with y becoming the new x and r becoming the new y.

We keep doing this until we reach a zero-remainder, in which case the condition y == 0 will be true and we return the last non-zero remainder x.