Verifying the binary algorithm for greatest common divisors
[examples
gcd
Isar
]
The Euclidean algorithm for the GCD is generally regarded as one of the first algorithms worthy of the name –– its main competitor perhaps the sieve of Eratosthenes, for generating prime numbers. Euclid’s algorithm requires repeated calculations of integer remainders, which was a heavy operation on early computers. (Some provided no hardware to perform multiplication and division.) One can replace the remainder operation by subtraction. That algorithm is far too slow, especially if one of the arguments is much smaller than the other, but it can be refined to achieve high performance by noting that any binary computer can efficiently divide by two. Let’s verify this algorithm, and along the way see how inductive definitions are done in Isabelle/HOL.
The binary GCD algorithm
The greatest common divisor of two natural numbers satisfies the following properties:
 The GCD of $x$ and 0 is $x$.
 If the GCD of $x$ and $y$ is $z$, then the GCD of $2x$ and $2y$ is $2z$.
 The GCD of $2x$ and $y$ is the same as that of $x$ and $y$ if $y$ is odd.
 The GCD of $x$ and $y$ is the same as that of $xy$ and $y$ if $y\le x$.
 The GCD of $x$ and $y$ is the same as the GCD of $y$ and $x$.
This system of rules corresponds precisely to an Isabelle/HOL inductive definition. Note that we are defining a 3argument relation rather than a 2argument function. That’s because frequently more than one of these cases is applicable, so it is not immediately obvious that they express a function.
inductive_set bgcd :: "(nat × nat × nat) set" where bgcdZero: "(u, 0, u) ∈ bgcd"  bgcdEven: "⟦ (u, v, g) ∈ bgcd ⟧ ⟹ (2*u, 2*v, 2*g) ∈ bgcd"  bgcdOdd: "⟦ (u, v, g) ∈ bgcd; ¬ 2 dvd v ⟧ ⟹ (2*u, v, g) ∈ bgcd"  bgcdStep: "⟦ (u  v, v, g) ∈ bgcd; v ≤ u ⟧ ⟹ (u, v, g) ∈ bgcd"  bgcdSwap: "⟦ (v, u, g) ∈ bgcd ⟧ ⟹ (u, v, g) ∈ bgcd"
The inductive definition is abstract, but the corresponding algorithm isn’t hard to see. If one operand is zero, then return the other; if both operands are even, then divide by two, obtain the GCD recursively and then double it; if only one of the operands is even then divide it by 2; if both of the operands are odd, then subtract the smaller from the larger. Testing whether an integer is even or odd, and dividing it by two, is a trivial binary shift.
Proving that the algorithm is correct
We must show that the computed “result”
($g$ in the triple $(x,y,g)$) really is the greatest common divisor.
First we simply prove that it is a common divisor.
The proof is by induction on the relation bgcd
,
which means reasoning separately on each of the five rules shown above.
Note that four of the five cases are trivial enough to be covered
by a single call to auto
at the end, with one case (subtraction)
handled specifically.
lemma bgcd_divides: "(x,y,g) ∈ bgcd ⟹ g dvd x ∧ g dvd y" proof (induct rule: bgcd.induct) case (bgcdStep u v g) with dvd_diffD show ?case by blast qed auto
So, it yields a common divisor. Let’s show that it is the greatest. And here it’s time to stress that greatest refers to divisibility, not magnitude.
lemma bgcd_greatest: "(x,y,g) ∈ bgcd ⟹ d dvd x ⟹ d dvd y ⟹ d dvd g" proof (induct arbitrary: d rule: bgcd.induct) case (bgcdEven u v g d) show ?case proof (cases "2 dvd d") case True thus ?thesis using bgcdEven by (force simp add: dvd_def) next case False thus ?thesis using bgcdEven by (simp add: coprime_dvd_mult_right_iff) qed next case (bgcdOdd u v g d) hence "coprime d 2" by fastforce thus ?case using bgcdOdd by (simp add: coprime_dvd_mult_right_iff) qed auto
As with the previous proof, it’s by induction on the relation we defined,
and most of the cases are trivially solved by auto
.
The case where both $u$ and $v$ are even is verified with separate cases
depending on whether $d$ (which is some other common divisor)
is even or not.
Proving uniqueness and existence
The two results just proved are enough to show that the relation bgcd
is indeed a function, in the sense that it yields at most one result:
lemma bgcd_unique: "(x,y,g) ∈ bgcd ⟹ (x,y,g') ∈ bgcd ⟹ g = g'" by (meson bgcd_divides bgcd_greatest gcd_nat.strict_iff_not)
Our final task is to show that this function is total, that some result is always determined. (Strictly speaking, we don’t have an operational semantics, so this is not quite the termination of an actual algorithm.) The proof is by complete induction, using the observation that the GCD of any given pair of values is ultimately determined by that of a smaller pair of values, possibly with the help of the swap rule. So when considering the GCD of $a$ and $b$, the induction will be on their sum.
lemma bgcd_defined_aux: "a+b ≤ n ⟹ ∃ g. (a, b, g) ∈ bgcd" proof (induction n arbitrary: a b rule: less_induct) case (less n a b) show ?case proof (cases b) case 0 thus ?thesis by (metis bgcdZero) next case (Suc b') then have *: "a + b' < n" using Suc_le_eq add_Suc_right less.prems by presburger show ?thesis proof (cases "b ≤ a") case True thus ?thesis by (metis bgcd.simps le_add1 le_add_diff_inverse less.IH [OF *]) next case False then show ?thesis by (metis less.IH [OF *] Suc Suc_leI bgcd.simps le_add_diff_inverse less_add_same_cancel2 nle_le zero_less_iff_neq_zero) qed qed qed
In the proof above, we reason by cases on whether $b=0$ or alternatively $b = b’+1$ for some $b’$ (where obviously $b’=b1$, a fact that doesn’t need to be used). We note that $a+b’<n$, allowing use of the induction hypothesis. It’s reasonable to ask, why not just do mathematical induction on $b$? And the answer is, I couldn’t get a proof that way, but maybe you will.^{1}
lemma bgcd_defined: "∃!g. (a, b, g) ∈ bgcd" using bgcd_defined_aux bgcd_unique by auto
We have finally established that our inductive definition uniquely identifies a result for every pair of operands. And as noted above, that result is the GCD.
theorem bgcd_defined: "∃!g. (a, b, g) ∈ bgcd" using bgcd_defined_aux bgcd_unique by auto
This example is based on an assignment set in 2010 for my late, lamented MPhil course, Interactive Formal Verification. The Isabelle theory file can be downloaded.

See the comments below for a simplification suggested by YawarRaza7349. ↩