7.3 Proofs and Algorithms III: Computing the Greatest Common Divisor

In the previous section, we studied some mathematical properties of the greatest common divisor. Now in this section, we’ll look at how to implement algorithms for calculating the greatest common divisor, and introduce a new form of Python loops along the way.

Naively searching for the gcd

In 4.6 Proofs and Programming I: Divisibility, we developed multiple Python implementations of the divisibility predicate. For example, here is the last implementation, which uses a remainder calculation:

def divides(d: int, n: int) -> bool:
    """Return whether d divides n."""
    if d == 0:
        return n == 0
    else:
        return n % d == 0

With this function in hand, we can implement a gcd function as follows: In this implementation, we use abs because m and n might be negative.

def naive_gcd(m: int, n: int) -> int:
    """Return the gcd of m and n."""
    if m == 0:
        return abs(n)
    elif n == 0:
        return abs(m)
    else:
        possible_divisors = range(1, min(abs(m), abs(n)) + 1)
        return max({d for d in possible_divisors if divides(d, m) and divides(d, n)})

GCD and remainders

The above implementation of naive_gcd is correct, and follows the mathematical definition of gcd quite closely. However, just as we saw when implementing divides, we can develop a better (i.e., faster) implementation of the function by using more advanced mathematical properties. First, recall the Quotient-Remainder Theorem:

Here is the key theorem we’ll use to improve our gcd implementation; this theorem allows us to break down the problem of finding the gcd of two numbers into a smaller one.

GCD, remainders, and the Euclidean Algorithm

The theorem we just stated suggests a possible way of computing the gcd of two numbers in an iterative (repeated) fashion. Let’s again use 24 and 16 as our example.

• Since the remainder \(24~\%~16\) is 8, we know that \(\gcd(24, 16) = \gcd(16, 8)\).
• Also, the remainder \(16~\%~8\) is 0, and so we know that \(\gcd(16, 8) = \gcd(8, 0)\).
• But the gcd of any positive integer \(n\) and 0 is simply \(n\) itself, Exercise: prove this using the definition of gcd! and so we know \(\gcd(8, 0) = 8\).
• This tells use that \(\gcd(24, 16) = 8\) as well!

Let’s formalize this in a high-level description of an algorithm before we write the code. This algorithm for computing the gcd of two numbers is known as the Euclidean Algorithm. This is named after the Greek mathematician Euclid, although he originally developed the algorithm using subtraction (\(a - b\)) rather than remainders (\(a~\%~b\)). To simplify the algorithm, we’ll only cover the version that works for non-negative integers in this section, but this algorithm can be extended to work for all integers.

Here is how we can visualize the changing values of x and y for the given 24 and 6 in our previous example: Note the similarity between this and the loop accumulation tables of Chapter 5.

The main question we must resolve to implement this algorithm in Python is how we achieve step 4: repeating the two previous steps until some condition is satisfied. We know how to use for loops to iterate over a collection of values. This allowed us to repeat a sequence of Python statements (i.e., the body of the for loop) on every iteration. Naturally, the for loop ends when the statements have been repeated for all elements in a collection or range.

But in the case of step 4, we would like to repeat code based on some condition: “Repeat steps 2 and 3 until the remainder is 0”. In these scenarios, we must use a different kind of loop in Python: the while loop.

The while loop

A while loop is a Python compound statement that looks very similar to an if statement:

while <condition>:
    <statement>
    ...

When the Python interpreter executes a while loop, it first executes the while loop condition. If the condition is True, the interpreter executes the statements in the loop body, and if the condition is False, the interpreter skips over the loop body. This is similar to how the Python interpreter executes if statements, but if there is one key difference: after executing the loop body, the Python interpreter checks the loop condition again. If the condition still evaluates to True, then the body is repeated. Let’s try an example:

>>> numbers = []
>>> number = 1
>>> while number < 100:
...     numbers.append(number)
...     number = number * 2
...
>>> numbers
[1, 2, 4, 8, 16, 32, 64]

Notice how number appears in both the while loop’s body and its condition. In the loop body, number is increasing at each iteration (we accumulated the values in the list numbers). Eventually, number refers to the value 128 and the while loop is done because 128 < 100 evaluates to False. Note that the number of iterations of our while loop is dependent on the initial value of number. Had we started with a value of, for example, 10, the loop would have only 4 iterations (not 6, as when number started with 2). Similarly, if number was initially some value greater than or equal to 100, then the while loop would never have executed its body (just as a for loop does not execute its body if given an empty collection).

Implementing the Euclidean Algorithm

Here is our (first) implementation of the Euclidean algorithm for computing the gcd of two numbers. We have labelled the code corresponding to the various steps in the Euclidean Algorithm described above.

def euclidean_gcd(a: int, b: int) -> int:
    """Return the gcd of a and b.

    Preconditions:
    - a >= 0
    - b >= 0
    """
    # Step 1: initialize x and y
    x = a
    y = b

    while y != 0:  # Step 4: repeat Steps 2 and 3 until y is 0
        # Step 2: calculate the remainder of x divided by y
        r = x % y

        # Step 3: reassign x and y
        x = y
        y = r

    # Step 5: x now refers to the gcd of a and b
    return x

How does this loop work? To understand it better, let’s see how this maps onto our original algorithm.

• Step 1, initializing x and y, occurs in the code before the while loop begins.
• Steps 2 and 3 are performed inside the loop body.
• Step 4, the repetition, is achieved by the while loop. One subtlety is that our original algorithm specified a stopping condition, “repeat until X”. When writing Python while loops, however, we must write a continuing condition, which is the negation of the stopping condition. So “until \(y = 0\)” becomes “while y != 0”.
• Step 5, the return value, is exactly what is specified by the algorithm.

Let’s see an example trace of the euclidean_gcd loop for the sample call euclidean_gcd(24, 16):

In this loop, we don’t have a typical accumulator pattern. Instead, both x and y are loop variables for the while loop, which illustrates one major difference between while loops and for loops. In a for loop, the loop variable is initialized and reassigned automatically by the Python interpreter to each element of the collection being looped over. In a while loop, the loop variable(s) must be initialized and reassigned explicitly in code that we write.

This difference makes while loops more flexible than for loops, as the programmer has full control over exactly how the loop variable changes. This is both a strength and a weakness! While loops can be used to express algorithms that are cumbersome or impossible to express with for loops, but at the cost of requiring the programmer to write more code to keep track of loop variables. Remember: the more code you write, the more potential there is for error. So a good rule of thumb is to use for loops where possible (when you have an explicit collection to loop over), and reserve while loops for situations that can’t be easily implemented with a for loop.

Parallel assignment

One subtlety of our loop body is the order in which the loop variables are updated. Suppose we had swapped the last two lines of the loop body:

    while y != 0:
        r = x % y
        y = r
        x = y

This is a really easy change to make, but also incorrect: because the statement y = r is executed first, the next statement x = y assigns x to the new value of y rather than its old one!

When performing reassignment of multiple variables, where the new variable values depend on the old ones, it is important to keep track of the reassignment order so that you don’t accidentally lose previous variable values. To avoid this problem altogether, Python has a neat feature called parallel assignment, in which multiple variables can be assigned in the same statment.

Here is how we can rewrite the loop body using parallel assignment:

    while y != 0:
        r = x % y
        x, y = y, r

The assignment statement x, y = y, r is evaluated as follows:

1. First, the right-hand side y, r is evaluated, producing two objects.
2. Then, each object is assigned to the corresponding variable on the left-hand side.

In parallel assignment, the right-hand side is fully evaluated before any variable reassignment occurs. This means that the assignment statement x, y = y, r has the same effect as y, x = r, y—order doesn’t matter, and so we can think of each variable assignment happening in parallel, without one affecting the other.

Parallel assignment is a very useful tool when reassigning variables, so please take advantage of it to help simplify your code and avoid the “update order” problem of variable reassignment. Here is how we can rewrite the euclidean_gcd using parallel assignment:

def euclidean_gcd(a: int, b: int) -> int:
    """Return the gcd of a and b.

    Preconditions:
    - a >= 0
    - b >= 0
    """
    # Step 1: initialize x and y
    x, y = a, b

    while y != 0:  # Step 4: repeat Steps 2 and 3 until y is 0
        # Step 2: calculate the remainder of x divided by y
        r = x % y

        # Step 3: reassign x and y
        x, y = y, r

    # Step 5: x now refers to the gcd of a and b
    return x

Documenting loop properties: loop invariants

Our implementation of euclidean_gcd doesn’t follow a typical pattern of code we’ve seen so far. If we didn’t know anything about the algorithm and were simply looking at the code, it would be quite mysterious why it works. To improve the readability of this code, we want some way of documenting what we know about the loop variables x and y inside the loop body.

Recall that the Euclidean Algorithm relies on one key property, that gcd(x, y) == gcd(y, x % y). At each loop iteration, x and y are updated so that x = y and y = x % y. The key property that we want to capture is that even though x and y change, their gcd doesn’t. Since x and y are initialized to a and b, another way to express this is that at every loop iteration, gcd(x, y) == gcd(a, b). We call this statement a loop invariant, which is a property about loop variables that must be true at the start and end of each loop iteration. This is similar to representation invariants, which are properties of instance attributes that must be true for every instance of a given data class.

By convention, we document loop invariants at the top of a loop body using an assert statement.

def euclidean_gcd(a: int, b: int) -> int:
    """Return the gcd of a and b.

    Preconditions:
    - a >= 0
    - b >= 0
    """
    # Step 1: initialize x and y
    x, y = a, b

    while y != 0:  # Step 4: repeat Steps 2 and 3 until y is 0
        assert math.gcd(x, y) == math.gcd(a, b)  # (NEW) Loop invariant

        # Step 2: calculate the remainder of x divided by y
        r = x % y

        # Step 3: reassign x and y
        x, y = y, r

    # Step 5: x now refers to the gcd of a and b
    return x

Because this loop invariant must be true at the start and end of each loop iteration, it is also true after the loop stops (i.e., when y == 0). In this case, the loop invariant tells us that gcd(x, 0) == gcd(a, b), and so we know that x == gcd(a, b), which is why x is returned.

Loop invariants are a powerful way to document properties of our code, to better enable us to reason about our code. But remember that loop invariants by themselves are just statements; the only way to know for sure whether a loop invariant is correct is to do a proof, much like the one we did at the beginning of this section.

The Extended Euclidean Algorithm

Recall the GCD characterization theorem from the previous section, which said that \(\gcd(m, n)\) is the smallest positive integer that can be written as a linear combination of \(m\) and \(n\). One thing that’s a bit unsatisfying about this theorem is that even though it says that \(\gcd(m, n)\) is a linear combination of \(m\) and \(n\), it does not tell us how to actually find that linear combination.

For example, we know \(\gcd(13, 10) = 1\). So by the GCD characterization theorem, there must exist integers \(p\) and \(q\) such that \(1 = 13p + 10q\). But what are \(p\) and \(q\)?You might spend some time “solving” this equation now, but the point of this example is to illustrate the difficulty of doing this in general!

It turns out that a variation of the Euclidean Algorithm called the Extended Euclidean Algorithm can be used to not just find the gcd of two numbers, but find the coefficients of the linear combination that equals that gcd. Here is the function specification for such an algorithm:

def extended_euclidean_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return the gcd of a and b, and integers p and q such that

    gcd(a, b) == p * a + b * q.

    Preconditions:
    - a >= 0
    - b >= 0

    >>> extended_euclidean_gcd(13, 10)
    (1, 7, -9)
    >>> 1 == 13 * 7 + 10 * -9  # This illustrates the linear combination
    True
    """

To get a sense of how we can implement this function, let’s first copy over the code for the Euclidean Algorithm from earlier, with a few of the comments removed

def extended_euclidean_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return the gcd of a and b, and integers p and q such that

    gcd(a, b) == p * a + b * q.

    Preconditions:
    - a >= 0
    - b >= 0

    >>> extended_euclidean_gcd(13, 10)
    (1, 7, -9)
    """
    x, y = a, b

    while y != 0:
        assert math.gcd(x, y) == math.gcd(a, b)  # Loop invariant

        r = x % y

        x, y = y, r

    return x, ..., ...  # Need to return "p" and "q" here!

In our implementation of the Euclidean Algorithm, each loop iteration makes the loop variables x and y smaller, while preserving the property gcd(x, y) == gcd(a, b). The key mathematical insight is that x and y are always linear combinations of a and b, at every loop iteration! This might sound surprising, so let’s double-check this informally:

• At the start, x = a, and a is certainly a combination of a and b: a == a * 1 + b * 0.
• At the start, y = b, and similarly b == a * 0 + b * 1.
• Inside the loop, x gets reassigned to y, and y gets reassigned to r = x % y.
  • If we assume y starts each loop iteration as a linear combination of a and b, then x will end each loop iteration as a lineary combination as well.
  • But what about x % y? From the definition of remainder, we know that x % y == x - q * y for some integer q (the quotient x // y in Python syntax). So this tells us x % y is a linear combination of x and y, and if we assume that x and y are both linear combinations of a and b, we can conclude that x % y is also a linear combination of a and b.A good exercise is to prove: for all integers \(a, b, c, d\), if \(c\) and \(d\) are both linear combinations of \(a\) and \(b\), then any linear combination of \(c\) and \(d\) is also a linear combination of \(a\) and \(b\).

This tells us that we can add the following two (English) loop invariants to our Extended Euclidean Algorithm:

def extended_euclidean_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return the gcd of a and b, and integers p and q such that

    gcd(a, b) == p * a + b * q.

    Preconditions:
    - a >= 0
    - b >= 0

    >>> extended_euclidean_gcd(13, 10)
    (1, 7, -9)
    """
    x, y = a, b

    while y != 0:
        assert math.gcd(x, y) == math.gcd(a, b)  # Loop invariant 1
        # Loop invariant 2: x is a linear combination of a and b
        # Loop invariant 3: y is a linear combination of a and b

        r = x % y

        x, y = y, r

    return x, ..., ...  # Need to return "p" and "q" here!

Okay, but it seems like we still have the same problem as before: just because we know that x and y are linear combinations of a and b doesn’t tell us how to find those combinations!

We need one more insight to make the Extended Euclidean Algorithm work: **we add four new loop variables px, qx, py, qy to represent the coefficients in the linear combinations for x and y. To make that concrete, here is yet another updated version of our code, with these four new loop variables. Pay special attention to our modified loop invariants 2 and 3, which now can be expressed in pure Python code.

def extended_euclidean_gcd(a: int, b: int) -> tuple[int, int, int]:
    """Return the gcd of a and b, and integers p and q such that

    gcd(a, b) == p * a + b * q.

    Preconditions:
    - a >= 0
    - b >= 0

    >>> extended_euclidean_gcd(13, 10)
    (1, 7, -9)
    """
    x, y = a, b

    # NEW: more loop variables
    px, qx = ..., ...
    py, qy = ..., ...

    while y != 0:
        assert math.gcd(x, y) == math.gcd(a, b)  # Loop invariant 1
        assert x == px * a + qx * b              # Loop invariant 2
        assert y == py * a + qy * b              # Loop invariant 3

        r = x % y

        x, y = y, r

        # NEW: update the new loop variables
        px, qx, py, qy = ..., ..., ..., ...

    return x, px, qx

With the introduction of these variables, we are now able to “fill in” the return statement: we return (x, px, qx), and Loop Invariant 2 tells us that x == px * a + qx * b.

So how do we fill in the remaining ...? Believe it or not, all of the required ideas were covered in our earlier “informal justification” that x and y are always linear combinations of a and b. We encourage you to think about this and try to complete this algorithm yourself, and bring your ideas to lecture when we discuss this algorithm. As a possible hint, here is a loop tracing table that demonstrates the correct values of all six (!) loop variables when calling extended_euclidean_gcd(100, 13):