Extended Euclidean algorithm without stack or recursion

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 not only to calculate the gcd (greatest common divisor) of 2 numbers, but gives us also a representation of the result in a form of a linear combination:

\gcd(a, b) = u \cdot a + v \cdot b \quad u,v \in \mathbb{Z}

gcd of more than 2 numbers can always be done by iteratively calculating the gcd of 2 numbers.

For example, let us calculate gcd(14, 5):

\begin{aligned}
14 &= 5 \cdot 2 + 4 \\
5 &= 4 \cdot 1 + 1 \\
4 &= 1 \cdot 4 + 0
\end{aligned}

So the greatest common divisor of 14 and 5 is 1.

We can find the linear combination coefficients by writing 1 in terms of 14 and 5:

\begin{aligned}
1 &= 5 - 4 \cdot 1 \\
&= 5 - (14 - 5 \cdot 2) \cdot 1 \\
&= 5 - 14 + 5 \cdot 2 \\
&= 3 \cdot 5 + (-1) \cdot 14
\end{aligned}

So in this case u = 3 and v = -1:

\gcd(14, 5) = (-1) \cdot 14 + 3 \cdot 5 = 1

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.

Recursive implementation

The formula

\gcd(a, b) =
	\begin{cases}
	b, & \text{if}\ a = 0 \\
	\gcd(b \bmod a, a), & \text{otherwise}
	\end{cases}

allows us to describe the algorithm in a functional way:

  • If a = 0, then the greatest common divisor is b. Coefficients u = 0 and v = 0.
  • Else, we make the problem simpler by calculating \gcd(b \bmod a, a). We can calculate the new coefficients based on the coefficients of the simpler problem.

So, how can we calculate u and v so that

\gcd(a, b) = u \cdot a + v \cdot b

by knowing u' and v' with:

\gcd(b \bmod a, a) = u' \cdot (b \bmod a) + v' \cdot a

In order to do that we can write b \bmod a in terms of initial a and b:

\begin{aligned}
	\gcd(b \bmod a, a)
    	&= u' \cdot (b \bmod a) + v' \cdot a \\
    	&= u' \cdot (b - \left\lfloor \frac{b}{a} \right\rfloor \cdot a) + v' \cdot a \\
    	&= u' \cdot b - u' \cdot \left\lfloor \frac{b}{a} \right\rfloor \cdot a + v' \cdot a \\
    	&= (v' - u' \cdot \left\lfloor \frac{b}{a} \right\rfloor) \cdot a + u' \cdot b
\end{aligned}

So the new linear combination coefficients are:

\begin{aligned}
    u &= v' - u' \cdot \left\lfloor \frac{b}{a} \right\rfloor \\
    v &= u'
\end{aligned}

With this formula we are now ready to implement the algorithm:

class GCD_Result: # Representation of the result
    def __init__(self, gcd, u, v):
        self.gcd = gcd
        self.u = u
        self.v = v

def extended_gcd(a, b):
    if a == 0:
        return GCD_Result(b, 0, 1)
    result = extended_gcd(b % a, a)
    u = result.u # save u'
    result.u = result.v - (b // a) * result.u # u = v' - u' * (b // a)
    result.v = u # v = u'
    return result

Non-recursive implementation

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 corresponding remainers as we calculate the remainders. Formally, we can define f finite sequence r_n:

\begin{aligned}
r_1 &= a \\
r_2 &= b \\
r_{n+2} &= r_n \bmod r_{n+1}
\end{aligned}

If r_{n+1} = 0, r_{n+2} is not defined. We can write each r_n as a linear combination of u and v. Now we are interested in how u and v change as we calculate remainders. To do this formally, we will need to define two new finite sequences u_n and v_n which will represent the linear combination coefficients:

r_n = u_n \cdot a + v_n \cdot b

By definition, r_1 = a and r_2 = b, so we can directly write the linear combination coefficients for r_1 and r_2:

\begin{aligned}
    u_1 &= 1 \\
    v_1 &= 0 \\
    u_2 &= 0 \\
    v_2 &= 1
\end{aligned}

Let q_n be the finite sequence of integer divisions in r_n:

r_n = r_{n+1} \cdot q_{n+2} + r_{n+2}

Now we can write u_n and v_n in terms of q_n:

\begin{aligned} 
    r_{n+2} &= r_n - r_{n+1} \cdot q_{n+2} \\
    &= u_n \cdot a + v_n \cdot b - r_{n+1} \cdot q_{n+2} \\
    &= u_n \cdot a + v_n \cdot b - (u_{n+1} \cdot a + v_{n+1} \cdot b) \cdot q_{n+2} \\
    &= u_n \cdot a + v_n \cdot b - u_{n+1} \cdot a \cdot q_{n+2} - v_{n+1} \cdot b \cdot q_{n+2} \\
    &= (u_n - u_{n+1} \cdot q_{n+2}) \cdot a + (v_n - v_{n+1} \cdot q_{n+2}) \cdot b
\end{aligned}

To get the formula for u_n and v_n we can just substitute n instead of n + 2:

\begin{aligned}
    u_n &= u_{n-2} - q_n \cdot u_{n-1} \\
    v_n &= v_{n-2} - q_n \cdot v_{n-1}
\end{aligned}

With this formula and the initial values of the u_n and v_n sequences we can now implement the extended Euclidean algorithm without recursion:

def extended_gcd(a, b):
    if a == 0:
        # The algorithm will work correctly without this check
        # But it will take one iteration of the inner loop
        return GCD_Result(b, 0, 1)

    unPrev = 1
    vnPrev = 0
    unCur = 0
    vnCur = 1

    while b != 0:
        # Calculate new element of the qn sequence
        qn = a // b
        
        # Calculate new element of the rn sequence
        newRemainder = a % b
        a = b
        b = newRemainder

        # Calculate new coefficients with the formula above
        unNew = unPrev - qn * unCur
        vnNew = vnPrev - qn * vnCur

        # Shift coefficients
        unPrev = unCur
        vnPrev = vnCur
        unCur = unNew
        vnCur = vnNew

    return GCD_Result(a, unPrev, vnPrev)

Example

We can visualize the finite sequences we defined and see how the algorithm works with a table. We will calculate \gcd(104, 47) and it’s linear combination coefficients:

\gcd(104, 47) = u \cdot 104 + v \cdot 47
r_nq_nu_nv_n
104-10
47-01
1021-2
74-49
315-11
12-1431
033320

At each step we first calculate the next element from the q_n sequence and then use it to calculate new linear combination coefficients u_n and v_n.

The result of the algorithm:

\gcd(104, 47) = -14 \cdot 104 + 31 \cdot 47 = 1

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 r_n sequence:

def extended_gcd(a, b):
    if a == 0: # Optional check
        return GCD_Result(b, 0, 1)

    if b == 0: # Without this check the first iteration will divide by zero
        return GCD_Result(a, 1, 0)

    unPrev = 1
    vnPrev = 0
    unCur = 0
    vnCur = 1

    while True:
        qn = a // b
        newR = a % b
        a = b
        b = newR

        if b == 0:
            return GCD_Result(a, unCur, vnCur)

        # Update coefficients
        unNew = unPrev - qn * unCur
        vnNew = vnPrev - qn * vnCur

        # Shift coefficients
        unPrev = unCur
        vnPrev = vnCur
        unCur = unNew
        vnCur = vnNew

Copyright © 2020 Alexander Mayorov. All rights reserved.