## Sunday, July 26, 2015

### Extended Euclidean Algorithm

Extended Euclidean Algorithm is an extension of Euclidean Algorithm which finds two things for integer $a$ and $b$:
1. It finds the value of $GCD(a,b)$
2. It finds two integers $x$ and $y$ such that, $ax + by = gcd(a,b)$.
The expression $ax + by = gcd(a,b)$ is known as Bezout's identity and the pair $(x,y)$ that satisfies the identity is called Bezout coefficients. We will look into Bezout's identity at the end of this post. For now just know the name.

## How It Works

In Euclidean Algorithm we worked with remainders only. In Extended Euclidean Algorithm (ext_gcd) we will use the quotient and few other extra variables. Suppose we are finding $GCD(240,46)$. Using Euclidean Algorithm the steps for finding $GCD$ would be:

$GCD(240,46) = GCD(46, 240 \: \% \: 46) = GCD(46,10)$
$GCD(46,10) = GCD(10, 46 \: \% \: 10) = GCD(10,6)$
$GCD(10,6) = GCD(6, 10 \: \% \: 6) = GCD(6,4)$
$GCD(6,4) = GCD(4, 6 \: \% \: 4) = GCD(4,2)$
$GCD(4,2) = GCD(2, 4 \: \% \: 2) = GCD(2,0)$
$GCD(2,0) = 2$

### Introducing $r$ and $q$

We will slowly move towards ext_gcd algorithm. Let us add two new variables. $r$ represents remainder and $q$ means quotient.

Let $r_0$ be $a$ and $r_1$ be $b$. At each step, we will calculate $r_i$. Let $q_i = \lfloor \frac{r_{i-2}}{ r_{i-1}} \rfloor$. Therefore, $r_i = r_{i-2} - q_i \times r_{i-1}$.

Then the above steps will look like the following:
Index i Quotient $q_i$ Remainder $r_i$
$0$ $240$
$1$ $46$
$2$ $240 / 46 = 5$ $240 - 5 \times 46 = 10$
$3$ $46 / 10 = 4$ $46 - 4 \times 10 = 6$
$4$ $10 / 6 = 1$ $10 - 1 \times 6 = 4$
$5$ $6 / 4 = 1$ $6 - 1 \times 4 = 2$
$6$ $4 / 2 = 2$ $4 - 2 \times 2 = 0$
This table is same as the calculation for Euclidean algorithm except for few extra details. Note that the line before last ( index $5$ ) contains the $GCD(a,b)$.

### Introducing $x$ and $y$

We are trying to express $GCD(a,b) = ax + by$. So the variable $x$ and $y$ will hold the coefficients. To be exact we will write each row as terms of $a$ and $b$, i.e, $r_i = ax_i + by_i$.

Initially, $(x_0, y_0) = (1,0)$ and $(x_1, y_1) = (0,1)$. But how do we calculate $(x_i,y_i)$?

We know that $r_i = r_{i-2} - q_i \times r_{i-1}$. We also claimed that  $r_i = ax_i + by_i$. By combining this two we get :

$r_i = ( ax_{i-2} + by_{i-2} ) - q_i \times ( ax_{i-1} + by_{i-1} )$
$r_i = ax_{i-2} - q_i \times ax_{i-1} + by_{i-2} - q_i \times by_{i-1}$
$r_i = a ( x_{i-2} - q_i \times x_{i-1} ) + b ( y_{i-2} - q_i \times y_{i-1})$
$r_i = a x_i + b y_i$
$\therefore x_i = x_{i-2} - q_i \times x_{i-1} \: \text{and} \: y_i = y_{i-2} - q_i \times y_{i-1}$.

Our new table enhanced with $x$ and $y$ will now look like:
Index i Quotient $q_i$ Remainder $r_i$ $x_i$ $y_i$
$0$ $240$ $1$ $0$
$1$ $46$ $0$ $1$
$2$ $240 / 46 = 5$ $240 - 5 \times 46 = 10$ $1 - 5 \times 0 = 1$ $0 - 5 \times 1 = -5$
$3$ $46 / 10 = 4$ $46 - 4 \times 10 = 6$ $0 - 4 \times 1 = -4$ $1 - 4 \times -5 = 21$
$4$ $10 / 6 = 1$ $10 - 1 \times 6 = 4$ $1 − 1 × −4 = 5$ $−5 − 1 × 21 = −26$
$5$ $6 / 4 = 1$ $6 - 1 \times 4 = 2$ $−4 − 1 × 5 = −9$ $21 − 1 × −26 = 47$
$6$ $4 / 2 = 2$ $4 - 2 \times 2 = 0$ $5 − 2 × -9 = 23$ $−26 − 2 × 47 = −120$
Our answer lies on the line before last. $240 \times -9 + 46 \times 47 = 2$.

So all we need to do now is implement these steps in code.

## Code

Even though we will be calculating many rows in ext_gcd algorithm, in order to calculate any row we just need information from previous two rows. So we can save memory by simply storing $r_{i-2}, r_{i-1}, x_{i-2}, x_{i-1}, y_{i-2}, y_{i-1}$.

In our code, $x2$ is $x_{i-2}$, $x1$ is $x_{i-1}$ and $x$ is $x_i$. Same style is followed for other variable.
int ext_gcd ( int A, int B, int *X, int *Y ){
int x2, y2, x1, y1, x, y, r2, r1, q, r;
x2 = 1; y2 = 0;
x1 = 0; y1 = 1;
for (r2 = A, r1 = B; r1 != 0; r2 = r1, r1 = r, x2 = x1, y2 = y1, x1 = x, y1 = y ) {
q = r2 / r1;
r = r2 % r1;
x = x2 - (q * x1);
y = y2 - (q * y1);
}
*X = x2; *Y = y2;
return r2;
}

In line $1$ we define the function. The function $ext\_ gcd()$ takes 4 parameters. The first two parameter $A$ and $B$ represents the two integers whose gcd we wish to find. The next two parameter $*X$ and $*Y$ are two integer pointers. Our function returns us $GCD(A,B)$ but we also want to find the coefficients $x,y$ such that $ax + by = GCD(A,B)$. So we extract those values using pointers.

In line $2$, we declare all necessary variables. We initiate some variables in line $3$ and $4$. Then a for loop starts in line $5$. We initiate few more variables in first section, $r2$ and $r1$. The loop runs till $r1$ becomes 0. In the last section of for loop, we make transitions of variables, such as $r2$ becomes $r1$ and $r1$ gets the new value $r$.

In side the for loop we calculate values needed for current row, $q, r, x, y$.

When the loop finishes, $r1$ contains $0$ and $r2$ is the row before it containing $GCD(A,B)$. So we send $x2,y2$ as coefficients.

## Complexity

Same as Euclidean Algorithm.

## Uniqueness of Solution

Using ext_gcd we found $(x,y)$ pair which satisfies $ax + by = gcd(a,b)$. But is it unique?

No. Once we find a pair $(x,y)$ using ext_gcd, we can generate infinite pairs of Bezout coefficients using the formula: $$( x + k \frac { b } { \text{gcd}(a,b)}, y - k \frac { a } { \text{gcd}(a,b)} )$$Using any integer value for $k$ we can generate a different pair of Bezout coefficients $(x,y)$ which will satisfy the Bezout's identity. Here is why it works:

$a ( x + \frac { kb } { \text{gcd}(a,b)} ) + b ( y - \frac { ka } { \text{gcd}(a,b)} )$
$ax + \frac { kab } { \text{gcd}(a,b)} + by - \frac { kab } { \text{gcd}(a,b)}$
$ax + by$

As you can see above, the terms with $k$ in them cancel each other out. They don't change the expression $ax + by$ in any way, therefore, there are infinite pairs of Bezout coefficients.

## Smallest Positive Integer of form $ax + by$

We showed that it is possible to write $ax + by = gcd(a,b)$, now we need to find a positive number smaller than $gcd(a,b)$ of form $ax + by$. Is that even possible?

No. $gcd(a,b)$ divides both $a$ and $b$. Therefore, if a number is of form $ax + by$ it will be divisible by $gcd(a,b)$ since $ax + by$ is divisible by $gcd(a,b)$. And the smallest positive number which is divisible by $gcd(a,b)$ is $gcd(a,b)$ itself.

So $gcd(a,b)$ is the smallest positive number of the form $ax + by$.

## Bezout's Identity

This was mentioned at the beginning of the post. Almost everything related to Bezout's Identity has been explained. I will still mention them once more for the sake of completeness.

Bézout's identity (also called Bézout's lemma) is a theorem in the elementary theory of numbers: let a and b be nonzero integers and let d be their greatest common divisor. Then there exist integers x and y such that $ax + by = d$

• the greatest common divisor d is the smallest positive integer that can be written as $ax + by$
• every integer of the form ax + by is a multiple of the greatest common divisor d.
-Wiki
Bezout's Lemma simply states that $ax + by = gcd(a,b)$ exists. We need to use Extended Euclidean Algorithm to find Bezout's Coefficients.

## Coding Pitfalls

Careful about the equation $ax + by = \text{gcd} (a,b)$. Here $\text{gcd}()$ means the result from Euclidean Algorithm, not what we mean mathematically. For example, when $a = 4$ and $b = -2$, Extended Euclidean Algorithm finds solution for $4x - 2y = -2$. According to Euclidean algorithm $gcd(4,-2)$ returns $-2$, though in common sense it should be $2$.

## Reference

#### 1 comment:

1. Good explanatory resource for extended euclid is hard to come by. Thanks :)

Leave comments for Queries, Bugs and Hugs.