gcd package:sbv
Generalized GCD, working for all integers. We simply call
nGCD with the absolute value of the arguments.
Symbolic GCD as our specification. Note that we cannot really
implement the GCD function since it is not symbolically terminating.
So, we instead uninterpret and axiomatize it below.
NB. The concrete part of the definition is only used in calls to
traceExecution and is not needed for the proof. If you don't
need to call
traceExecution, you can simply ignore that part
and directly uninterpret. In that case, we simply use Prelude's
version.
Computing GCD symbolically, and generating C code for it. This example
illustrates symbolic termination related issues when programming with
SBV, when the termination of a recursive algorithm crucially depends
on the value of a symbolic variable. The technique we use is to
statically enforce termination by using a recursion depth counter.
We define three different versions of the GCD algorithm: (1) Regular
version using the modulus operator, (2) the more basic version using
subtraction, and (3) the so called binary GCD. We prove that the
modulus based algorithm correct, i.e., that it calculates the
greatest-common-divisor of its arguments. We then prove that the other
two variants are equivalent to this version, thus establishing their
correctness as well.
Proof of correctness of an imperative GCD (greatest-common divisor)
algorithm, using weakest preconditions. The termination measure here
illustrates the use of lexicographic ordering. Also, since symbolic
version of GCD is not symbolically terminating, this is another
example of using uninterpreted functions and axioms as one writes
specifications for WP proofs.
Proof
>>> runTP gcdAdd
Lemma: dvdSum1 Q.E.D.
Lemma: dvdSum2 Q.E.D.
Lemma: gcdDivides Q.E.D.
Lemma: gcdLargest Q.E.D.
Lemma: gcdAdd
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Step: 4 Q.E.D.
Step: 5 Q.E.D.
Step: 6 Q.E.D.
Step: 7 Q.E.D.
Result: Q.E.D.
[Proven] gcdAdd :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Generalized version that works on arbitrary integers.
Instead of proving
gcdBin correct, we'll simply show that it
is equivalent to
gcd, hence it has all the properties we
already established.
Proof
>>> runTP gcdBinEquiv
Lemma: gcdEvenEven Q.E.D.
Lemma: gcdOddEven Q.E.D.
Lemma: gcdAdd Q.E.D.
Lemma: commutative Q.E.D.
Inductive lemma (strong): nGCDBinEquiv
Step: Measure is non-negative Q.E.D.
Step: 1 (5 way case split)
Step: 1.1 Q.E.D.
Step: 1.2 Q.E.D.
Step: 1.3.1 Q.E.D.
Step: 1.3.2 Q.E.D.
Step: 1.3.3 Q.E.D.
Step: 1.4.1 Q.E.D.
Step: 1.4.2 Q.E.D.
Step: 1.4.3 Q.E.D.
Step: 1.5 (3 way case split)
Step: 1.5.1 Q.E.D.
Step: 1.5.2.1 Q.E.D.
Step: 1.5.2.2 Q.E.D.
Step: 1.5.2.3 Q.E.D.
Step: 1.5.2.4 Q.E.D.
Step: 1.5.2.5 Q.E.D.
Step: 1.5.2.6 Q.E.D.
Step: 1.5.3.1 Q.E.D.
Step: 1.5.3.2 Q.E.D.
Step: 1.5.3.3 Q.E.D.
Step: 1.5.3.4 Q.E.D.
Step: 1.5.Completeness Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdBinEquiv
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Result: Q.E.D.
[Proven] gcdBinEquiv :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Putting it all together: GCD divides both arguments, and its maximal.
Proof
>>> runTP gcdCorrect
Lemma: gcdDivides Q.E.D.
Lemma: gcdMaximal Q.E.D.
Lemma: gcdCorrect
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Result: Q.E.D.
[Proven] gcdCorrect :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
GCD of two numbers divide these numbers. This is part one of the
proof, where we are not concerned with maximality. Our goal is to show
that the calculated gcd divides both inputs.
Proof
>>> runTP gcdDivides
Lemma: dvdAbs Q.E.D.
Lemma: helper
Step: 1 Q.E.D.
Result: Q.E.D.
Inductive lemma (strong): dvdNGCD
Step: Measure is non-negative Q.E.D.
Step: 1 (2 way case split)
Step: 1.1 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdDivides Q.E.D.
[Proven] gcdDivides :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Proof
>>> runTP gcdEvenEven
Lemma: modEE Q.E.D.
Inductive lemma (strong): nGCDEvenEven
Step: Measure is non-negative Q.E.D.
Step: 1 (2 way case split)
Step: 1.1 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.2.3 Q.E.D.
Step: 1.2.4 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdEvenEven
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Step: 4 Q.E.D.
Result: Q.E.D.
[Proven] gcdEvenEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Additionally prove that GCD is really maximum, i.e., it is the largest
in the regular sense. Note that we have to make an exception for
gcd 0 0 since by definition the GCD is
0, which is
clearly not the largest divisor of
0 and
0. (Since
any number is a GCD for the pair
(0, 0), there is no
maximum.)
Proof
>>> runTP gcdLargest
Lemma: gcdMaximal Q.E.D.
Lemma: gcdZero Q.E.D.
Lemma: gcdNonNegative Q.E.D.
Lemma: gcdLargest
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Result: Q.E.D.
[Proven] gcdLargest :: Ɐa ∷ Integer → Ɐb ∷ Integer → Ɐx ∷ Integer → Bool
Maximality. Any divisor of the inputs divides the GCD.
Proof
>>> runTP gcdMaximal
Lemma: dvdAbs Q.E.D.
Lemma: eDiv Q.E.D.
Lemma: helper
Step: 1 (x `dvd` a && x `dvd` b) Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Result: Q.E.D.
Inductive lemma (strong): mNGCD
Step: Measure is non-negative Q.E.D.
Step: 1 (2 way case split)
Step: 1.1 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdMaximal
Step: 1 (2 way case split)
Step: 1.1.1 Q.E.D.
Step: 1.1.2 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
[Proven] gcdMaximal :: Ɐa ∷ Integer → Ɐb ∷ Integer → Ɐx ∷ Integer → Bool
Proof
>>> runTP gcdNonNegative
Inductive lemma (strong): nonNegativeNGCD
Step: Measure is non-negative Q.E.D.
Step: 1 (2 way case split)
Step: 1.1 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: nonNegative Q.E.D.
[Proven] nonNegative :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Proof
>>> runTP gcdOddEven
Lemma: gcdDivides Q.E.D.
Lemma: gcdLargest Q.E.D.
Lemma: dvdMul Q.E.D.
Lemma: dvdOddThenOdd Q.E.D.
Lemma: dvdEvenWhenOdd Q.E.D.
Lemma: gcdOddEven
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Step: 4 Q.E.D.
Step: 5 Q.E.D.
Step: 6 Q.E.D.
Step: 7 Q.E.D.
Step: 8 Q.E.D.
Result: Q.E.D.
[Proven] gcdOddEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Generalized version of subtraction based GCD, working over all
integers.
Instead of proving
gcdSub correct, we'll simply show that it
is equivalent to
gcd, hence it has all the properties we
already established.
Proof
>>> runTP gcdSubEquiv
Lemma: commutative Q.E.D.
Lemma: gcdAdd Q.E.D.
Inductive lemma (strong): nGCDSubEquiv
Step: Measure is non-negative Q.E.D.
Step: 1 (5 way case split)
Step: 1.1 Q.E.D.
Step: 1.2 Q.E.D.
Step: 1.3 Q.E.D.
Step: 1.4.1 Q.E.D.
Step: 1.4.2 Q.E.D.
Step: 1.4.3 Q.E.D.
Step: 1.5.1 Q.E.D.
Step: 1.5.2 Q.E.D.
Step: 1.5.3 Q.E.D.
Step: 1.5.4 Q.E.D.
Step: 1.5.5 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdSubEquiv
Step: 1 Q.E.D.
Step: 2 Q.E.D.
Step: 3 Q.E.D.
Result: Q.E.D.
[Proven] gcdSubEquiv :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
Proof
>>> runTP gcdZero
Inductive lemma (strong): nGCDZero
Step: Measure is non-negative Q.E.D.
Step: 1 (2 way case split)
Step: 1.1 Q.E.D.
Step: 1.2.1 Q.E.D.
Step: 1.2.2 Q.E.D.
Step: 1.Completeness Q.E.D.
Result: Q.E.D.
Lemma: gcdZero Q.E.D.
[Proven] gcdZero :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
The state for the sum program, parameterized over a base type
a.
This call will generate the required C files. The following is the
function body generated for
sgcd. (We are not showing the
generated header,
Makefile, and the driver programs for
brevity.) Note that the generated function is a constant time
algorithm for GCD. It is not necessarily fastest, but it will take
precisely the same amount of time for all values of
x and
y.
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdbool.h>
#include "sgcd.h"
SWord8 sgcd(const SWord8 x, const SWord8 y)
{
const SWord8 s0 = x;
const SWord8 s1 = y;
const SBool s3 = s1 == 0;
const SWord8 s4 = (s1 == 0) ? s0 : (s0 % s1);
const SWord8 s5 = s3 ? s0 : s4;
const SBool s6 = 0 == s5;
const SWord8 s7 = (s5 == 0) ? s1 : (s1 % s5);
const SWord8 s8 = s6 ? s1 : s7;
const SBool s9 = 0 == s8;
const SWord8 s10 = (s8 == 0) ? s5 : (s5 % s8);
const SWord8 s11 = s9 ? s5 : s10;
const SBool s12 = 0 == s11;
const SWord8 s13 = (s11 == 0) ? s8 : (s8 % s11);
const SWord8 s14 = s12 ? s8 : s13;
const SBool s15 = 0 == s14;
const SWord8 s16 = (s14 == 0) ? s11 : (s11 % s14);
const SWord8 s17 = s15 ? s11 : s16;
const SBool s18 = 0 == s17;
const SWord8 s19 = (s17 == 0) ? s14 : (s14 % s17);
const SWord8 s20 = s18 ? s14 : s19;
const SBool s21 = 0 == s20;
const SWord8 s22 = (s20 == 0) ? s17 : (s17 % s20);
const SWord8 s23 = s21 ? s17 : s22;
const SBool s24 = 0 == s23;
const SWord8 s25 = (s23 == 0) ? s20 : (s20 % s23);
const SWord8 s26 = s24 ? s20 : s25;
const SBool s27 = 0 == s26;
const SWord8 s28 = (s26 == 0) ? s23 : (s23 % s26);
const SWord8 s29 = s27 ? s23 : s28;
const SBool s30 = 0 == s29;
const SWord8 s31 = (s29 == 0) ? s26 : (s26 % s29);
const SWord8 s32 = s30 ? s26 : s31;
const SBool s33 = 0 == s32;
const SWord8 s34 = (s32 == 0) ? s29 : (s29 % s32);
const SWord8 s35 = s33 ? s29 : s34;
const SBool s36 = 0 == s35;
const SWord8 s37 = s36 ? s32 : s35;
const SWord8 s38 = s33 ? s29 : s37;
const SWord8 s39 = s30 ? s26 : s38;
const SWord8 s40 = s27 ? s23 : s39;
const SWord8 s41 = s24 ? s20 : s40;
const SWord8 s42 = s21 ? s17 : s41;
const SWord8 s43 = s18 ? s14 : s42;
const SWord8 s44 = s15 ? s11 : s43;
const SWord8 s45 = s12 ? s8 : s44;
const SWord8 s46 = s9 ? s5 : s45;
const SWord8 s47 = s6 ? s1 : s46;
const SWord8 s48 = s3 ? s0 : s47;
return s48;
}
The symbolic GCD algorithm, over two 8-bit numbers. We define sgcd
a 0 to be a for all a, which implies sgcd 0
0 = 0. Note that this is essentially Euclid's algorithm, except
with a recursion depth counter. We need the depth counter since the
algorithm is not symbolically terminating, as we don't have a
means of determining that the second argument (b) will
eventually reach 0 in a symbolic context. Hence we stop after 12
iterations. Why 12? We've empirically determined that this algorithm
will recurse at most 12 times for arbitrary 8-bit numbers. Of course,
this is a claim that we shall prove below.
We have:
>>> prove sgcdIsCorrect
Q.E.D.
nGCD is the version of GCD that works on non-negative
integers.
Ideally, we should make this function local to gcd, but then
we can't refer to it explicitly in our proofs.
Note on maximality: Note that, by definition gcd 0 0 = 0.
Since any number divides 0, there is no greatest common
divisor for the pair (0, 0). So, maximality here is meant to
be in terms of divisibility. That is, any divisor of a and
b will also divide their gcd.
nGCDBin is the binary GCD algorithm that works on
non-negative numbers.