Proving that euclidian algorithm complexity is at most logarithmic

I’m helping lately with C++ teaching for students with some mathematical background; obviously, you can’t avoid the topic of recursion at some point and it comes with the stereotypical exercises about factorial, Fibonacci sequences and computation of the Greatest Common Denominator.

The students are asked to write a naive (eg not tail recursive) recursive Fibonacci algorithm and to witness that it triggers a stack overflow for index like 60. We then ask them to write the GCD algorithm and to explain why it is robust even with big input.

It boils down to the argument that every two step the size of the input is halved, that is if \(a\) and \(b\) are natural numbers then \(a \mod b \leq \frac{a}{2} \). The proof is quite simply and involves a case analysis so I thought it would be a nice opportunity to write a coq proof.

I’m using mathcomp and ssrnat type for this task. Manipulating equations in Coq is cumbersome since it involves issuing all intermediate steps including the most trivial ones like reassociating operators or commutating operands. That’s why vanilla Coq provides a couple of tactics that help proving that some equation hold, namely:

  • lia (previously omega) which stands for linear integer arithmetic. It tries to prove equations or inequalities where both left hand side and right hand side are of the form \(\sum C_i x_i \geq 0\) where the \(C_i\) are some integer constants.
  • lra which stands for linear rational arithmetic, which is the same, except that it works for real numbers.
  • nia which stands for nonlinear integer arithmetic. I’m not sure which kind of equations it covers, according to the doc 1 it uses lia under the hood.
  • nra which is the same as nia but for rational number.

These tactics are not available to mathcomp natively but fortunately latest Coq version (from 8.11 I think) provides a Zify module which allows to add support for custom type for these tactics. Kazuhiko Sakaguchi wrote a lib, mczify (for MathComp Zify), which add supports for at least ssrbool and ssrnat; this lib is not on opam yet but we can build it and add it to our opam installation quite easily:

  • First you need to clone the repository: https://github.com/math-comp/mczify with the branch corresponding to your coq version (coq-8.11 or coq-8.12 at the time of writing, master doesn’t work with coq 8.12) and add an opam file whose content is given below.
  • Then you can call « opam install . » to make it build mczify and install it in our opam folder.
opam-version: "2.0"
maintainer: "pi8027@gmail.com"
homepage: "https://github.com/math-comp/mczify"
bug-reports: "https://github.com/math-comp/mczify/issues"
dev-repo: "git+https://github.com/math-comp/mczify.git"
license: "CECILL-C"
authors: [
  "Kazuhiko Sakaguchi"
]
build: [
  [make "INSTMODE=global"]
  [make "-j%{jobs}%"]
]
install: [
  [make "install"]
]
depends: [
  "coq" { ((>= "8.10" &amp; < "8.13~") | = "dev") }
  "coq-mathcomp-field"       {(>= "1.11.0" &amp; < "1.12~")}
  "coq-mathcomp-finmap"      {(>= "1.5.0" &amp; < "1.6~")}
]
synopsis: "An analysis library for mathematical components"
description: """
This repository contains an experimental library for real analysis for
the Coq proof-assistant and using the Mathematical Components library.

It is inspired by the Coquelicot library.
"""
tags: [
  "category:Mathematics/Real Calculus and Topology"
  "keyword: analysis"
  "keyword: topology"
  "keyword: real numbers"
  "logpath: mathcomp.analysis"
  "date:2020-08-11"
]

Now we can check that we correctly installed zify by importing it as well as all the ssrnat and div libs:

From mathcomp Require Import ssreflect ssrnat div ssrbool.
From mathcomp Require Import zify.

Now we can start proving our proof script. Formulating the lemma isn’t complicated. I added some hypotheses (\( b \neq 0\) and \( a \ge b\)) that I found necessary while writing the proof. They really don’t change the scope of what we’re trying to prove though.

Lemma half: forall a b : nat, 0 < b -> a >= b -> a %% b <= a %/ 2.

Now the proof. We will follow the proof by case here. The first step is to move all our hypotheses from the goal to the context:

move => a b H0 H1.

And then do the case analysis. While this may sound trivial in classical logic, we’re working with constructive mathematic, where the excluded middle lemma doesn’t hold in general. This is however the case for our example : the comparison function on nat (\(\leq\)) is actually a decision procedure, meaning that there is an algorithm that can compute a result for \(a \leq b\) where \(a\) and \(b\) are natural number. Note that this isn’t true for every type: for instance, it is known2 that equality on the real type isn’t decidable and so inequality isn’t either.

This means that the case tactic requires us to provide a proof that there is actually two cases (and thus as the expression is decidable). For reference in vanilla Coq one would look for something of the form :

forall (a b:nat), {a <= b} + {a > b}

In mathcomp however the lemma is a little less obvious : it’s leqP whose type is :

Lemma leqP m n : leq_xor_gtn m n (minn n m) (minn m n) (maxn n m) (maxn m n) (m <= n) (n < m).

It’s a common trick used in several place of mathcomp where some lemma use an intermediate inductive type which embeds various consequences of a case analysis. Here we see what’s happening if we look at what leq_xor_gtn embeds :

Variant leq_xor_gtn m n : nat -> nat -> nat -> nat -> bool -> bool -> Set :=
  | LeqNotGtn of m <= n : leq_xor_gtn m n m m n n true false
  | GtnNotLeq of n < m  : leq_xor_gtn m n n n m m false true.

This reads as : there is two case (LeqNotGtn and GtnNotLeq), and for instance in the first case (min n m) and (min m n) is m, (max n m) and (max m n) is n, (m <= n) is true and (m < n) is false.

So now we can do our case analysis and put the case in the context :

case (leqP b (a %/2)) => H.

First case is trivial: it’s \(b \le \frac{a}{2}\). We just need to notice that \(a \mod b \le b\) which is provided by ltn_pmod lemma :

    apply: leq_trans.
    apply: ltnW.
    apply: ltn_pmod.
    by []. by [].

Second case is a little more involving: we will use the fact that the quotient of \(a\) per \(b\) is at least \(1\) so \(a\mod b \leq a – b\) which in turns is less than \(\frac{a}{2}\) due to the hypothesis.

First let’s prove that \(a / b\leq b\), it relies on the lemma leq_divRL which proves : (m <= n %/ d) = (m * d <= n).

    have: 1 <= a %/b.
        rewrite leq_divRL.
        by rewrite mul1n. by [].

Now let’s transform the goal into \( a – (a / b) * b \leq \frac{a}{2}\) ; for that we will use the lemma divn_eq (which expand \(a\) into \(a = (a / b) * b + a \mod b\)) but it’s not enough : we want to substitute \(a \mod b \) so we need to have an expression of the form \( a\mod b = \cdots\). This can be tedious to do manually and that’s where the tactic presented earlier comes in handy. Since it’s just substracting a term on both side of the equation, lia is enough.

    have: a %% b = a -  a %/ b * b.
        move: (divn_eq a b).
        lia.
    move ->.

Now all the elements are in place. An opportunistic call to the tactic nia (since this time there is a little more than moving integers expression between sides of an equation) finishes the proof:

    nia.

For reference here is the complete script proof:

From mathcomp Require Import ssreflect ssrnat div ssrbool.
From mathcomp Require Import zify.

Lemma half: forall a b : nat, 0 < b -> a >= b -> a %% b <= a %/ 2.
Proof.
    move => a b H0 H1.
    case (leqP b (a %/2)) => H.
        apply: leq_trans.
        apply: ltnW.
        apply: ltn_pmod.
        by []. by [].

    have: 1 <= a %/b.
        Check leq_divRL.
        rewrite leq_divRL.
        by rewrite mul1n. by [].

    have: a %% b = a -  a %/ b * b.
        move: (divn_eq a b).
        lia.
    move ->.
    nia.
Qed.
  1. https://coq.inria.fr/refman/addendum/micromega.html#coq:tacn.nia
  2. https://en.wikipedia.org/wiki/Richardson%27s_theorem

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *