Verifying the binary algorithm for greatest common divisors

22 Feb 2023

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

This system of rules corresponds precisely to an Isabelle/HOL inductive definition. Note that we are defining a 3-argument relation rather than a 2-argument 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’=b-1$, 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.

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