Typical implementation of the extended Euclidean algorithm on the internet will just iteratively calculate modulo until 0 is reached. However, sometimes you also need to calculate the linear combination coefficients for the greatest common divisor.
Extended Euclidean algorithm
The extended Euclidean algorithm allows us to not only calculate the greatest common divisor of two numbers, but gives us also a representation of the result in a form of a linear combination:
gcd of more than 2 numbers can always be done by iteratively calculating the gcd of 2 numbers.
For example, let’s calculate :
So the greatest common divisor of and is .
We can find the linear combination coefficients by writing in terms of and :
So in this case and :
We can calculate the linear combination coefficients by doing back substitution. But it is not so easy to implement this without recursion, because the back substitution is done when we are climbing out of the recursive calls. We will implement the algorithm recursively first.
allows us to describe the algorithm in a functional way:
- If , then the greatest common divisor is . Coefficients and .
- Else, we make the problem simpler by calculating . We can calculate the new coefficients based on the coefficients of the simpler problem.
So, how can we calculate and so that
by knowing and with:
In order to do that we can write in terms of initial and :
So the new linear combination coefficients are:
With this formula we are now ready to implement the algorithm:
def extended_gcd(a: int, b: int) -> (int, int, int): if a == 0: return b, 0, 1 (g, u, v) = extended_gcd(b % a, a) new_u = v - (b // a) * u new_v = u return g, new_u, new_v
The recursion in the algorithm above cannot be easily eliminated because the function is not tail-recursive.
In order to implement the algorithm with a loop we need to define a sequence of division remainders and then update the coefficients as we calculate the remainders. Formally, we can define a finite sequence :
If , is not defined. We can write each as a linear combination of and . Now we are interested in how and change as we calculate remainders. To do this formally, we will need to define two new finite sequences and which will represent the linear combination coefficients:
By definition, and , so we can directly write the linear combination coefficients for and :
Let be the finite sequence of integer divisions in :
Now we can write and in terms of :
To get the formula for and we can just substitute instead of :
With this formula and the initial values of the and sequences we can now implement the extended Euclidean algorithm without recursion:
def extended_gcd(a: int, b: int) -> (int, int, int): if a == 0: # The algorithm will work correctly without this check # But it will take one iteration of the inner loop return b, 0, 1 un_prev = 1 vn_prev = 0 un_cur = 0 vn_cur = 1 while b != 0: # Calculate new element of the qn sequence qn = a // b # Calculate new element of the rn sequence new_r = a % b a = b b = new_r # Calculate new coefficients with the formula above un_new = un_prev - qn * un_cur vn_new = vn_prev - qn * vn_cur # Shift coefficients un_prev = un_cur vn_prev = vn_cur un_cur = un_new vn_cur = vn_new return a, un_prev, vn_prev
We can visualize the finite sequences we defined and see how the algorithm works with a table. We will calculate and it’s linear combination coefficients:
At each step we first calculate the next element from the sequence and then use it to calculate new linear combination coefficients and .
The result of the algorithm:
Improvement of the non-recusive solution
As we see in the example above, we don’t need to calculate the last row of the table because we aren’t interested in the linear combination that forms zero. We can terminate the algorithm directly after calculating the new element of the sequence:
def extended_gcd(a: int, b: int) -> (int, int, int): if a == 0: # Optional check return b, 0, 1 if b == 0: # Without this check the first iteration will divide by zero return a, 1, 0 un_prev = 1 vn_prev = 0 un_cur = 0 vn_cur = 1 while True: qn = a // b new_r = a % b a = b b = new_r if b == 0: return a, un_cur, vn_cur # Update coefficients un_new = un_prev - qn * un_cur vn_new = vn_prev - qn * vn_cur # Shift coefficients un_prev = un_cur vn_prev = vn_cur un_cur = un_new vn_cur = vn_new