GNATprove Tips and Tricks: Proving the Ghost Common Divisor (GCD)

by Yannick Moy in Formal Verification – April 6, 2017

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:

   function Divides (A, B : in Positive) return Boolean is (B mod A = 0) with Ghost;

   function GCD (A, B : in Positive) return Positive with
     Post => Divides (GCD'Result, A)
       and then Divides (GCD'Result, B);

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:

   function GCD (A, B : in Positive) return Positive is
   begin
      return 1;
   end 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:

   function GCD (A, B : in Positive) return Positive with
     Post => Divides (GCD'Result, A)
       and then Divides (GCD'Result, B)
       and then (for all X in GCD'Result + 1 .. Integer'Min (A, B) =>
                   not (Divides (X, A) and Divides (X, B)));

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:

   function GCD (A, B : in Positive) return Positive is
      C : Positive := Integer'Min (A, B);
   begin
      while C > 1 loop
         exit when A mod C = 0 and B mod C = 0;
         pragma Loop_Invariant (for all X in C .. Integer'Min (A, B) =>
                                  not (Divides (X, A) and Divides (X, B)));
         C := C - 1;
      end loop;

      return C;
   end GCD;

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:

   function GCD (A, B : in Positive) return Positive is
      C : Positive := Integer'Min (A, B);
   begin
      if A mod C = 0 and B mod C = 0 then
         return C;

      else
         C := C / 2;

         while C > 1 loop
            exit when A mod C = 0 and B mod C = 0;
            pragma Loop_Invariant (for all X in C .. Integer'Min (A, B) =>
                                     not (Divides (X, A) and Divides (X, B)));
            C := C - 1;
         end loop;

         return C;
      end if;
   end GCD;

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:

math_simple_half.ads:9:20: medium: postcondition might fail, cannot prove not (Divides (X, A) and Divides (X, B)

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:

         ...
         C := C / 2;
         for J in C + 1 .. Integer'Min (A, B) - 1 loop
            pragma Loop_Invariant (for all X in C + 1 .. J =>
                                     not Divides (X, Integer'Min (A, B)));
         end loop;

         while C > 1 loop
         ...

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:

math_simple_half.adb:24:38: medium: loop invariant might fail after first iteration, cannot prove not Divides (X, Integer'min)

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:

   procedure Lemma_Not_Divisor (Arg1, Arg2 : Positive) with
     Ghost,
     Global => null,
     Pre  => Arg1 in Arg2 / 2 + 1 .. Arg2 - 1,
     Post => not Divides (Arg1, Arg2)
   is
   begin
      null;
   end Lemma_Not_Divisor;

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:

         ...
         C := C / 2;
         for J in C + 1 .. Integer'Min (A, B) - 1 loop
            Lemma_Not_Divisor (J, Integer'Min (A, B));
            pragma Loop_Invariant (for all X in C + 1 .. J =>
                                     not Divides (X, Integer'Min (A, B)));
         end loop;

         while C > 1 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:

math_simple_half.adb:32:20: medium: range check might fail

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):

   procedure Lemma_Not_Divisor (Arg1, Arg2 : Positive) with
     Ghost,
     Global => null,
     Pre  => Arg1 in Arg2 / 2 + 1 .. Arg2 - 1,
     Post => not Divides (Arg1, Arg2)
   is
   begin
      null;
   end Lemma_Not_Divisor;

   function GCD (A, B : in Positive) return Positive is
      C : Positive := Integer'Min (A, B);
   begin
      if A mod C = 0 and B mod C = 0 then
         return C;

      else
         C := C / 2;
         for J in C + 1 .. Integer'Min (A, B) - 1 loop
            Lemma_Not_Divisor (J, Integer'Min (A, B));
            pragma Loop_Invariant (for all X in C + 1 .. J =>
                                     not Divides (X, Integer'Min (A, B)));
         end loop;

         while C > 1 loop
            exit when A mod C = 0 and B mod C = 0;
            pragma Loop_Invariant (C > 1);
            pragma Loop_Invariant (for all X in C .. Integer'Min (A, B) =>
                                     not (Divides (X, A) and Divides (X, B)));
            C := C - 1;
         end loop;

         return C;
      end if;
   end GCD;

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

   function GCD (A, B : in Positive) return Positive is
      An : Positive := A;
      Bn : Natural := B;
      C  : Positive;
   begin
      while Bn /= 0 loop
         C  := An;
         An := Bn;
         Bn := C mod Bn;
      end loop;

      return An;
   end 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:

         pragma Loop_Invariant
           (for all X in Positive =>
              (Divides (X, A) and Divides (X, B))
                =
              (Divides (X, An) and (Bn = 0 or else Divides (X, Bn))));

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:

   procedure Lemma_Same_Divisor_Mod (A, B : Positive) with
     Ghost,
     Global => null,
     Post => (for all X in Positive =>
                (Divides (X, A) and Divides (X, B))
                  =
                (Divides (X, B) and then (Divides (B, A) or else Divides (X, A mod B))))

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:

   procedure Lemma_Divisor_Mod (A, B, X : Positive) with
     Ghost,
     Global => null,
     Pre  => Divides (X, A) and then Divides (X, B) and then not Divides (B, A),
     Post => Divides (X, A mod B)
   is
   begin
      null;
   end Lemma_Divisor_Mod;

   procedure Lemma_Divisor_Transitive (A, B, X : Positive) with
     Ghost,
     Global => null,
     Pre  => Divides (B, A) and then Divides (X, B),
     Post => Divides (X, A)
   is
   begin
      null;
   end Lemma_Divisor_Transitive;

   procedure Lemma_Divisor_Mod_Inverse (A, B, X : Positive) with
     Ghost,
     Global => null,
     Pre  => not Divides (B, A) and then Divides (X, A mod B) and then Divides (X, B),
     Post => Divides (X, A)
   is
   begin
      null;
   end Lemma_Divisor_Mod_Inverse;

   procedure Lemma_Same_Divisor_Mod (A, B : Positive) with
     Ghost,
     Global => null,
     Post => (for all X in Positive =>
                (Divides (X, A) and Divides (X, B))
                  =
                (Divides (X, B) and then (Divides (B, A) or else Divides (X, A mod B))))
   is
   begin
      for X in Positive loop
         if Divides (X, A) and then Divides (X, B) and then not Divides (B, A) then
            Lemma_Divisor_Mod (A, B, X);
         end if;
         if Divides (B, A) and then Divides (X, B) then
            Lemma_Divisor_Transitive (A, B, X);
         end if;
         if not Divides (B, A) and then Divides (X, A mod B) and then Divides (X, B) then
            Lemma_Divisor_Mod_Inverse (A, B, X);
         end if;
         pragma Loop_Invariant
           (for all Y in 1 .. X =>
              (Divides (Y, A) and Divides (Y, B))
                =
              (Divides (Y, B) and then (Divides (B, A) or else Divides (Y, A mod B))));
      end loop;
   end Lemma_Same_Divisor_Mod;

   function GCD (A, B : in Positive) return Positive is
      An : Positive := A;
      Bn : Natural := B;
      C  : Positive;
   begin
      while Bn /= 0 loop
         C  := An;
         An := Bn;
         Bn := C mod Bn;
         Lemma_Same_Divisor_Mod (C, An);
         pragma Loop_Invariant
           (for all X in Positive =>
              (Divides (X, A) and Divides (X, B))
                =
              (Divides (X, An) and (Bn = 0 or else Divides (X, Bn))));
      end loop;

      pragma Assert (Divides (An, An));
      pragma Assert (Divides (An, A));
      pragma Assert (Divides (An, B));

      pragma Assert (for all X in An + 1 .. Integer'Min (A, B) =>
                       not (Divides (X, A) and Divides (X, B)));
      return An;
   end GCD;

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:

   function Divides (A, B : in Positive) return Boolean is
     (for some C in Positive => A * C = B)
   with Ghost;

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).

comments powered by Disqus