Euclid's algorithm for computing the greatest common divisor of two numbers is one of the first ones we learn in school (source: myself), and also one of the first algorithms that humans devised (source: WikiPedia). So it's quite appealing to try to prove it with an automatic proving toolset like SPARK. It turns out that proving it automatically is not so easy, just like understanding why it works is not so easy. To prove it, I am going to use a fair amount of ghost code, to convey the necessary information to GNATprove, the SPARK proof tool.

Let's start with the specification of what the algorithm should do: computing the greatest common divisor of two positive numbers. The algorithm can be described for natural numbers or even negative numbers, but we'll keep things simple by using positive inputs. Given that Euclid himself only considered positive numbers, that should be enough for us here. A common divisor C of two integers A and B should divide both of them, which we usually express efficiently in software as (A mod C = 0) and (B mod C = 0). This gives the following contract for GCD:

Note that we're already using some ghost code here: the function Divides is marked as Ghost because it is only used in specifications. The problem with the contract above is that it is satisfied by a very simple implementation that does not actually compute the GCD:

And GNATprove manages to prove that easily. So we need to specify the GCD function more precisely, by stating that it returns not only a common divisor of A and B, but the greatest of them:

Here, I am relying on the property that the GCD is lower than both A and B, to only look for a greater divisor between the value returned (GCD'Result) and the minimum of A and B (Integer'Min (A, B)). That contract uniquely defines the GCD, and indeed GNATprove cannot prove that the previous incorrect implementation of GCD implements this contract.

Let's implement now the GCD, starting with a very simple linear search that looks for the greatest common divisor of A and B:

We start at the minimum of A and B, and check if this value divides A and B. If so, we're done and we exit the loop. Otherwise we decrease this value by 1 and we continue. This is not efficient, but this computes the GCD as desired. In order to prove that implementation, we need a simple loop invariant that states that no common divisor of A and B has been encountered in the loop so far. GNATprove proves easily that this implementation implements the contract of GCD (with --level=2 -j0 in less than a second).

Can we increase the efficiency of the algorithm? If we had to do this computation by hand, we would try (Integer'Min (A, B) / 2) rather than (Integer'Min (A, B) - 1) as second possible value for the GCD, because we know that no value in between has a chance of dividing both A and B. Here is an implementation of that idea:

I kept the same loop invariant as before, because we need that loop invariant to prove the postcondition of GCD. But GNATprove does not manage to prove that this invariant holds the first time that it enters the loop:

Indeed, we know that the GCD cannot lie between (Integer'Min (A, B) / 2) and Integer'Min (A, B) in that case, but provers don't. We have to provide the necessary information for GNATprove to deduce this fact. This is where we are going to use ghost code more extensively. As we need to prove a quantified property, we are going to establish it progressively, though a loop, accumulating the knowledge gathered in a loop invariant:

With the help of this ghost code, the loop invariant in our original loop is now fully proved! But the loop invariant of the ghost loop itself is not proved, for iterations beyond the initial one:

The property we need here is that, given a value J between (Integer'Min (A, B) / 2 + 1) and (Integer'Min (A, B) - 1), J cannot divide Integer'Min (A, B). We are going to establish this property in a separate ghost lemma, that we'll call inside the loop to get this property. The lemma is an abstraction of the property just stated:

Even better, by isolating the desired property from its context of use, GNATprove can now prove it! Here, it is the underlying prover Z3 which manages to prove the postcondition of Lemma_Not_Divisor. It remains to call this lemma inside the loop:

We're almost done. The only check that GNATprove does not manage to prove is the range check on the decrement to C in our original loop:

It looks trivial though. The loop test states that (C > 1), so it should be easy to prove that C can be decremented and remain positive, right? Well, it helps here to understand how proof is carried on loops, by creating an artificial loop starting from the loop invariant in an iteration and ending in the same loop invariant at the next iteration. Because the loop invariant is not the first statement of the loop here, the loop test (C > 1) is not automatically added to the loop invariant (which would be incorrect in general, if the first statement modified C), so provers do not see this hypothesis when trying to prove the range check later. All that is needed here is to explicitly carry the loop test inside the loop invariant. Here is the final implementation of the more efficient GCD, which is proved easily by GNATprove (with --level=2 -j0 in less than two seconds):

It's time to try to prove Euclid's variant of the algorithm, still with the same contract on GCD:

The key insight here to justify the large jumps when decreasing the value of C is that, at each iteration of the loop, any common divisor between A and B is still a common divisor of An and Bn. And conversely. So the greatest common divisor of A and B turns out to be also the greatest common divisor of An and Bn. Which is just An when Bn turns out to be zero. Here is the loop invariant that expresses that property:

In order to prove this loop invariant, we need to state a ghost lemma like before, that helps prove the property needed at a particular iteration of the loop:

Because this lemma is more complex, it is not proved automatically like the previous one. Here, we need to implement it with ghost code, extracting three further lemmas that help proving this one. Here is the complete code for this implementation:

The implementation of GCD as well as the main lemma Lemma_Same_Divisor_Mod are proved easily by GNATprove (with --level=2 -j0 in less than ten seconds). The three elementary lemmas Lemma_Divisor_Mod, Lemma_Divisor_Transitive and Lemma_Divisor_Mod_Inverse are not proved automatically. Currently I only reviewed them and convinced myself that they were true. A better solution in the future will be to prove them in Coq, which should not be difficult given that there are similar existing lemmas in Coq, and include them in the SPARK lemma library.

All the code for these versions of GCD can be found in the GitHub repository of SPARK 2014, including another version of GCD where I changed the definition of Divides to reflect more accurately the mathematical notion of divisibility:

This is also proved automatically by GNATprove, and this proof also requires the introduction of suitable lemmas that relate the mathematical notion of divisibility to its programming counterpart (B mod A = 0).