Extended Euclidean algorithm without stack or recursion

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)=ua+vbu,vZ\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’s calculate gcd(14,5)\gcd(14, 5): 14=52+45=41+14=14+0\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 1414 and 55 is 11.

We can find the linear combination coefficients by writing 11 in terms of 1414 and 55: 1=541=5(1452)1=514+52=35+(1)14\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=3u = 3 and v=1v = -1: gcd(14,5)=(1)14+35=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)={b,if a=0gcd(bmoda,a),otherwise\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=0a = 0, then the greatest common divisor is bb. Coefficients u=0u = 0 and v=0v = 0.
  • Else, we make the problem simpler by calculating gcd(bmoda,a)\gcd(b \bmod a, a). We can calculate the new coefficients based on the coefficients of the simpler problem.

So, how can we calculate uu and vv so that gcd(a,b)=ua+vb\gcd(a, b) = u \cdot a + v \cdot b

by knowing uu' and vv' with: gcd(bmoda,a)=u(bmoda)+va\gcd(b \bmod a, a) = u' \cdot (b \bmod a) + v' \cdot a

In order to do that we can write bmodab \bmod a in terms of initial aa and bb: gcd(bmoda,a)=u(bmoda)+va=u(bbaa)+va=ububaa+va=(vuba)a+ub\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: u=vubav=u\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:

def extended_gcd(a, b):
    if a == 0:
        return (b, 0, 1)
    (g, u, v) = extended_gcd(b % a, a)
    newU = v - (b // a) * u
    newV = u
    return (g, newU, newV)

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 coefficients as we calculate the remainders. Formally, we can define a finite sequence rnr_n: r1=ar2=brn+2=rnmodrn+1\begin{aligned} r_1 &= a \\ r_2 &= b \\ r_{n+2} &= r_n \bmod r_{n+1} \end{aligned}

If rn+1=0r_{n+1} = 0, rn+2r_{n+2} is not defined. We can write each rnr_n as a linear combination of uu and vv. Now we are interested in how uu and vv change as we calculate remainders. To do this formally, we will need to define two new finite sequences unu_n and vnv_n which will represent the linear combination coefficients: rn=una+vnbr_n = u_n \cdot a + v_n \cdot b

By definition, r1=ar_1 = a and r2=br_2 = b, so we can directly write the linear combination coefficients for r1r_1 and r2r_2: u1=1v1=0u2=0v2=1\begin{aligned} u_1 &= 1 \\ v_1 &= 0 \\ u_2 &= 0 \\ v_2 &= 1 \end{aligned}

Let qnq_n be the finite sequence of integer divisions in rnr_n: rn=rn+1qn+2+rn+2r_n = r_{n+1} \cdot q_{n+2} + r_{n+2}

Now we can write unu_n and vnv_n in terms of qnq_n: rn+2=rnrn+1qn+2=una+vnbrn+1qn+2=una+vnb(un+1a+vn+1b)qn+2=una+vnbun+1aqn+2vn+1bqn+2=(unun+1qn+2)a+(vnvn+1qn+2)b\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 unu_n and vnv_n we can just substitute nn instead of n+2n + 2: un=un2qnun1vn=vn2qnvn1\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 unu_n and vnv_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 (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 (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)\gcd(104, 47) and it’s linear combination coefficients: gcd(104,47)=u104+v47\gcd(104, 47) = u \cdot 104 + v \cdot 47

rnr_nqnq_nunu_nvnv_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 qnq_n sequence and then use it to calculate new linear combination coefficients unu_n and vnv_n.

The result of the algorithm: gcd(104,47)=14104+3147=1\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 rnr_n sequence:

def extended_gcd(a, b):
    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)

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

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

        if b == 0:
            return (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.