\documentclass[dvipdfmx]{book}
\newcommand{\VolumeName}{Volume 13: Proving Axiom Correct}
\usepackage{bbold}
\usepackage{amsmath}
\input{bookheader.tex}
\mainmatter
\setcounter{chapter}{0} % Chapter 1
Ultimately we would like Axiom to be able to prove that an
algorithm generates correct results. There are many steps
between here and that goal, including proving one Axiom
algorithm correct through all of the levels from Spad code,
to the Lisp code, to the C code, to the machine code; a 
daunting task of its own. 

The proof of a single Axiom algorithm is done with an eye toward
automating the process. Automated machine proofs are not possible
in general but will exist for known algorithms. 

\begin{quote}
{\bf Writing is nature's way of letting you know how sloppy your
thinking is
} -- Guindon\cite{Lamp02}
\end{quote}

\begin{quote}
{\bf Mathematics is nature's way of letting you know how sloppy
your writing is.
} -- Leslie Lamport\cite{Lamp02}
\end{quote}

\begin{quote}
{\bf 
The existence of the computer is giving impetus to the discovery of
algorithms that generate proofs. I can still hear the echos of the
collective sigh of relief that greeted the announcement in 1970 that
there is no general algorithm to test for integer solutions to
polynomial Diophantine equations; Hilbert's tenth problem has no
solution. Yet, as I look at my own field, I see that creating
algorithms that generate proofs constitutes some of the most important
mathematics being done. The all-purpose proof machine may be dead, but
tightly targeted machines are thriving.}
-- Dave Bressoud \cite{Bres93}
\end{quote}

\begin{quote}
{\bf In contrast to humans, computers are good at performing formal
processes. There are people working hard on the project of actually
formalizing parts of mathematics by computer, with actual formally
correct formal deductions. I think this is a very big but very
worthwhile project, and I am confident that we will learn a lot from
it. The process will help simplify and clarify mathematics. In not too
many years, I expect that we will have interactive computer programs
that can help people compile significant chunks of formally complete
and correct mathematics (based on a few perhaps shaky but at least
explicit assumptions) and that they will become part of the standard
mathematicians's working environment.}
-- William P. Thurston \cite{Thur94}
\end{quote}

\begin{quote}
{\bf Our basic premise is that the ability to construct and modify programs
will not improve without a new and comprehensive look at the entire
programming process. Past theoretical research, say, in the logic of
programs, has tended to focus on methods for reasoning about
individual programs; little has been done, it seems to us, to develop
a sound understanding of the process of programming -- the process by
which programs evolve in concept and in practice. At present, we lack
the means to describe the techniques of program construction and
improvement in ways that properly link verification, documentation and
adaptability.}

-- Scherlis and Scott (1983) in \cite{Maso86}
\end{quote}

\begin{quote}
{\bf ...constructive mathematics provides a way of viewing the language of
logical propositions as a {\sl specification} language for
programs. An ongoing thrust of work in computer science has been to
develop program specification languages and formalisms for
systematically deriving programs from specifications. For constructive
mathematics to provide such a methodology, techniques are needed for
systematically extracting programs from constructive proofs. Early work
in this field includes that of Bishop and Constable. What
distinguished Martin-L\"of's '82 type theory was that the method it
suggested for program synthesis was exceptionally simple: a direct
correspondence was set up between the constructs of mathematical
logic, and the constructs of a functional programming
language. Specifically, every proposition was considered to be
isomorphic to a type expression, and the proof of a proposition would
suggest precisely how to construct an inhabitant of the type, which
would be a term in a functional programming language. The term that
inhabits the type corresponding to a proposition is often referred to as
the {\sl computational content} of the proposition.}

-- Paul Bernard Jackson\cite{Jack95}

\end{quote}

\begin{quote}
Q: Why bother doing proofs about programming languages? They are almost
always boring if the definitions are right.

A: The definitions are almost always wrong.
\end{quote}

\begin{quote}
Type theory is nothin short of a grand unified theory of computation
unified with mathematics so ultimately there is no difference between
math and the code.

-- Robert Harper\cite{Harp13}
\end{quote}

\chapter{Here is a problem}

Proving programs correct involves working with a second programming
language, the proof language, that is well-founded on some theory.
Proofs (programs), can be reduced (compiled) in this new language
to the primitive constructs (machine language).

The ideal case would be that the programming language used, such as
Spad, can be iosmorphic, or better yet, syntactically the same as
the proof language. Unfortunately that is not (yet?) the case with Spad.

The COQ system language, Gallina, is the closest match to Spad.

\section{Setting up the problem}

The GCD function will be our first example of a proof.

The goal is to prove that Axiom's implementation of 
the Euclidean GCD algorithm is correct.

We need to be clear about what is to be proven. In this case, we
need to show that, given GCD(a,b),
\begin{enumerate}
\item {\bf GCD} is a function from $a\times b \Rightarrow c$
\item {\bf a} and {\bf b} are elements of the correct type
\item {\bf c}, the resullt, is the correct type
\item the meaning of {\bf divisor}
\item the meaning of a {\bf common divisor}
\item {\bf GCD} terminates
\end{enumerate}

We next need to set up the things we know in "the global environment",
generally referred to as {\bf E} in Coq.

Axiom's GCD is categorically defined to work over any Euclidean domain.
This means that the axioms of a Euclidean domain are globally available.
In fact, this is stronger than we need since
\begin{itemize}
\item commutative rings $\subset$ integral domains
\item integral domains $\subset$ integrally closed domains
\item integrally closed domains  $\subset$ GCD domains
\item GCD domains $\subset$ unique factorization domains
\item unique factorization domains $\subset$ principlal ideal domains
\item principlal ideal domains $\subset$ Euclidean domains
\end{itemize}

A Euclidean function on $R$ is a function $f$ from $\mathbb{R}\\\{0\}$
to the non-negative integers satisfying the following fundamental 
division-with-remainder property\cite{WikiED}:

$D(a,b)$ = set of common divisors of $a$ and $b$.

$\gcd(a,b) = \max{D(a,b)}$


\section{Axiom NNI GCD}

NonNegativeInteger inherits {\tt gcd} from {\tt Integer} up the ``add chain''
since it is a subtype of {\tt Integer}. {\tt Integer} has {\tt EuclideanDomain}
as an ancestor\cite{Book103}:
\begin{verbatim}
(1) -> getAncestors "Integer"

   (1)
   {AbelianGroup, AbelianMonoid, AbelianSemiGroup, Algebra, BasicType,
    BiModule, CancellationAbelianMonoid, CharacteristicZero, CoercibleTo,
    CombinatorialFunctionCategory, CommutativeRing, ConvertibleTo,
    DifferentialRing, EntireRing, EuclideanDomain, GcdDomain,
    IntegerNumberSystem, IntegralDomain, LeftModule, LeftOreRing,
    LinearlyExplicitRingOver, Module, Monoid, OpenMath, OrderedAbelianGroup,
    OrderedAbelianMonoid, OrderedAbelianSemiGroup,
    OrderedCancellationAbelianMonoid, OrderedIntegralDomain, OrderedRing,
    OrderedSet, PatternMatchable, PrincipalIdealDomain, RealConstant,
    RetractableTo, RightModule, Ring, Rng, SemiGroup, SetCategory, StepThrough,
    UniqueFactorizationDomain}
                                                            Type: Set(Symbol)
\end{verbatim}

From category {\tt EuclideanDomain} (EUCDOM) we find the implementation of 
the Euclidean GCD algorithm\cite{Book102}:
\begin{verbatim}
      gcd(x,y) ==                --Euclidean Algorithm
         x:=unitCanonical x
         y:=unitCanonical y
         while not zero? y repeat
            (x,y):= (y,x rem y)
            y:=unitCanonical y   -- this doesn't affect the
                                 -- correctness of Euclid's algorithm,
                                 -- but
                                 -- a) may improve performance
                                 -- b) ensures gcd(x,y)=gcd(y,x)
                                 --    if canonicalUnitNormal
         x
\end{verbatim}
The {\tt unitCanonical} function comes from the category 
{\tt IntegralDomain} (INTDOM) where we find:
\begin{verbatim}
    unitNormal: % -> Record(unit:%,canonical:%,associate:%)
        ++ unitNormal(x) tries to choose a canonical element
        ++ from the associate class of x.
        ++ The attribute canonicalUnitNormal, if asserted, means that
        ++ the "canonical" element is the same across all associates of x
        ++ if \spad{unitNormal(x) = [u,c,a]} then
        ++ \spad{u*c = x}, \spad{a*u = 1}.
    unitCanonical: % -> %
        ++ \spad{unitCanonical(x)} returns \spad{unitNormal(x).canonical}.
\end{verbatim}
implemented as
\begin{verbatim}
      UCA ==> Record(unit:%,canonical:%,associate:%)
      if not (% has Field) then
        unitNormal(x) == [1$%,x,1$%]$UCA -- the non-canonical definition
      unitCanonical(x) == unitNormal(x).canonical -- always true
      recip(x) == if zero? x then "failed" else _exquo(1$%,x)
      unit?(x) == (recip x case "failed" => false; true)
      if % has canonicalUnitNormal then
         associates?(x,y) ==
           (unitNormal x).canonical = (unitNormal y).canonical
       else
         associates?(x,y) ==
           zero? x => zero? y
           zero? y => false
           x exquo y case "failed" => false
           y exquo x case "failed" => false
           true
\end{verbatim}

Coq proves the following GCD function:
\begin{verbatim}
  Fixpoint gcd a b :=
    match a with
      | 0 => b
      | S a' => gcd (b mod (S a')) (S a')
    end.
\end{verbatim}

This can be translated directly to working Spad code:
\begin{verbatim}
  GCD(x:NNI,y:NNI):NNI ==
    zero? x => y
    GCD(y rem x,x)
\end{verbatim}
with the test case results of:
\begin{verbatim}
(1) -> GCD(2415,945)
   Compiling function mygcd2 with type (NonNegativeInteger,
      NonNegativeInteger) -> NonNegativeInteger 

   (1)  105
                                                        Type: PositiveInteger
(2) -> GCD(0,945)

   (2)  945
                                                        Type: PositiveInteger
(3) -> GCD(2415,0)

   (3)  2415
                                                        Type: PositiveInteger
(4) -> GCD(17,15)

   (4)  1
                                                        Type: PositiveInteger
\end{verbatim}

\section{Mathematics}

From Buchberger\cite{Buch97},

Define ``divides''
\[ t\vert a \Longleftrightarrow \exists u (t \cdot u = a)\]

Define ``greatest common divisor''
\[ {\rm GCD}(a,b) = \forall t\ max(t\vert a \land t\vert b)\]

Theorem:
\[ (t\vert a \land t\vert b) \Longleftrightarrow t\vert (a-b) \land t\vert b\]

Euclids' Algorithm
\[ a > b \Rightarrow {\rm GCD}(a,b) = {\rm GCD}(a-b,b)\]

By the definition of GCD we need to show that
\[\forall t\ max(t\vert a \land t\vert b) =
  \forall t\ max(t\vert (a-b) \land t\vert b)\]

Thus we need to show that
\[(t\vert a \land t\vert b) \Longleftrightarrow (t\vert (a-b) \land t\vert b)\]

Let $t$ be arbitrary but fixed and assume
\begin{equation}\label{eq1}(t\vert a \land t\vert b)\end{equation}

We have to show
\begin{equation}\label{eq2}t\vert (a-b)\end{equation}

and
\begin{equation}\label{eq3}t\vert b\end{equation}

Equation \ref{eq3} follows propositionally. For equation \ref{eq2},
by definition of ``divides'', we have to find a $w$ such that
\begin{equation}\label{eq4}t \cdot w = a-b\end{equation}

From \ref{eq1}, by definition of ``divides'', we know that for certain
$u$ and $v$
\[t \cdot u = a\]

and
\[t \cdot v - b\]

Hence,
\[ a-b = t \cdot u - t \cdot v\]

But
\[t \cdot u - t \cdot v = t \cdot (u - v)\]

So we need to find
\[w = u - v\]

and 
\[\textrm{Find w such that }t \cdot u - t \cdot v = t \cdot w\]


\section{Approaches}
There are several systems that could be applied to approach the proof.

The plan is to initially look at Coq and ACL2.  Coq seems to be
applicable at the Spad level. ACL2 seems to be applicable at the Lisp
level. Both levels are necessary for a proper proof.

Coq is very close to Spad in spirit so we can use it for the
high-level proofs.

ACL2 is a Lisp-level proof technology which can be used to prove
the Spad-to-Lisp level.

There is an LLVM to ACL2 translator which can be used to move from
the GCL Lisp level to the hardware since GCL compiles to C.
In particular, the "Vellvm: Verifying the LLVM" \cite{Zdan14}
project is important.

Quoting from Hardin \cite{Hard14}
\begin{quote}
LLVM is a register-based intermediate in Static Single Assignment
(SSA) form. As such, LLVM supports any number of registers, each of
which is only assigned once, statically (dynamically, of course, a
given register can be assigned any number of times). Appel has
observed that ``SSA form is a kind of functional programming''; this
observation, in turn, inspired us to build a translator from LLVM to
the applicative subset of Common Lisp accepted by the ACL2 theorem
prover. Our translator produces an executable ACL2 specification that
is able to efficiently support validation via testing, as the
generated ACL2 code features tail recursion, as well as in-place
updates via ACL2's single-threaded object (stobj) mechanism. In order
to ease the process of proving properties about these translated
functions, we have also developed a technique for reasoning about
tail-recursive ACL2 functions that execute in-place, utilizing a
formally proven ``bridge'' to primitive-recursive versions of those
functions operating on lists.
\end{quote}

{\center{\includegraphics{ps/v13llvmtoacl2.eps}}}

Hardin \cite{Hard13} describes the toolchain thus:
\begin{quote}
Our translation toolchain architecture is shown in Figure 1. The left
side of tthe figure depicts a typical compiler frontend producing LLVM
intermediate code. LLVM output can be produced either as a binary
``bitcode'' (.bc) file, or as text (.ll file). We chose to parse the
text form, producing an abstract syntax tree (AST) representation of
the LLVM program. Our translator then converts the AST to ACL2
source. The ACL2 source file can then be admitted into an ACL2
session, along with conjectures that one wishes to prove about the
code, which ACL2 processes mostly automatically. In addition to
proving theorems about the translated LLVM code, ACL2 can also be used
to execute test vectors at reasonable speed.
\end{quote}

Note that you can see the intermediate form from clang with
\begin{verbatim}
clang -O4 -S -emit-llvm foo.c
\end{verbatim}

Both Coq and the Hardin translator use OCAML \cite{OCAM14} so we will have to
learn that language.

\chapter{Theory}
The proof of the Euclidean algorithm has been known since Euclid.
We need to study an existing proof and use it to guide our use of
Coq along the same lines, if possible. Some of the ``obvious''
natural language statements may require Coq lemmas.

From WikiProof \cite{Wiki14a} we quote:

Let \[a, b \in \Z\] and $a \ne 0$ or $b \ne 0$.

The steps of the algorithm are:
\begin{enumerate}
\item Start with $(a,b)$ such that $\abs{a} \ge \abs{b}$.
If $b = 0$ then the task is complete and the GCD is $a$.
\item if $b \ne 0$ then you take the remainder $r$ of $a/b$.
\item set $a \leftarrow b$, $b \leftarrow r$ (and thus
$\abs{a} \ge \abs{b}$ again).
\item repeat these steps until $b = 0$
\end{enumerate}
Thus the GCD of $a$ and $b$ is the value of the variable $a$ 
at the end of the algorithm.

The proof is:

Suppose \[a, b \in \Z\] and $a or b \ne 0$.

From the {\bf division theorem}, $a = qb + r$ 
where $0 \le r \le \abs{b}$

From {\bf GCD with Remainder}, the GCD of $a$ and $b$ is also the GCD
of $b$ and $r$.

Therefore we may search instead for the $gcd(b,r)$.

Since $\abs{r} \ge \abs{b}$ and \[b \in \Z\],
we will reach $r = 0$ after finitely many steps.

At this point, $gcd(r,0) = r$ from {\bf GCD with Zero}.

We quote the {\bf Division Theorem} proof \cite{Wiki14b}:

For every pair of integers $a$, $b$ where $b \ne 0$, there exist unique
integers $q,r$ such that $a = qb + r$ and $0 \le r \le \abs{b}$.

\subsection{Hoare's axioms and gcd proof}
From Hoare\cite{Hoar69}
\begin{tabular}{rll}
A1 & $x+y = y+x$ & addition is commutative\\
A2 & $x\times y = y\times x$ & multiplication is commutative\\
A3 & $(x+y)+z = x+(y+z)$ &  addition is associative\\
A4 & $(x\times y)\times z = x\times (y\times z)$ &
multiplication is associative\\
A5 & $x\times (y+z) = x\times y + x\times z$ &
multiplication distributes through addition\\
A6 & $y \le x \to (x-y) + y = x$ & addition cancels subtraction\\
A7 & $x + 0 = x$ &\\
A8 & $x\times 0 = 0$ &\\
A9 & $x\times 1 = x$ &\\
\end{tabular}

{\bf D0 Axiom of Assignment}
\[\vdash P_0 \{x:=f\} P\]
where
\begin{itemize}
\item $x$ is a variable identifer
\item $f$ is an expression
\item $P_0$ is obtained from $P$ by substituting $f$ for all
occurences of $x$
\end{itemize}

\section{The Division Algorithm}
From Judson \cite{Juds15},

An Application of the Principle of Well-Ordering that we will use
often is the division algorithm.

{\bf Theorem 2.9 Division Algorithm} Let $a$ and $b$ be integers, with
$b > 0$. Then there exists unique integers $q$ and $r$ such that
\[a=bq+r\]
where $0 \le r < b$.

{\bf Proof}

Let $a$ and $b$ be integers. If $b = ak$ for some integer $k$, we write
$a \vert b$. An integer $d$ is called a {\sl common divisor} of $a$ and
\index{common divisor}
$b$ if $d \vert a$ and $d \vert b$. The {\sl greatest common divisor}
\index{greatest common divisor}
of integers $a$ and $b$ is a positive integer $d$ such that $d$ is
a common divisor of $a$ and $b$ and if $d^{'}$ is any other common
divisor of $a$ and $b$, then $d^{'} \vert d$. We write $d=gcd(a,b)$;
for example, $gcd(24,36)=12$ and $gcd(120,102)=6$. We say that two
integers $a$ and $b$ are {\sl relatively prime} if $gcd(a,b)=1$.
\index{relatively prime}

{\bf Theorem 2.10} Let $a$ and $b$ be nonzero integers. Then there
exist integers $r$ and $s$ such that
\[gcd(a,b)=ar+bs\]
Furthermore, the greatest common divisor of $a$ and $b$ is unique.

{\bf Proof}

{\bf Collary 2.11} Let $a$ and $b$ be two integers that are relatively
prime. Then there exist integers $r$ and $s$ such that
\[ar+bs=1\]

{\bf The Euclidean Algorithm}
\index{The Euclidean Algorithm}

Among other things, Theorem 2.10 allows us to compute the greatest 
common divisor of two integers.

{\bf Example 2.1.2} Let us compute the greatest common divisor of 945 and
2415. First observe that
\[
\begin{array}{rcl}
2415 & = & 945 \cdot 2 + 525\\
945 & = & 525 \cdot 1 + 420\\
525 & = & 420 \cdot 1 + 105\\
420 & = & 105 \cdot 4 + 0\\
\end{array}
\]

Reversing our steps, 105 divides 420, 105 divides 525, 105 divides 945,
and 105 divides 2415. Hence, 105 divides both 945 and 2415. If $d$ were
another common divisor of 945 and 2415, then $d$ would also have to
divide 105. Therefore, $gcd(945,2415)=105$.

If we work backward through the above sequence of equations, we can also
obtain numbers $r$ and $s$ such that
\[945r + 2415s = 105\]

\[
\begin{array}{rcl}
105 & = & 525 + (-1)\cdot 420\\
105 & = & 525 + (-1)\cdot [945+(-1)\cdot 525]\\
105 & = & 2\cdot 525+(-1)\cdot 945\\
105 & = & 2\cdot[2415+(-2)\cdot 945]+(-1)\cdot 945\\
105 & = & 2*2415+(-5)\cdot 945
\end{array}
\]
So $r=-5$ and $s-2$. Notice the $r$ and $s$ are not unique, since
$r=41$ and $s=-16$ would also work.

To compute $gcd(a,b)=d$, we are using repeated divisions to obtain
a decreasing sequence of positive integers
$r_1 > r_2 > \ldots > r_n = d$; that is

\[
\begin{array}{rcl}
b & = & aq_1+r_1\\
a & = & r_1q_2+r_2\\
r_1 & = & r_2q_3+r_3\\
\vdots\\
r_{n-2} & = & r_{n-1}q_n+r_n\\
r_{n-1} & = & r_nq_{n+1}
\end{array}
\]

To find $r$ and $s$ such that $ar+bs=d$, we begin with the last equation
and substitute results obtained from the previous equations:

\[
\begin{array}{rcl}
d & = & r_n\\
d & = & r_{n-2} - r_{n-1}q_n\\
d & = & r_{n-2} - q_n(r_{n-3}-q_{n-1}r_{n-2}\\
d & = & -q_nr_{n-3}+(1+q_nq_{n-1})r_{n-2}\\
\vdots\\
d & = & ra+sb
\end{array}
\]

\chapter{Software Details}
\section{Installed Software}
Install CLANG, LLVM
\begin{verbatim}
http://llvm.org/releases/download.html
\end{verbatim}
Install OCAML
\begin{verbatim}
sudo apt-get install ocaml
\end{verbatim}

An OCAML version of gcd would be written 
\begin{verbatim}
let rec gcd a b = if b = 0 then a else gcd b (a mod b)
val gcd : int -> int -> int = <fun>
\end{verbatim}

\chapter{Temporal Logic of Actions (TLA)}
\begin{quote}
{\bf Sloppiness is easier than precision and rigor}
-- Leslie Lamport\cite{Lamp14a}
\end{quote}

Leslie Lamport\cite{Lamp14}\cite{Lamp16} on $21^{st}$ Century Proofs.

A method of writing proofs is described that makes it harder to prove
things that are not true. The method, based on hierarchical
structuring, is simple and practical. The author's twenty years of
experience writing such proofs is discussed.

Lamport points out that proofs need rigor and precision.
Structure and Naming are important. Every step of the proof
names the facts it uses. 

Quoting from \cite{Lamp16}:

Broadly speaking, a TLA+ proof is a collection of {\sl claims},
arranged in a hierarchical structure which we describe below, where
each claim has an {\sl assertion} tht is either {\sl unjustified} or
justified by a collection of {\sl cited facts}. The purpose of TLAPS
is to check the user-provided proofs of theorems, that is, to check
that the hierarchy of claims indeed establishes the truth of the theorem
if the claims were true, and then to check that the assertion of every
justified claim indeed is implied {\sl by} its cited facts. If a TLA+
theorem has a proof with no unjustified claims, then, as a result of
checking the proof, TLAPS verifies the truth of the theorem.

\section{The algorithm}

The well-known Euclidean algorithm can be written in the PlusCal 
language as follows:
\begin{verbatim}
--algorithm Euclid {
  variables x \in 1..M, y \in 1..N, x0 = x, y0 = y;
  {
    while (x # y) {
      if (x < y) { y := y - x; }
      else { x := x-y; }
    };
    assert x = GCD(x0, y0) /\ y = GCD(x0, y0)
  }
\end{verbatim}

The PlusCal translator translates this algorithm into a TLA+ specification
that we could prove correct. However, in this tutorial, we shall write a
somewhat simpler specification of Euclid's algorithm directly in TLA+.

\subsection{Creating a new TLA+ module}

In order to get the definitions of arithmetic operators ($+$, $-$, etc.),
we shall make this specification {\sl extend} the {\tt Integers} 
standard module.

\begin{verbatim}
--------------------- Module Euclid ----------------------
EXTENDS Integers
\end{verbatim}

\subsection{Definitions}

We shall then define the GCD of two integers. For that purpose, let us
define the predicate ``p divides q'' as follows: p divides q iff there
exists some integer d in the interval 1..q such that q is equal to p
times d.

\begin{verbatim}
p | q == \E d \in 1..q : q = p * d
\end{verbatim}

We then define the set of divisors of an integer q as the sets of integers
which both belong to the interval 1..q and divide q:
\begin{verbatim}
Divisors(q) == {d \in 1..q : d | q}
\end{verbatim}

We define the maximum of a set S as one of the elements of this set which
is greater than or equal to all the other elements:
\begin{verbatim}
Maximum(S) == CHOOSE x \in S : \A y \in S : x >= y
\end{verbatim}

And finally, we define the GCD of two integers p and q to be the maximum
of the intersection of Divisors(p) and Divisors(a):
\begin{verbatim}
GCD(p,q) == Maximum(Divisors(p) \cap Divisors(q))
\end{verbatim}

For convenience, we shall also define the set of all positive integers as:
\begin{verbatim}
Number = Nat \ {0}
\end{verbatim}

\subsection{Constants and variables}
We then define the two constants and two variables needed to describe
the Euclidean algorithm, where M and N are the values whose GCD is to
be computed:
\begin{verbatim}
CONSTANTS M, N
VARIABLES x, y
\end{verbatim}

\subsection{The specification}
We define the initial state of the Euclidean algorithm as follows:
\begin{verbatim}
Init == (x = M) /\ (y = N)
\end{verbatim}

In the Euclidean algorithm, two actions can be performed:
\begin{itemize}
\item set the value of y to y - x if x $<$ y
\item set the value of x to x - y if x $>$ y
\end{itemize}

These actions are again written as a definition of {\tt Next}, which
specifies the next-state relation. In TLA+, a primed variable refers
to its value at the next state of the algorithm.
\begin{verbatim}
Next == \/ /\ x < y
           /\ y' = y - x
           /\ x' = x
        \/ /\ y < x
           /\ x' = x-y
           /\ y' = y
\end{verbatim}

The specification of the algorithm asserts that the variables have the
correct initial values and, in each execution step, either a {\tt Next}
action is performed or x and y keep the same values:
\begin{verbatim}
Spec == Init /\ [][Next]_<<x,y>>
\end{verbatim}

(For reasons that are irrelevant to this algorithm, TLA specifications
always allow {\sl stuttering steps} that leave all the variables
unchanged.)

We want to prove that the algorithm always satisfies the following
property:
\begin{verbatim}
ResultCorrect == (x = y) => x = GCD(M, N)
\end{verbatim}

Hence we want to prove the following theorem named Correctness:
\begin{verbatim}
THEOREM Correctness == Spec => []ResultCorrect
\end{verbatim}

\subsection{Summary}

\begin{verbatim}
--------------------- Module Euclid ----------------------
EXTENDS Integers

p | q == \E d \in 1..q : q = p * d
Divisors(q) == {d \in 1..q : d | q}
Maximum(S) == CHOOSE x \in S : \A y \in S : x >= y
GCD(p,q) == Maximum(Divisors(p) \cap Divisors(q))
Number == Nat \ {0}

CONSTANTS M, N
VARIABLES x, y

Init == (x = M) /\ (y = N)

Next == \/ /\ x < y
           /\ y' = y - x
           /\ x' = x
        \/ /\ y < x
           /\ x' = x-y
           /\ y' = y

Spec == Init /\ [][Next]_<<x,y>>

ResultCorrect == (x = y) => x = GCD(M,N)

THEOREM Correctness == Spec => []ResultCorrect
\end{verbatim}

\section{A simple proof}
\subsection{The invariant}

Intuitively, the theorem Correctness holds because the implementation
guarantees the following {\sl invariant}
\begin{verbatim}
InductiveInvariant == /\ x \in Number
                      /\ y \in Number
                      /\ GCD(x, y) = GCD(M, N)
\end{verbatim}

That is, {\tt InductiveInvariant} holds for the initial state (i.e.,
the state specified by {\tt Init}) and is preserved by the 
next-state relation {\tt [Next]\_}$<<x,y>>$

\subsection{Checking proofs}

First we need to assume that constants M and N are not equal to zero
\begin{verbatim}
ASSUME NumberAssumption == M \in Number /\ N \in Number
\end{verbatim}

Let us then prove that {\tt InductiveInvariant} holds for the initial state.
\begin{verbatim}
THEOREM InitProperty == Init => InductiveInvariant
\end{verbatim}

To check whether TLAPS can prove that theorem by itself, we declare
its proof obvious.
\begin{verbatim}
THEOREM InitProperty == Init => InductiveInvariant
  OBVIOUS
\end{verbatim}

We now ask TLAPS to prove that theorem. But TLAPS does not know how to
prove the proof obligation corresponding to that proof. It prints
that obligation and reports failures to three backends, Zenon, Isabelle,
and SMT. The default behavior of TLAPS is to send obligations first to
an SMT solver (by default CVC3), then if that fails to the automatic
prover Zenon, then if Zenon fails to Isabelle (with the tactic ``auto'').

\subsection{Using facts and definitions}
The obligation cannot be proved because TLAPS treats the symbols
{\tt Init} and {\tt InductiveInvariant} as opaque identifiers unless
it is explicitly instructed to expand their definitions using the
directive {\tt DEF}. The main purpose of this treantment of definitions
is to make proof-checking tractable, because expanding definitions can
arbitrarily increase the size of expressions. Explicit use of definitions
is also a good hint to the (human) reader to look only at the listed
definitions to understand a proof step. In tha precise case, we can ask
TLAPS to expand definitions of {\tt Init} and {\tt InductiveInvariant},
by replacing the proof {\tt OBVIOUS} by the proof\\
{\tt BY DEF Init, InductiveInvariant}. In the obligations sent to the
backends, the definitions of {\tt Init} and {\tt InductiveInvariant}
have been expanded. 

Unfortunately, none of the back-ends could prove that obligation. As with
{\tt definitions}, we have to specify which facts are {\sl usable}. In this
case, we have to make the fact {\tt NumberAssumption} usable by changing
the proof to
\begin{verbatim}
THEOREM InitProperty == Init => InductiveInvariant
  BY NumberAssumption DEF Init, InductiveInvariant
\end{verbatim}

The general form of a {\tt BY} proof is:
\[ {\tt BY\ } e_1,\ldots,e_m {\tt \ DEF\ } d_1,\ldots,d_n\]
which claims that the assertion follows by assuming $e_1,\ldots,e_m$
and expanding the definitions $d_1,\ldots,d_n$. It is the job of TLAPS
to then check this claim, and also to check that the cited facts
$e_1,\ldots,e_m$ are indeed true. 

Finally, SMT succeeds in proving that obligation.
\begin{verbatim}
--------------------- Module Euclid ----------------------
EXTENDS Integers

p | q == \E d \in 1..q : q = p * d
Divisors(q) == {d \in 1..q : d | q}
Maximum(S) == CHOOSE x \in S : \A y \in S : x >= y
GCD(p,q) == Maximum(Divisors(p) \cap Divisors(q))
Number == Nat \ {0}

CONSTANTS M, N
VARIABLES x, y

Init == (x = M) /\ (y = N)

Next == \/ /\ x < y
           /\ y' = y - x
           /\ x' = x
        \/ /\ y < x
           /\ x' = x-y
           /\ y' = y

Spec == Init /\ [][Next]_<<x,y>>

ResultCorrect == (x = y) => x = GCD(M,N)

InductiveInvariant == /\ x \in Number
                      /\ y \in Number
                      /\ GCD(x, y) = GCD(M, N)

ASSUME NumberAssumption == M \in Number /\ N \in Number

THEOREM InitProperty == Init => InductiveInvariant
  BY NumberAssumption DEF Init, InductiveInvariant

THEOREM Correctness == Spec => []ResultCorrect
\end{verbatim}


\section{Divisibility Definition}
In Shoup\cite{Sho08} we find the divisibility definition.

Given the integers, $a$ and $b$
\[ a {\rm \ divides\ } b \implies az = b {\rm \ for\ some\ } z\]
so or all $a$,$b$, and $c$
\[ a | a,\quad 1 | a,\quad and \quad a | 0\]
because $a\cdot 1 = a$, $1\cdot a = a$, and $a\cdot 0 = 0$
\[ 0 | a \iff a = 0\]
\[a | b \iff -a | b \iff a | -b\]
\[ a | b {\rm \ and\ } a | c \implies a | (b+c)\]
\[a | b {\rm \ and\ } b | c \implies a | c\]
\[a | b {\rm \ and\ } b \ne 0 \implies 1 \le |a| \le |b|\]
\[az = b \ne 0 {\rm \ and\ } a \ne 0 {\rm \ and\ }z \ne 0 
\implies |a| \ge 1 {\rm \ and\ } |z| \ge 1\]
\[ a | b {\rm \ and\ } b | a \implies a = \pm b\]
proof: \[ a | b \implies |a| \le |b|; b | a \implies |b| \le |a|; 
{\rm \ therefore\ } |a|=|b| \implies a = \pm b\]
\[ a | 1 \iff a = \pm 1\]

\chapter{COQ proof of GCD}
\section{Basics of the Calculus of Constructions}
Coquend\cite{Coqu86}\cite{Wiki17} 
defines the Calculus of Constructions which can
be considered an extension of the Curry-Howard Isomorphism. The components
are
\subsection{Terms}
A {\sl term} in the calculus of constructions is constructed using the
following rules:
\begin{itemize}
\item {\bf T} is a term (also called {\sl Type})
\item {\bf P} is a term (also called {\sl Prop}, the type of all propositions)
\item Variables ($x,y,\ldots$) are terms
\item if {\bf A} and {\bf B} are terms, then so are
\begin{itemize}
\item {\bf $(A,B)$}
\item {\bf $(\lambda x: A,B)$}
\item {\bf $(\forall x: A,B)$}
\end{itemize}
\end{itemize}

The calculus of constructions has five kinds of objects:
\begin{enumerate}
\item {\sl proofs}, which are terms whose types are {\sl propositions}
\item {\sl propositions}, which are also known as {\sl small types}
\item {\sl predicates}, which are functions that return propositions
\item {\sl large types}, which are the types of predicates. {\bf P} is
an example of a large type)
\item {\bf T} itself, which is the type of large types.
\end{enumerate}

\subsection{Judgements}

The calculus of constructions allows proving {\bf typing judgements}
\[x_1:A_1, x_2:A_2, \ldots \vdash t:B\]
which can be read as the implication
\[{\rm\ if\ variables\ }x_1, x_2, \ldots, 
{\rm\ have\ types\ }A_1, A_2,\ldots,
{\rm\ then\ term\ }t{\rm\ has\ type\ }B\]

The valid judgements for the calculus of constructions are derivable
from a set of inference rules. In the following, we use $\Gamma$
to mean a sequence of type assignments
$x_1:A_1$, $x_2:A_2,\ldots$,
and we use {\bf K} to mean either {\bf P} or {\bf T}.
We shall write {\bf $A:B:C$} to mean 
"{\bf $A$} has type {\bf $B$}, and {\bf $B$} has type {\bf $C$}".
We shall write {\bf $B(x:=N)$} to mean the result of substituting
the term {\bf $N$} for the variable {\bf $x$} in the term {\bf $B$}.

An inference rule is written in the form
\[\frac{{\bf \Gamma} \vdash {\bf A : B}}
{{\bf \Gamma^\prime} \vdash {\bf C : D}}\]
which means
\[{\rm\ if\ }{\bf \Gamma} \vdash {\bf A : B}
{\rm\ is\ a\ valid\ judgement,\ then\ so\ is\ }
{\bf \Gamma^\prime} \vdash {\bf C : D}\]

\subsection{Inference Rules}

In Frade\cite{Frad08} we find:

\[\begin{array}{lcl}

({\rm axiom}) &
() \vdash s_1 : s_2 
& {\rm\ if\ }(s_1,s_2)\in A\\
&&\\

({\rm start}) &
\displaystyle{\frac{\Gamma \vdash A : s}{\Gamma,x : A \vdash x : A}}
& {\rm\ if\ }x\notin {\rm\ dom}(\Gamma) \\
&&\\

({\rm weakening}) &
\displaystyle{\frac{\Gamma \vdash M : A\quad\quad\Gamma\vdash B : s}
{\Gamma,x : B \vdash M : A}}
& {\rm\ if\ }x\notin {\rm\ dom}(\Gamma) \\
&&\\

({\rm product}) &
\displaystyle{\frac{\Gamma\vdash A : s_1\quad\quad\Gamma,x : A\vdash B : s_2}
{\Gamma\vdash(\prod x:A. B) : s_3}}
& {\rm\ if\ }(s_1,s_2,s_3) \in \mathbb{R}\\
&&\\

({\rm application}) &
\displaystyle{\frac{\Gamma\vdash M:(\prod x : A. B)\quad\quad\Gamma \vdash N:A}
{\Gamma \vdash MN : B[x:=N]}} & \\
&&\\

({\rm abstraction}) & 
\displaystyle{\frac{\Gamma,x:A\vdash M:B\quad\quad\Gamma\vdash(\prod x:A. B):s}
{\Gamma\vdash\lambda x : A.M : (\prod x:A.B)}} & \\
&&\\

({\rm conversion}) &
\displaystyle{\frac{\Gamma\vdash M : A\quad\quad\Gamma\vdash B : s}
{\Gamma\vdash M : B}}
& {\rm\ if\ }A =_\beta B\\

\end{array}\]

\subsection{Defining Logical Operators}

\[\begin{array}{cccc}
A\Rightarrow B & \equiv & \forall x:A.B & (x\notin B)\\
A\wedge B & \equiv & \forall C:P.(A\Rightarrow B\Rightarrow C)\Rightarrow C &\\
A\vee B & \equiv & \forall C:P.(A\Rightarrow C)\Rightarrow(B\Rightarrow C)
\Rightarrow C &\\
\neg A & \equiv & \forall C : P.(A\Rightarrow C) &\\
\exists x : A. B & \equiv & \forall C : P.(\forall x : A.(B\Rightarrow C))
\Rightarrow C &\\
\end{array}\]

\subsection{Defining Types}

The basic data types used in computer science can be defined within the
Calculus of Constructions:

{\bf Booleans}

\[\forall A : P.A \Rightarrow A \Rightarrow A\]

{\bf Naturals}

\[\forall A : P.(A \Rightarrow A) \Rightarrow (A \Rightarrow A)\]

{\bf Product} $A\times B$

\[A \wedge B\]

{\bf Disjoint Union} $A + B$

\[A \vee B\]

Note that Booleans and Naturals are defined in the same way as in
Church encoding. However additional problems raise from propositional
extensionality and proof irrelevance.

\section{Why does COQ have Prop?}

From a stackexchange post\cite{Stac17} we find the question:

"Coq has a type {\tt Prop} of proof irrelevant propositions which are discarded
during extraction. What are the reasons for having this if we use Coq only
for proofs? Prop is impredicative, however, Coq automatically infers
universe indexes and we can use {\tt Type(i)} instead everywhere. 
It seems {\tt Prop} complicates everything a lot."

Prop is very useful for program extraction because it allows us to 
delete parts of code that are useless. For example, to extract a sorting
algorithm we would prove the statement ``for every list $l$ there
is a list $k$ such that $k$ is ordered and $k$ is a permutation of $l$''.
If we write this down in Coq and extract without using Prop, we will get:
\begin{enumerate}
\item ``for all $l$ there is a $k$'' which gives us a map {\bf sort} which
takes lists to lists,
\item ``such that $k$ is ordered'' will give a function {\bf verify} which
runs through $k$ and checks that it is sorted, and
\item ``$k$ is a permutation of $l$ will give a permutation {\bf p1} which
takes $l$ to $k$. Note that {\bf p1} is not just a mapping, but also the
inverse mapping together with programs verifying that the two maps really
are inverses.
\end{enumerate}

While the extra stuff is not totally useless, in many applications we want
to get ride of it and keep just {\bf sort}. This can be accomplished if we
use {\tt Prop} to state ``$k$ is ordered'' and ``$k$ is a permutation of $l$'',
but {\sl not} ``for all $l$ there is $k$''.

In general, a common way to extract code is to consider a statement of the
form
\[\forall x : A. \exists y : B. \phi(x,y)\]
where $x$ is input, $y$ is output, and $\phi(x,y)$ explains what it means
for $y$ to be a correct output. (In the above example $A$ and $B$ are the
types of lists and $\phi(l,k)$ is ''$k$ is ordered and $k$ is a permutation
of $l$.'') if $\phi$ is in {\tt Prop} then extraction gives a map 
$f : A\Rightarrow B$ such that $\phi(x,f(x))$ holds for all $x\in A$. 
If $\phi$ is in {\tt Set} then we also get a function $g$ such that $g(x)$
is the proof that $\phi(x,f(x))$ holds, for all $x\in A$. Often the proof
is computationally useless and we prefer to get rid of it, especially when
it is nested deeply inside some other statement. {\tt Prop} gives use the
possibility to do so.

There is a question whether we could avoid {\tt Prop} altogether by
automatically optimizing away ``useless extracted code''. To some extent
we can do that, for instance all code extracted from the negative fragment
of logic (stuff build from the empty type, unit type, products) is useless
as it just shuffles around the unit. But there are genuine design decisions
one has to make when using {\tt Prop}. Here is a simple example, where
$\sum$ means that we are in {\tt Type} and $\exists$ means we are in 
{\tt Prop}. If we extract from
\[\prod_{n:N}\sum_{b:[0,1]}\sum_{k:N}\quad n=2\cdot k+b\]
we will gat a program which decomposes $n$ into its lowest bit $b$ and 
the remaining bits $k$, i.e., it computes everything. If we extract from
\[\prod_{n:N}
\sum_{b:[0,1]}
\underset{k:N}{\scalebox{2}{\raisebox{-2pt}{\ensuremath{\exists}}}}
\quad n=2\cdot k+b\]

then the program will only compute the lowest bit $b$. The machine cannot
tell which is the correct one, the user has to tell it what he wants.

\section{Source code of COQ GCD Proof}
This is the proof of GCD\cite{Coqu16a} in the COQ\cite{Coqu16} sources:
\begin{verbatim}
Library Coq.ZArith.Znumtheory

Require Import ZArith_base.
Require Import ZArithRing.
Require Import Zcomplements.
Require Import Zdiv.
Require Import Wf_nat.

For compatibility reasons, this Open Scope isn't local as it should

Open Scope Z_scope.

This file contains some notions of number theory upon Z numbers:

    a divisibility predicate Z.divide
    a gcd predicate gcd
    Euclid algorithm euclid
    a relatively prime predicate rel_prime
    a prime predicate prime
    properties of the efficient Z.gcd function


Notation Zgcd := Z.gcd (compat "8.3").
Notation Zggcd := Z.ggcd (compat "8.3").
Notation Zggcd_gcd := Z.ggcd_gcd (compat "8.3").
Notation Zggcd_correct_divisors := Z.ggcd_correct_divisors (compat "8.3").
Notation Zgcd_divide_l := Z.gcd_divide_l (compat "8.3").
Notation Zgcd_divide_r := Z.gcd_divide_r (compat "8.3").
Notation Zgcd_greatest := Z.gcd_greatest (compat "8.3").
Notation Zgcd_nonneg := Z.gcd_nonneg (compat "8.3").
Notation Zggcd_opp := Z.ggcd_opp (compat "8.3").

The former specialized inductive predicate Z.divide is now a generic existential predicate.

Notation Zdivide := Z.divide (compat "8.3").

Its former constructor is now a pseudo-constructor.

Definition Zdivide_intro a b q (H:b=q*a) : Z.divide a b := ex_intro _ q H.

Results concerning divisibility

Notation Zdivide_refl := Z.divide_refl (compat "8.3").
Notation Zone_divide := Z.divide_1_l (compat "8.3").
Notation Zdivide_0 := Z.divide_0_r (compat "8.3").
Notation Zmult_divide_compat_l := Z.mul_divide_mono_l (compat "8.3").
Notation Zmult_divide_compat_r := Z.mul_divide_mono_r (compat "8.3").
Notation Zdivide_plus_r := Z.divide_add_r (compat "8.3").
Notation Zdivide_minus_l := Z.divide_sub_r (compat "8.3").
Notation Zdivide_mult_l := Z.divide_mul_l (compat "8.3").
Notation Zdivide_mult_r := Z.divide_mul_r (compat "8.3").
Notation Zdivide_factor_r := Z.divide_factor_l (compat "8.3").
Notation Zdivide_factor_l := Z.divide_factor_r (compat "8.3").

Lemma Zdivide_opp_r a b : (a | b) -> (a | - b).

Lemma Zdivide_opp_r_rev a b : (a | - b) -> (a | b).

Lemma Zdivide_opp_l a b : (a | b) -> (- a | b).

Lemma Zdivide_opp_l_rev a b : (- a | b) -> (a | b).

Theorem Zdivide_Zabs_l a b : (Z.abs a | b) -> (a | b).

Theorem Zdivide_Zabs_inv_l a b : (a | b) -> (Z.abs a | b).

Hint Resolve Z.divide_refl Z.divide_1_l Z.divide_0_r: zarith.
Hint Resolve Z.mul_divide_mono_l Z.mul_divide_mono_r: zarith.
Hint Resolve Z.divide_add_r Zdivide_opp_r Zdivide_opp_r_rev Zdivide_opp_l
  Zdivide_opp_l_rev Z.divide_sub_r Z.divide_mul_l Z.divide_mul_r
  Z.divide_factor_l Z.divide_factor_r: zarith.

Auxiliary result.

Lemma Zmult_one x y : x >= 0 -> x * y = 1 -> x = 1.

Only 1 and -1 divide 1.

Notation Zdivide_1 := Z.divide_1_r (compat "8.3").

If a divides b and b divides a then a is b or -b.

Notation Zdivide_antisym := Z.divide_antisym (compat "8.3").
Notation Zdivide_trans := Z.divide_trans (compat "8.3").

If a divides b and b<>0 then |a| <= |b|.

Lemma Zdivide_bounds a b : (a | b) -> b <> 0 -> Z.abs a <= Z.abs b.

Z.divide can be expressed using Z.modulo.

Lemma Zmod_divide : forall a b, b<>0 -> a mod b = 0 -> (b | a).

Lemma Zdivide_mod : forall a b, (b | a) -> a mod b = 0.

Z.divide is hence decidable

Lemma Zdivide_dec a b : {(a | b)} + {~ (a | b)}.

Theorem Zdivide_Zdiv_eq a b : 0 < a -> (a | b) -> b = a * (b / a).

Theorem Zdivide_Zdiv_eq_2 a b c :
 0 < a -> (a | b) -> (c * b) / a = c * (b / a).

Theorem Zdivide_le: forall a b : Z,
 0 <= a -> 0 < b -> (a | b) -> a <= b.

Theorem Zdivide_Zdiv_lt_pos a b :
 1 < a -> 0 < b -> (a | b) -> 0 < b / a < b .

Lemma Zmod_div_mod n m a:
 0 < n -> 0 < m -> (n | m) -> a mod n = (a mod m) mod n.

Lemma Zmod_divide_minus a b c:
 0 < b -> a mod b = c -> (b | a - c).

Lemma Zdivide_mod_minus a b c:
 0 <= c < b -> (b | a - c) -> a mod b = c.

Greatest common divisor (gcd).
There is no unicity of the gcd; hence we define the predicate Zis_gcd a b g expressing that g is a gcd of a and b. (We show later that the gcd is actually unique if we discard its sign.)

Inductive Zis_gcd (a b g:Z) : Prop :=
 Zis_gcd_intro :
  (g | a) ->
  (g | b) ->
  (forall x, (x | a) -> (x | b) -> (x | g)) ->
  Zis_gcd a b g.

Trivial properties of gcd

Lemma Zis_gcd_sym : forall a b d, Zis_gcd a b d -> Zis_gcd b a d.

Lemma Zis_gcd_0 : forall a, Zis_gcd a 0 a.

Lemma Zis_gcd_1 : forall a, Zis_gcd a 1 1.

Lemma Zis_gcd_refl : forall a, Zis_gcd a a a.

Lemma Zis_gcd_minus : forall a b d, Zis_gcd a (- b) d -> Zis_gcd b a d.

Lemma Zis_gcd_opp : forall a b d, Zis_gcd a b d -> Zis_gcd b a (- d).

Lemma Zis_gcd_0_abs a : Zis_gcd 0 a (Z.abs a).

Hint Resolve Zis_gcd_sym Zis_gcd_0 Zis_gcd_minus Zis_gcd_opp: zarith.

Theorem Zis_gcd_unique: forall a b c d : Z,
 Zis_gcd a b c -> Zis_gcd a b d -> c = d \/ c = (- d).

Extended Euclid algorithm.
Euclid's algorithm to compute the gcd mainly relies on the following property.

Lemma Zis_gcd_for_euclid :
  forall a b d q:Z, Zis_gcd b (a - q * b) d -> Zis_gcd a b d.

Lemma Zis_gcd_for_euclid2 :
  forall b d q r:Z, Zis_gcd r b d -> Zis_gcd b (b * q + r) d.

We implement the extended version of Euclid's algorithm, i.e. the one computing Bezout's coefficients as it computes the gcd. We follow the algorithm given in Knuth's "Art of Computer Programming", vol 2, page 325.

Section extended_euclid_algorithm.

  Variables a b : Z.

The specification of Euclid's algorithm is the existence of u, v and d such that ua+vb=d and (gcd a b d).

  Inductive Euclid : Set :=
    Euclid_intro :
    forall u v d:Z, u * a + v * b = d -> Zis_gcd a b d -> Euclid.

The recursive part of Euclid's algorithm uses well-founded recursion of non-negative integers. It maintains 6 integers u1,u2,u3,v1,v2,v3 such that the following invariant holds: u1*a+u2*b=u3 and v1*a+v2*b=v3 and gcd(u3,v3)=gcd(a,b).

  Lemma euclid_rec :
    forall v3:Z,
      0 <= v3 ->
      forall u1 u2 u3 v1 v2:Z,
        u1 * a + u2 * b = u3 ->
        v1 * a + v2 * b = v3 ->
        (forall d:Z, Zis_gcd u3 v3 d -> Zis_gcd a b d) -> Euclid.

We get Euclid's algorithm by applying euclid_rec on 1,0,a,0,1,b when b>=0 and 1,0,a,0,-1,-b when b<0.

  Lemma euclid : Euclid.

End extended_euclid_algorithm.

Theorem Zis_gcd_uniqueness_apart_sign :
  forall a b d d':Z, Zis_gcd a b d -> Zis_gcd a b d' -> d = d' \/ d = - d'.

Bezout's coefficients

Inductive Bezout (a b d:Z) : Prop :=
  Bezout_intro : forall u v:Z, u * a + v * b = d -> Bezout a b d.

Existence of Bezout's coefficients for the gcd of a and b

Lemma Zis_gcd_bezout : forall a b d:Z, Zis_gcd a b d -> Bezout a b d.

gcd of ca and cb is c gcd(a,b).

Lemma Zis_gcd_mult :
  forall a b c d:Z, Zis_gcd a b d -> Zis_gcd (c * a) (c * b) (c * d).

Relative primality

Definition rel_prime (a b:Z) : Prop := Zis_gcd a b 1.

Bezout's theorem: a and b are relatively prime if and only if there exist u and v such that ua+vb = 1.

Lemma rel_prime_bezout : forall a b:Z, rel_prime a b -> Bezout a b 1.

Lemma bezout_rel_prime : forall a b:Z, Bezout a b 1 -> rel_prime a b.

Gauss's theorem: if a divides bc and if a and b are relatively prime, then a divides c.

Theorem Gauss : forall a b c:Z, (a | b * c) -> rel_prime a b -> (a | c).

If a is relatively prime to b and c, then it is to bc

Lemma rel_prime_mult :
  forall a b c:Z, rel_prime a b -> rel_prime a c -> rel_prime a (b * c).

Lemma rel_prime_cross_prod :
  forall a b c d:Z,
    rel_prime a b ->
    rel_prime c d -> b > 0 -> d > 0 -> a * d = b * c -> a = c /\ b = d.

After factorization by a gcd, the original numbers are relatively prime.

Lemma Zis_gcd_rel_prime :
  forall a b g:Z,
    b > 0 -> g >= 0 -> Zis_gcd a b g -> rel_prime (a / g) (b / g).

Theorem rel_prime_sym: forall a b, rel_prime a b -> rel_prime b a.

Theorem rel_prime_div: forall p q r,
 rel_prime p q -> (r | p) -> rel_prime r q.

Theorem rel_prime_1: forall n, rel_prime 1 n.

Theorem not_rel_prime_0: forall n, 1 < n -> ~ rel_prime 0 n.

Theorem rel_prime_mod: forall p q, 0 < q ->
 rel_prime p q -> rel_prime (p mod q) q.

Theorem rel_prime_mod_rev: forall p q, 0 < q ->
 rel_prime (p mod q) q -> rel_prime p q.

Theorem Zrel_prime_neq_mod_0: forall a b, 1 < b -> rel_prime a b -> a mod b <> 0.

Primality

Inductive prime (p:Z) : Prop :=
  prime_intro :
    1 < p -> (forall n:Z, 1 <= n < p -> rel_prime n p) -> prime p.

The sole divisors of a prime number p are -1, 1, p and -p.

Lemma prime_divisors :
  forall p:Z,
    prime p -> forall a:Z, (a | p) -> a = -1 \/ a = 1 \/ a = p \/ a = - p.

A prime number is relatively prime with any number it does not divide

Lemma prime_rel_prime :
  forall p:Z, prime p -> forall a:Z, ~ (p | a) -> rel_prime p a.

Hint Resolve prime_rel_prime: zarith.

As a consequence, a prime number is relatively prime with smaller numbers

Theorem rel_prime_le_prime:
 forall a p, prime p -> 1 <= a < p -> rel_prime a p.

If a prime p divides ab then it divides either a or b

Lemma prime_mult :
  forall p:Z, prime p -> forall a b:Z, (p | a * b) -> (p | a) \/ (p | b).

Lemma not_prime_0: ~ prime 0.

Lemma not_prime_1: ~ prime 1.

Lemma prime_2: prime 2.

Theorem prime_3: prime 3.

Theorem prime_ge_2 p : prime p -> 2 <= p.

Definition prime' p := 1<p /\ (forall n, 1<n<p -> ~ (n|p)).

Lemma Z_0_1_more x : 0<=x -> x=0 \/ x=1 \/ 1<x.

Theorem prime_alt p : prime' p <-> prime p.

Theorem square_not_prime: forall a, ~ prime (a * a).

Theorem prime_div_prime: forall p q,
 prime p -> prime q -> (p | q) -> p = q.

we now prove that Z.gcd is indeed a gcd in the sense of Zis_gcd.

Notation Zgcd_is_pos := Z.gcd_nonneg (compat "8.3").

Lemma Zgcd_is_gcd : forall a b, Zis_gcd a b (Z.gcd a b).

Theorem Zgcd_spec : forall x y : Z, {z : Z | Zis_gcd x y z /\ 0 <= z}.

Theorem Zdivide_Zgcd: forall p q r : Z,
 (p | q) -> (p | r) -> (p | Z.gcd q r).

Theorem Zis_gcd_gcd: forall a b c : Z,
 0 <= c -> Zis_gcd a b c -> Z.gcd a b = c.

Notation Zgcd_inv_0_l := Z.gcd_eq_0_l (compat "8.3").
Notation Zgcd_inv_0_r := Z.gcd_eq_0_r (compat "8.3").

Theorem Zgcd_div_swap0 : forall a b : Z,
 0 < Z.gcd a b ->
 0 < b ->
 (a / Z.gcd a b) * b = a * (b/Z.gcd a b).

Theorem Zgcd_div_swap : forall a b c : Z,
 0 < Z.gcd a b ->
 0 < b ->
 (c * a) / Z.gcd a b * b = c * a * (b/Z.gcd a b).

Notation Zgcd_comm := Z.gcd_comm (compat "8.3").

Lemma Zgcd_ass a b c : Z.gcd (Z.gcd a b) c = Z.gcd a (Z.gcd b c).

Notation Zgcd_Zabs := Z.gcd_abs_l (compat "8.3").
Notation Zgcd_0 := Z.gcd_0_r (compat "8.3").
Notation Zgcd_1 := Z.gcd_1_r (compat "8.3").

Hint Resolve Z.gcd_0_r Z.gcd_1_r : zarith.

Theorem Zgcd_1_rel_prime : forall a b,
 Z.gcd a b = 1 <-> rel_prime a b.

Definition rel_prime_dec: forall a b,
 { rel_prime a b }+{ ~ rel_prime a b }.

Definition prime_dec_aux:
 forall p m,
  { forall n, 1 < n < m -> rel_prime n p } +
  { exists n, 1 < n < m /\ ~ rel_prime n p }.

Definition prime_dec: forall p, { prime p }+{ ~ prime p }.

Theorem not_prime_divide:
 forall p, 1 < p -> ~ prime p -> exists n, 1 < n < p /\ (n | p).
 
\end{verbatim}

\chapter{LEAN proof of GCD}
This is the proof of GCD\cite{Avig14} in the LEAN\cite{Avig16} sources:
\begin{verbatim}
/-
Copyright (c) 2014 Jeremy Avigad. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Authors: Jeremy Avigad, Leonardo de Moura

Definitions and properties of gcd, lcm, and coprime.
-/
import .div
open eq.ops well_founded decidable prod

namespace nat

/- gcd -/

private definition pair_nat.lt : nat × nat → nat × nat → Prop := measure pr₂
private definition pair_nat.lt.wf : well_founded pair_nat.lt :=
intro_k (measure.wf pr₂) 20  -- we use intro_k to be able to execute gcd efficiently in the kernel

local attribute pair_nat.lt.wf [instance]      -- instance will not be saved in .olean
local infixl ` ≺ `:50 := pair_nat.lt

private definition gcd.lt.dec (x y₁ : nat) : (succ y₁, x % succ y₁) ≺ (x, succ y₁) :=
!mod_lt (succ_pos y₁)

definition gcd.F : Π (p₁ : nat × nat), (Π p₂ : nat × nat, p₂ ≺ p₁ → nat) → nat
| (x, 0)      f := x
| (x, succ y) f := f (succ y, x % succ y) !gcd.lt.dec

definition gcd (x y : nat) := fix gcd.F (x, y)

theorem gcd_zero_right [simp] (x : nat) : gcd x 0 = x := rfl

theorem gcd_succ [simp] (x y : nat) : gcd x (succ y) = gcd (succ y) (x % succ y) :=
well_founded.fix_eq gcd.F (x, succ y)

theorem gcd_one_right (n : ℕ) : gcd n 1 = 1 :=
calc gcd n 1 = gcd 1 (n % 1)  : gcd_succ
         ... = gcd 1 0        : mod_one

theorem gcd_def (x : ℕ) : Π (y : ℕ), gcd x y = if y = 0 then x else gcd y (x % y)
| 0        := !gcd_zero_right
| (succ y) := !gcd_succ ⬝ (if_neg !succ_ne_zero)⁻¹


theorem gcd_self : Π (n : ℕ), gcd n n = n
| 0         := rfl
| (succ n₁) := calc
    gcd (succ n₁) (succ n₁) = gcd (succ n₁) (succ n₁ % succ n₁) : gcd_succ
                      ...   = gcd (succ n₁) 0                     : mod_self

theorem gcd_zero_left : Π (n : ℕ), gcd 0 n = n
| 0         := rfl
| (succ n₁) := calc
    gcd 0 (succ n₁) = gcd (succ n₁) (0 % succ n₁) : gcd_succ
                ... = gcd (succ n₁) 0               : zero_mod

theorem gcd_of_pos (m : ℕ) {n : ℕ} (H : n > 0) : gcd m n = gcd n (m % n) :=
gcd_def m n ⬝ if_neg (ne_zero_of_pos H)

theorem gcd_rec (m n : ℕ) : gcd m n = gcd n (m % n) :=
by_cases_zero_pos n
  (calc
          m = gcd 0 m       : gcd_zero_left
        ... = gcd 0 (m % 0) : mod_zero)
  (take n, assume H : 0 < n, gcd_of_pos m H)

theorem gcd.induction {P : ℕ → ℕ → Prop}
                   (m n : ℕ)
                   (H0 : ∀m, P m 0)
                   (H1 : ∀m n, 0 < n → P n (m % n) → P m n) :
                 P m n :=
induction (m, n) (prod.rec (λm, nat.rec (λ IH, H0 m)
   (λ n₁ v (IH : ∀p₂, p₂ ≺ (m, succ n₁) → P (pr₁ p₂) (pr₂ p₂)),
      H1 m (succ n₁) !succ_pos (IH _ !gcd.lt.dec))))

theorem gcd_dvd (m n : ℕ) : (gcd m n ∣ m) ∧ (gcd m n ∣ n) :=
gcd.induction m n
  (take m, and.intro (!one_mul ▸ !dvd_mul_left) !dvd_zero)
  (take m n (npos : 0 < n), and.rec
     (assume (IH₁ : gcd n (m % n) ∣ n) (IH₂ : gcd n (m % n) ∣ (m % n)),
    have H : (gcd n (m % n) ∣ (m / n * n + m % n)), from
      dvd_add (dvd.trans IH₁ !dvd_mul_left) IH₂,
    have H1 : (gcd n (m % n) ∣ m), from !eq_div_mul_add_mod⁻¹ ▸ H,
    show (gcd m n ∣ m) ∧ (gcd m n ∣ n), from !gcd_rec⁻¹ ▸ (and.intro H1 IH₁)))

theorem gcd_dvd_left (m n : ℕ) : gcd m n ∣ m := and.left !gcd_dvd

theorem gcd_dvd_right (m n : ℕ) : gcd m n ∣ n := and.right !gcd_dvd

theorem dvd_gcd {m n k : ℕ} : k ∣ m → k ∣ n → k ∣ gcd m n :=
gcd.induction m n (take m, imp.intro)
  (take m n (npos : n > 0)
    (IH : k ∣ n → k ∣ m % n → k ∣ gcd n (m % n))
    (H1 : k ∣ m) (H2 : k ∣ n),
    have H3 : k ∣ m / n * n + m % n, from !eq_div_mul_add_mod ▸ H1,
    have H4 : k ∣ m % n, from nat.dvd_of_dvd_add_left H3 (dvd.trans H2 !dvd_mul_left),
    !gcd_rec⁻¹ ▸ IH H2 H4)

theorem gcd.comm (m n : ℕ) : gcd m n = gcd n m :=
dvd.antisymm
  (dvd_gcd !gcd_dvd_right !gcd_dvd_left)
  (dvd_gcd !gcd_dvd_right !gcd_dvd_left)

theorem gcd.assoc (m n k : ℕ) : gcd (gcd m n) k = gcd m (gcd n k) :=
dvd.antisymm
  (dvd_gcd
    (dvd.trans !gcd_dvd_left !gcd_dvd_left)
    (dvd_gcd (dvd.trans !gcd_dvd_left !gcd_dvd_right) !gcd_dvd_right))
  (dvd_gcd
    (dvd_gcd !gcd_dvd_left (dvd.trans !gcd_dvd_right !gcd_dvd_left))
    (dvd.trans !gcd_dvd_right !gcd_dvd_right))

theorem gcd_one_left (m : ℕ) : gcd 1 m = 1 :=
!gcd.comm ⬝ !gcd_one_right

theorem gcd_mul_left (m n k : ℕ) : gcd (m * n) (m * k) = m * gcd n k :=
gcd.induction n k
  (take n, calc gcd (m * n) (m * 0) = gcd (m * n) 0 : mul_zero)
  (take n k,
    assume H : 0 < k,
    assume IH : gcd (m * k) (m * (n % k)) = m * gcd k (n % k),
    calc
      gcd (m * n) (m * k) = gcd (m * k) (m * n % (m * k)) : !gcd_rec
                      ... = gcd (m * k) (m * (n % k))     : mul_mod_mul_left
                      ... = m * gcd k (n % k)             : IH
                      ... = m * gcd n k                   : !gcd_rec)

theorem gcd_mul_right (m n k : ℕ) : gcd (m * n) (k * n) = gcd m k * n :=
calc
  gcd (m * n) (k * n) = gcd (n * m) (k * n) : mul.comm
                  ... = gcd (n * m) (n * k) : mul.comm
                  ... = n * gcd m k         : gcd_mul_left
                  ... = gcd m k * n         : mul.comm

theorem gcd_pos_of_pos_left {m : ℕ} (n : ℕ) (mpos : m > 0) : gcd m n > 0 :=
pos_of_dvd_of_pos !gcd_dvd_left mpos

theorem gcd_pos_of_pos_right (m : ℕ) {n : ℕ} (npos : n > 0) : gcd m n > 0 :=
pos_of_dvd_of_pos !gcd_dvd_right npos

theorem eq_zero_of_gcd_eq_zero_left {m n : ℕ} (H : gcd m n = 0) : m = 0 :=
or.elim (eq_zero_or_pos m)
  (assume H1, H1)
  (assume H1 : m > 0, absurd H⁻¹ (ne_of_lt (!gcd_pos_of_pos_left H1)))

theorem eq_zero_of_gcd_eq_zero_right {m n : ℕ} (H : gcd m n = 0) : n = 0 :=
eq_zero_of_gcd_eq_zero_left (!gcd.comm ▸ H)

theorem gcd_div {m n k : ℕ} (H1 : k ∣ m) (H2 : k ∣ n) :
  gcd (m / k) (n / k) = gcd m n / k :=
or.elim (eq_zero_or_pos k)
  (assume H3 : k = 0, by subst k; rewrite *nat.div_zero)
  (assume H3 : k > 0, (nat.div_eq_of_eq_mul_left H3 (calc
        gcd m n = gcd m (n / k * k)             : nat.div_mul_cancel H2
            ... = gcd (m / k * k) (n / k * k) : nat.div_mul_cancel H1
            ... = gcd (m / k) (n / k) * k     : gcd_mul_right))⁻¹)

theorem gcd_dvd_gcd_mul_left (m n k : ℕ) : gcd m n ∣ gcd (k * m) n :=
dvd_gcd (dvd.trans !gcd_dvd_left !dvd_mul_left) !gcd_dvd_right

theorem gcd_dvd_gcd_mul_right (m n k : ℕ) : gcd m n ∣ gcd (m * k) n :=
!mul.comm ▸ !gcd_dvd_gcd_mul_left

theorem gcd_dvd_gcd_mul_left_right (m n k : ℕ) : gcd m n ∣ gcd m (k * n) :=
dvd_gcd  !gcd_dvd_left (dvd.trans !gcd_dvd_right !dvd_mul_left)

theorem gcd_dvd_gcd_mul_right_right (m n k : ℕ) : gcd m n ∣ gcd m (n * k) :=
!mul.comm ▸ !gcd_dvd_gcd_mul_left_right

/- lcm -/

definition lcm (m n : ℕ) : ℕ := m * n / (gcd m n)

theorem lcm.comm (m n : ℕ) : lcm m n = lcm n m :=
calc
  lcm m n = m * n / gcd m n : rfl
      ... = n * m / gcd m n : mul.comm
      ... = n * m / gcd n m : gcd.comm
      ... = lcm n m           : rfl

theorem lcm_zero_left (m : ℕ) : lcm 0 m = 0 :=
calc
  lcm 0 m = 0 * m / gcd 0 m : rfl
      ... = 0 / gcd 0 m     : zero_mul
      ... = 0                 : nat.zero_div

theorem lcm_zero_right (m : ℕ) : lcm m 0 = 0 := !lcm.comm ▸ !lcm_zero_left

theorem lcm_one_left (m : ℕ) : lcm 1 m = m :=
calc
  lcm 1 m = 1 * m / gcd 1 m : rfl
      ... = m / gcd 1 m     : one_mul
      ... = m / 1           : gcd_one_left
      ... = m                 : nat.div_one

theorem lcm_one_right (m : ℕ) : lcm m 1 = m := !lcm.comm ▸ !lcm_one_left

theorem lcm_self (m : ℕ) : lcm m m = m :=
have H : m * m / m = m, from
  by_cases_zero_pos m !nat.div_zero (take m, assume H1 : m > 0, !nat.mul_div_cancel H1),
calc
  lcm m m = m * m / gcd m m : rfl
      ... = m * m / m       : gcd_self
      ... = m                 : H

theorem dvd_lcm_left (m n : ℕ) : m ∣ lcm m n :=
have H : lcm m n = m * (n / gcd m n), from nat.mul_div_assoc _ !gcd_dvd_right,
dvd.intro H⁻¹

theorem dvd_lcm_right (m n : ℕ) : n ∣ lcm m n :=
!lcm.comm ▸ !dvd_lcm_left

theorem gcd_mul_lcm (m n : ℕ) : gcd m n * lcm m n = m * n :=
eq.symm (nat.eq_mul_of_div_eq_right (dvd.trans !gcd_dvd_left !dvd_mul_right) rfl)

theorem lcm_dvd {m n k : ℕ} (H1 : m ∣ k) (H2 : n ∣ k) : lcm m n ∣ k :=
or.elim (eq_zero_or_pos k)
  (assume kzero : k = 0, !kzero⁻¹ ▸ !dvd_zero)
  (assume kpos : k > 0,
    have mpos : m > 0, from pos_of_dvd_of_pos H1 kpos,
    have npos : n > 0, from pos_of_dvd_of_pos H2 kpos,
    have gcd_pos : gcd m n > 0, from !gcd_pos_of_pos_left mpos,
    obtain p (km : k = m * p), from exists_eq_mul_right_of_dvd H1,
    obtain q (kn : k = n * q), from exists_eq_mul_right_of_dvd H2,
    have ppos : p > 0, from pos_of_mul_pos_left (km ▸ kpos),
    have qpos : q > 0, from pos_of_mul_pos_left (kn ▸ kpos),
    have H3 : p * q * (m * n * gcd p q) = p * q * (gcd m n * k), from
    calc
      p * q * (m * n * gcd p q)
            = m * p * (n * q * gcd p q)       : by rewrite [*mul.assoc, *mul.left_comm q,
                                                             mul.left_comm p]
        ... = k * (k * gcd p q)               : by rewrite [-kn, -km]
        ... = k * gcd (k * p) (k * q)         : by rewrite gcd_mul_left
        ... = k * gcd (n * q * p) (m * p * q) : by rewrite [-kn, -km]
        ... = k * (gcd n m * (p * q))         : by rewrite [*mul.assoc, mul.comm q, gcd_mul_right]
        ... = p * q * (gcd m n * k)           : by rewrite [mul.comm, mul.comm (gcd n m), gcd.comm,
                                                             *mul.assoc],
    have H4 : m * n * gcd p q = gcd m n * k,
      from !eq_of_mul_eq_mul_left (mul_pos ppos qpos) H3,
    have H5 : gcd m n * (lcm m n * gcd p q) = gcd m n * k,
      from !mul.assoc ▸ !gcd_mul_lcm⁻¹ ▸ H4,
    have H6 : lcm m n * gcd p q = k,
      from !eq_of_mul_eq_mul_left gcd_pos H5,
    dvd.intro H6)

theorem lcm.assoc (m n k : ℕ) : lcm (lcm m n) k = lcm m (lcm n k) :=
dvd.antisymm
  (lcm_dvd
    (lcm_dvd !dvd_lcm_left (dvd.trans !dvd_lcm_left !dvd_lcm_right))
    (dvd.trans !dvd_lcm_right !dvd_lcm_right))
  (lcm_dvd
    (dvd.trans !dvd_lcm_left !dvd_lcm_left)
    (lcm_dvd (dvd.trans !dvd_lcm_right !dvd_lcm_left) !dvd_lcm_right))

/- coprime -/

definition coprime [reducible] (m n : ℕ) : Prop := gcd m n = 1

lemma gcd_eq_one_of_coprime {m n : ℕ} : coprime m n → gcd m n = 1 :=
λ h, h

theorem coprime_swap {m n : ℕ} (H : coprime n m) : coprime m n :=
!gcd.comm ▸ H

theorem dvd_of_coprime_of_dvd_mul_right {m n k : ℕ} (H1 : coprime k n) (H2 : k ∣ m * n) : k ∣ m :=
have H3 : gcd (m * k) (m * n) = m, from
  calc
    gcd (m * k) (m * n) = m * gcd k n : gcd_mul_left
                    ... = m * 1       : H1
                    ... = m           : mul_one,
have H4 : (k ∣ gcd (m * k) (m * n)), from dvd_gcd !dvd_mul_left H2,
H3 ▸ H4

theorem dvd_of_coprime_of_dvd_mul_left {m n k : ℕ} (H1 : coprime k m) (H2 : k ∣ m * n) : k ∣ n :=
dvd_of_coprime_of_dvd_mul_right H1 (!mul.comm ▸ H2)

theorem gcd_mul_left_cancel_of_coprime {k : ℕ} (m : ℕ) {n : ℕ} (H : coprime k n) :
   gcd (k * m) n = gcd m n :=
have H1 : coprime (gcd (k * m) n) k, from
  calc
    gcd (gcd (k * m) n) k
         = gcd (k * gcd 1 m) n : by rewrite [-gcd_mul_left, mul_one, gcd.comm, gcd.assoc]
     ... = 1                   : by rewrite [gcd_one_left, mul_one, ↑coprime at H, H],
dvd.antisymm
  (dvd_gcd (dvd_of_coprime_of_dvd_mul_left H1 !gcd_dvd_left) !gcd_dvd_right)
  (dvd_gcd (dvd.trans !gcd_dvd_left !dvd_mul_left) !gcd_dvd_right)

theorem gcd_mul_right_cancel_of_coprime (m : ℕ) {k n : ℕ} (H : coprime k n) :
   gcd (m * k) n = gcd m n :=
!mul.comm ▸ !gcd_mul_left_cancel_of_coprime H

theorem gcd_mul_left_cancel_of_coprime_right {k m : ℕ} (n : ℕ) (H : coprime k m) :
   gcd m (k * n) = gcd m n :=
!gcd.comm ▸ !gcd.comm ▸ !gcd_mul_left_cancel_of_coprime H

theorem gcd_mul_right_cancel_of_coprime_right {k m : ℕ} (n : ℕ) (H : coprime k m) :
   gcd m (n * k) = gcd m n :=
!gcd.comm ▸ !gcd.comm ▸ !gcd_mul_right_cancel_of_coprime H

theorem coprime_div_gcd_div_gcd {m n : ℕ} (H : gcd m n > 0) :
  coprime (m / gcd m n) (n / gcd m n) :=
calc
  gcd (m / gcd m n) (n / gcd m n) = gcd m n / gcd m n : gcd_div !gcd_dvd_left !gcd_dvd_right
     ... = 1 : nat.div_self H

theorem not_coprime_of_dvd_of_dvd {m n d : ℕ} (dgt1 : d > 1) (Hm : d ∣ m) (Hn : d ∣ n) :
  ¬ coprime m n :=
assume co : coprime m n,
have d ∣ gcd m n, from dvd_gcd Hm Hn,
have d ∣ 1, by rewrite [↑coprime at co, co at this]; apply this,
have d ≤ 1, from le_of_dvd dec_trivial this,
show false, from not_lt_of_ge `d ≤ 1` `d > 1`

theorem exists_coprime {m n : ℕ} (H : gcd m n > 0) :
  exists m' n', coprime m' n' ∧ m = m' * gcd m n ∧ n = n' * gcd m n :=
have H1 : m = (m / gcd m n) * gcd m n, from (nat.div_mul_cancel !gcd_dvd_left)⁻¹,
have H2 : n = (n / gcd m n) * gcd m n, from (nat.div_mul_cancel !gcd_dvd_right)⁻¹,
exists.intro _ (exists.intro _ (and.intro (coprime_div_gcd_div_gcd H) (and.intro H1 H2)))

theorem coprime_mul {m n k : ℕ} (H1 : coprime m k) (H2 : coprime n k) : coprime (m * n) k :=
calc
  gcd (m * n) k = gcd n k : !gcd_mul_left_cancel_of_coprime H1
            ... = 1       : H2

theorem coprime_mul_right {k m n : ℕ} (H1 : coprime k m) (H2 : coprime k n) : coprime k (m * n) :=
coprime_swap (coprime_mul (coprime_swap H1) (coprime_swap H2))

theorem coprime_of_coprime_mul_left {k m n : ℕ} (H : coprime (k * m) n) : coprime m n :=
have H1 : (gcd m n ∣ gcd (k * m) n), from !gcd_dvd_gcd_mul_left,
eq_one_of_dvd_one (H ▸ H1)

theorem coprime_of_coprime_mul_right {k m n : ℕ} (H : coprime (m * k) n) : coprime m n :=
coprime_of_coprime_mul_left (!mul.comm ▸ H)

theorem coprime_of_coprime_mul_left_right {k m n : ℕ} (H : coprime m (k * n)) : coprime m n :=
coprime_swap (coprime_of_coprime_mul_left (coprime_swap H))

theorem coprime_of_coprime_mul_right_right {k m n : ℕ} (H : coprime m (n * k)) : coprime m n :=
coprime_of_coprime_mul_left_right (!mul.comm ▸ H)

theorem comprime_one_left : ∀ n, coprime 1 n :=
λ n, !gcd_one_left

theorem comprime_one_right : ∀ n, coprime n 1 :=
λ n, !gcd_one_right

theorem exists_eq_prod_and_dvd_and_dvd {m n k : nat} (H : k ∣ m * n) :
  ∃ m' n', k = m' * n' ∧ m' ∣ m ∧ n' ∣ n :=
or.elim (eq_zero_or_pos (gcd k m))
 (assume H1 : gcd k m = 0,
    have H2 : k = 0, from eq_zero_of_gcd_eq_zero_left H1,
    have H3 : m = 0, from eq_zero_of_gcd_eq_zero_right H1,
    have H4 : k = 0 * n, from H2 ⬝ !zero_mul⁻¹,
    have H5 : 0 ∣ m, from H3⁻¹ ▸ !dvd.refl,
    have H6 : n ∣ n, from !dvd.refl,
    exists.intro _ (exists.intro _ (and.intro H4 (and.intro H5 H6))))
  (assume H1 : gcd k m > 0,
    have H2 : gcd k m ∣ k, from !gcd_dvd_left,
    have H3 : k / gcd k m ∣ (m * n) / gcd k m, from nat.div_dvd_div H2 H,
    have H4 : (m * n) / gcd k m = (m / gcd k m) * n, from
      calc
        m * n / gcd k m = n * m / gcd k m   : mul.comm
                      ... = n * (m / gcd k m) : !nat.mul_div_assoc !gcd_dvd_right
                      ... = m / gcd k m * n   : mul.comm,
    have H5 : k / gcd k m ∣ (m / gcd k m) * n, from H4 ▸ H3,
    have H6 : coprime (k / gcd k m) (m / gcd k m), from coprime_div_gcd_div_gcd H1,
    have H7 : k / gcd k m ∣ n, from dvd_of_coprime_of_dvd_mul_left H6 H5,
    have H8 : k = gcd k m * (k / gcd k m), from (nat.mul_div_cancel' H2)⁻¹,
    exists.intro _ (exists.intro _ (and.intro H8 (and.intro !gcd_dvd_right H7))))

end nat

\end{verbatim}

\chapter{Formal Pre- and Post-conditions}
In Boldo\cite{Bold11} we find an effort to verify floating point software
using preconditions, postconditions, and assertions. Quoting:

``These conjectures can be described formally by annotations as follows.
\begin{verbatim}
/*@ requires \abs(x) <= 0x1p-5;
  @ ensures \abs(\result - \cos(x)) <= 0x1p-23;
  @*/
float my_cosine(float x) {
  //@ assert \abs(1.0 - x*x*0.5 - \cos(x)) < 0x1p-24;
  return 1.0f - x * x * 0.5f;
}
\end{verbatim}
The {\sl precondition}, introduced by {\bf requires}, states that we expect
argument {\sl x} in the interval [-1/32; 1/32]. The {\sl postcondition},
introduced by {\bf ensures}, states that the distance between the value
returned by the function, denoted by the keyword {\bf \verb|\result|}, and the
model of the program, which is here the true mathematical cosine function
denoted by {\bf \verb|\cos|} in ACSL, is not greater than $2^{-23}$. It is
important to notice that in annotations the operators like $+$ or $*$
denote operations on real numbers and not on floating-point numbers. In
particular, there is no rounding error and no overflow in annotations,
unlike in the early Leavens' proposal. The C variables of type {\tt float},
like {\tt x} and {\tt \verb|\result|} in this example, are interpreted as the
real number they represent. Thus, the last annotation, given as an
assertion inside the code, is a way to make explicit the reasoning we made
above, making the total error the sum of the method error and the rounding
error: it states that the method error is less than $2^{-24}$. Again, it
is thanks to the choice of having exact operations in the annotations that
we are able to state a property of the method error.''

In Boldo\cite{Bold07,Bold07a} we find 'search in an array' annotated:
\begin{verbatim}
/*@ requires \valid_range(t,0,n-1)
  @ ensures
  @  (0 <= \result < n => t[\result] == v) &&
  @  (\result == n =>
  @    \forall int i; 0 <= i < n => t[i] != v) */
int index(int t[], int n, int v) {
  int i = 0;
  /*@ invariant 0 <= i &&
    @  \forall int k; 0 <= k <i => t[k] != v
    @ variant n - i */
  while (i < n) {
    if (t[i] == v) break;
    i++;
  }
  return i;
}
\end{verbatim}

\chapter{Types and Signatures}

We need to start from a base of the existing types in Common Lisp,
eventually providing Axiom combinations or specializations.
Common Lisp has these standard type specifier symbols.
\begin{center}
{\bf Common Lisp Type Hierarchy}\cite{Pfei12}\\
\includegraphics[scale=0.5]{ps/v13cltypehierarchy.eps}
\end{center}

Axiom adds these types:
\begin{itemize}
\item Command = String
\end{itemize}

\chapter{COQ nat vs Axiom NNI}

COQ's nat domain includes a proof of GCD.

We would like to show an isomorphism between types in Coq and
types in Axiom. Having such an isomorphism will make lemmas
available and simplify future proofs.

Note that Coq's {\tt nat} domain stops at O (a symbolic 0)
as does Axiom's NNI. The Axiom interpreter will promote a
subtraction to Integer whereas Coq will not.

COQ's nat domain\cite{COQnat} is

\subsection{Library Coq.Init.Nat}

\begin{verbatim}
Require Import Notations Logic Datatypes.

Local Open Scope nat_scope.
\end{verbatim}

{\bf Peano natural numbers, definitions of operations}

This file is meant to be used as a whole module, without importing it,
leading to qualified definitions (e.g. Nat.pred)

\begin{verbatim}
Definition t := nat.
\end{verbatim}

{\bf Constants}

\begin{verbatim}
Definition zero := 0.
Definition one := 1.
Definition two := 2.
\end{verbatim}

{\bf Basic operations}

\begin{verbatim}
Definition succ := S.

Definition pred n :=
  match n with
    | 0 => n
    | S u => u
  end.

Fixpoint add n m :=
  match n with
  | 0 => m
  | S p => S (p + m)
  end

where "n + m" := (add n m) : nat_scope.

Definition double n := n + n.

Fixpoint mul n m :=
  match n with
  | 0 => 0
  | S p => m + p * m
  end

where "n * m" := (mul n m) : nat_scope.

\end{verbatim}

Note that Axiom's NNI domain will be automatically promoted to Integer
when the subtraction result is negative. Coq returns O when this occurs.

\begin{verbatim}
Truncated subtraction: n-m is 0 if n<=m

Fixpoint sub n m :=
  match n, m with
  | S k, S l => k - l
  | _, _ => n
  end

where "n - m" := (sub n m) : nat_scope.

\end{verbatim}
{\bf Comparisons}
\begin{verbatim}

Fixpoint eqb n m : bool :=
  match n, m with
    | 0, 0 => true
    | 0, S _ => false
    | S _, 0 => false
    | S n', S m' => eqb n' m'
  end.

Fixpoint leb n m : bool :=
  match n, m with
    | 0, _ => true
    | _, 0 => false
    | S n', S m' => leb n' m'
  end.

Definition ltb n m := leb (S n) m.

Infix "=?" := eqb (at level 70) : nat_scope.
Infix "<=?" := leb (at level 70) : nat_scope.
Infix "<?" := ltb (at level 70) : nat_scope.

Fixpoint compare n m : comparison :=
  match n, m with
   | 0, 0 => Eq
   | 0, S _ => Lt
   | S _, 0 => Gt
   | S n', S m' => compare n' m'
  end.

Infix "?=" := compare (at level 70) : nat_scope.

\end{verbatim}
{\bf Minimum, maximum}
\begin{verbatim}

Fixpoint max n m :=
  match n, m with
    | 0, _ => m
    | S n', 0 => n
    | S n', S m' => S (max n' m')
  end.

Fixpoint min n m :=
  match n, m with
    | 0, _ => 0
    | S n', 0 => 0
    | S n', S m' => S (min n' m')
  end.

\end{verbatim}
{\bf Parity tests}
\begin{verbatim}

Fixpoint even n : bool :=
  match n with
    | 0 => true
    | 1 => false
    | S (S n') => even n'
  end.

Definition odd n := negb (even n).

\end{verbatim}
{\bf Power}
\begin{verbatim}

Fixpoint pow n m :=
  match m with
    | 0 => 1
    | S m => n * (n^m)
  end

where "n ^ m" := (pow n m) : nat_scope.

\end{verbatim}
{\bf Euclidean division}
\begin{verbatim}
This division is linear and tail-recursive. In divmod, y is the 
predecessor of the actual divisor, and u is y minus the real remainder

Fixpoint divmod x y q u :=
  match x with
    | 0 => (q,u)
    | S x' => match u with
                | 0 => divmod x' y (S q) y
                | S u' => divmod x' y q u'
              end
  end.

Definition div x y :=
  match y with
    | 0 => y
    | S y' => fst (divmod x y' 0 y')
  end.

Definition modulo x y :=
  match y with
    | 0 => y
    | S y' => y' - snd (divmod x y' 0 y')
  end.

Infix "/" := div : nat_scope.
Infix "mod" := modulo (at level 40, no associativity) : nat_scope.

\end{verbatim}
{\bf Greatest common divisor} 
\begin{verbatim}
We use Euclid algorithm, which is normally not structural, but Coq is 
now clever enough to accept this (behind modulo there is a subtraction, 
which now preserves being a subterm)

Fixpoint gcd a b :=
  match a with
   | O => b
   | S a' => gcd (b mod (S a')) (S a')
  end.

\end{verbatim}
{\bf Square}
\begin{verbatim}

Definition square n := n * n.

\end{verbatim}
{\bf Square root}
\begin{verbatim} 

The following square root function is linear (and tail-recursive). 
With Peano representation, we can't do better. For faster algorithm, 
see Psqrt/Zsqrt/Nsqrt...  We search the square root of 
n = k + p^2 + (q - r) with q = 2p and 0<=r<=q. We start with
p=q=r=0, hence looking for the square root of n = k. Then we
progressively decrease k and r. When k = S k' and r=0, it means we can
use (S p) as new sqrt candidate, since (S k')+p^2+2p = k'+(S
p)^2. When k reaches 0, we have found the biggest p^2 square contained
in n, hence the square root of n is p.

Fixpoint sqrt_iter k p q r :=
  match k with
    | O => p
    | S k' => match r with
                | O => sqrt_iter k' (S p) (S (S q)) (S (S q))
                | S r' => sqrt_iter k' p q r'
              end
  end.

Definition sqrt n := sqrt_iter n 0 0 0.

\end{verbatim}
{\bf Log2} 
\begin{verbatim}
This base-2 logarithm is linear and tail-recursive.  In
log2_iter, we maintain the logarithm p of the counter q, while r is
the distance between q and the next power of 2, more precisely q + S r
= 2^(S p) and r<2^p. At each recursive call, q goes up while r goes
down. When r is 0, we know that q has almost reached a power of 2, and
we increase p at the next call, while resetting r to q.  Graphically
(numbers are q, stars are r) :

                    10
                  9
                8
              7   *
            6       *
          5           ...
        4
      3   *
    2       *
  1   *       *
0   *   *       *

We stop when k, the global downward counter reaches 0. At that moment,
q is the number we're considering (since k+q is invariant), and p its
logarithm.

Fixpoint log2_iter k p q r :=
  match k with
    | O => p
    | S k' => match r with
                | O => log2_iter k' (S p) (S q) q
                | S r' => log2_iter k' p (S q) r'
              end
  end.

Definition log2 n := log2_iter (pred n) 0 1 0.

Iterator on natural numbers

Definition iter (n:nat) {A} (f:A->A) (x:A) : A :=
 nat_rect (fun _ => A) x (fun _ => f) n.

Bitwise operations We provide here some bitwise operations for unary
numbers. Some might be really naive, they are just there for
fullfiling the same interface as other for natural representations. As
soon as binary representations such as NArith are available, it is
clearly better to convert to/from them and use their ops.

Fixpoint div2 n :=
  match n with
  | 0 => 0
  | S 0 => 0
  | S (S n') => S (div2 n')
  end.

Fixpoint testbit a n : bool :=
 match n with
   | 0 => odd a
   | S n => testbit (div2 a) n
 end.

Definition shiftl a := nat_rect _ a (fun _ => double).
Definition shiftr a := nat_rect _ a (fun _ => div2).

Fixpoint bitwise (op:bool->bool->bool) n a b :=
 match n with
  | 0 => 0
  | S n' =>
    (if op (odd a) (odd b) then 1 else 0) +
    2*(bitwise op n' (div2 a) (div2 b))
 end.

Definition land a b := bitwise andb a a b.
Definition lor a b := bitwise orb (max a b) a b.
Definition ldiff a b := bitwise (fun b b' => andb b (negb b')) a a b.
Definition lxor a b := bitwise xorb (max a b) a b.
\end{verbatim}

\chapter{Binary Power in COQ by Casteran and Sozeau}

From Casteran and Sozeau\cite{Cast16}:

\begin{verbatim}
(* About integer powers (monomorphic version) *)

Set Implicit Arguments.
Require Import ZArith.
Require Import Div2.
Require Import Program.
Open Scope Z_scope.
\end{verbatim}

Let us consider a simple arithmetic operation: raising some integer $x$
to the $n$-th power, where $n$ is a natural number. The following
function definition is a direct translation of the mathematical
concept:
\begin{verbatim}
Fixpoint power (a:Z)(n:nat) :=
  match n with 0%nat => 1
             | S p =>  a * power a p
  end.

Eval vm_compute in power 2 40.
 = 1099511627776 : Z
\end{verbatim}

This definition can be considered as a very naive way of programming,
since computing $x^n$ requires $n$ multiplications. Nevertheless, this
definition is very simple to read, and everyone can admit that it is
correct with respect to the mathematical definition. Thus, we can consider
it as a {\sl specification}: when we write more efficient but less readable
functions for exponentiation, we should be able to prove their correctness
by proving in Coq their equivalence with the naive power function.

The following function allows one to compute $x^n$, with a number of
multiplications proportional to $\log_2(n)$:
\begin{verbatim}
Program
Fixpoint binary_power_mult (acc x:Z) (n:nat) {measure (fun i=>i) n} : Z
  (* acc * (power x n) *) :=
  match n with 
    | 0%nat => acc
    | _ => if Even.even_odd_dec n
           then binary_power_mult acc (x * x) (div2 n)
           else binary_power_mult (acc * x) (x * x) (div2 n)
  end.

Solve Obligations with program_simpl; intros; apply lt_div2; auto with arith.

Definition binary_power (x:Z)(n:nat) := binary_power_mult 1 x n.

Eval vm_compute in binary_power 2 40.
 = 1099511627776 : Z

Goal binary_power 2 234 = power 2 234.
reflexivity.
Qed.
\end{verbatim}

We want now to {\sl prove} {\tt binary\_power's} correctness, i.e. that
this function and the naive {\tt power} function are pointwise equivalent.

Proving this equivalence in Coq may require a lot of work. Thus it is not
worth at all writing a proof dedicated only to powers of integers. In fact, 
the correctness of {\tt binary\_power} with respect to {\tt power} holds
in any structure composed of an associative binary operation on some domain, 
that admits a neutral element. For instance, we can compute powers of square
matrices using the most efficient of both algorithms.

Thus, let us throw away our previous definition, and try to define them in
a more generic framework.

\section{On Monoids}

{\bf Definition 2.1} {\sl A monoid is a mathematical structure composed of}
\begin{itemize}
\item {\sl a carrier A}
\item {\sl a binary, associative operation} $\circ$ {\sl on A}
\item {\sl a neutral element} $1 \in A$ {\sl for} $\circ$
\end{itemize}

Such a mathematical structure can be defined in Coq as a type class.
\cite{Soze08}. In the following definition, parameterized by a type $A$
(implicit), a binary operation {\tt dot} and a neutral element {\tt unit},
three fields describe the properties that {\tt dot} and {\tt unit} must
satisfy.
\begin{verbatim}
Class Monoid {A:Type}(dot : A -> A -> A)(one : A) : Prop := {
  dot_assoc : forall x y z:A, dot x (dot y z) = dot (dot x y) z;
  unit_left : forall x, dot one x = x;
  unit_right : forall x, dot x one = x }.
\end{verbatim}

Note that other definitions could have been given for representing this
mathematical structure. 

From an implementational point of view, such a type class is just a record
type, i.e. an inductive type with a single construtor {\tt Build\_Monoid}
\begin{verbatim}
Print Monoid.

Record Monoid (A:Type)(dot : A -> A -> A)(one : A) : Prop := Build_Monoid
  { dot_assoc : forall x y z:A, dot x (dot y z) = dot (dot x y) z;
    one_left : forall x, dot one x = x;
    one_right : forall x, dot x one = x }

For Monoid: Argument A is implicit and maximally inserted
For Build_Monoid: Argument A is implicit
For Monoid: Argument scopes are [type_scope _ _]
For Build_Monoid: Argument scopes are [type_scope _ _ _ _ _]
\end{verbatim}

Nevertheless, implementation of type classes by M. Sozeau provides several
specific tools --- dedicated tactics for instance --, and we advise the
reader not to replace the {\tt Class} keyword with {\tt Record} or
{\tt Inductive}.

WIth the command {\tt About}, we can see the polymorphic type of the 
fields of the class {\tt Monoid}:
\begin{verbatim}
About one_left

one_left:
forall (A : Type) (dot : A -> A -> A) (one : A),
Monoid dot one -> forall x : A, dot one x = x

Arguments A, dot, one, Monoid are implicit and maximally inserted
Argument scopes are [type_scope _ _ _ _]
one_left is transparent
\end{verbatim}

\subsection{Classes and Instances}

Members of a given class are called {\sl instances} of this class. 
Instances are defined to the Coq system through the {\tt Instance}
keyword. Our first example is a definition of the monoid structure
on the set $\mathbb{Z}$ of integers, provided with integer multiplication,
with 1 as the neutral element. Thus we give these parameters to the
{\tt Monoid} class (note that $\mathbb{Z}$ is implicitly given).
\begin{verbatim}
Instance ZMult : Monoid Zmult 1
\end{verbatim}

For this instance to be created, we need to prove that the binary
operation {\tt Zmult} is associative and admits 1 as the neutral
element. Applying the constructor {\tt Build\_Monoid} -- for instance
with the tactic {\tt split} -- generates three subgoals.
\begin{verbatim}
split.
3 subgoals
 =================================================
  forall x y z : Z, x * (y * z) = x * y * z

subgoal 2 is:
  forall x : Z, 1 * x = x
subgoal 3 is:
  forall x : Z, x * 1 = x
\end{verbatim}

Each subgoal is easily solved by {\tt intros; ring}.

When the proof is finished, we register our instance with a simple {\tt Qed}.
Note that we used {\tt Qed} because we consider a class of sort {\tt Prop}.
In some cases where instances must store some information constants,
ending an instance construction with {\tt Defined} may be necessary.

\begin{verbatim}
Check Zmult.
ZMult : Monoid Zmult 1
\end{verbatim}

We explained on the preceding page why it is better to use the {\tt Class}
keyword than {\tt Record} or {\tt Inductive}. For the same reason, the 
definition of an instance of some class should be written using 
{\tt Instance} and not {\tt Lemma}, {\tt Theorem}, {\tt Example}, etc.
nor {\tt Definition}.

\subsection{A generic definition of {\tt power}}

We are now able to give a definition of the function {\tt power} than can
be applied with any instance of class {\tt Monoid}:

A first definition could be
\begin{verbatim}
Fixpoint power {A:Type}{dot:A->A->A}{one:A}{M: Monoid dot one}
               (a:A)(n:nat) :=
  match n with 0:nat => one
             | S p => dot a (power a p)
  end.

Compute power 2 10.
 = 1024 : Z
\end{verbatim}

Happily, we can make the declaration of the three first arguments
implicit, by using the {\tt Generalizable Variables} command:
\begin{verbatim}
Reset power.

Generalizable Variables A dot one.

Fixpoint power `{M: Monoid A dot one}(a:A)(n:nat) :=
  match n with 0%nat => one
             | S p => dot a (power a p)
  end.

Compute power 2 10.
 = 1024 : Z
\end{verbatim}

The variables {\tt A dot one} appearing in the binder for {\tt M} are
implicitly bound before the binder for {\tt M} and their types are
inferred from the {\tt Monoid A dot one} type. This syntactic sugar
helps abbreviate bindings for classes with parameters. The resulting
internal Coq term is exactly the same as the first definition above.

\subsection{Instance Resolution}

The attentitive reader has certainly noticed that in the last
computation, the binary operation {\tt Zmult} and the neutral element 1
need not to be given explicitly. The mechanism that allows Coq to infer
all the arguments needed by the {\tt power} funtion to be applied is called
{\sl instance resolution}.

In order to understand how it operates, let's have a look at {\tt power}'s
type:
\begin{verbatim}
About power.
power :
forall (A : Type) (dot : A -> A -> A) (one : A),
Monoid dot one -> A -> nat -> A

Arguments A, dot, one, M are implicit and maximally inserted

Compute power 2 100.
= 1267650600228229401496703205376 : Z

Set Printing Implicit.
Check power 2 100.
@power Z Zmult 1 Zmult 2 100 : Z
Unset Printing Implicit.
\end{verbatim}

We see that the {\sl instance} {\tt ZMult} has been inferred from the
type of 2. We are in the simple case where only one monoid of carrier {\tt Z}
has been declared as an instance of the {\tt Monoid} class.

The implementation of type classes in Coq can retrieve the instance
{\tt ZMult} from the type {\tt Z}, then filling the arguments {\tt ZMult}
and 1 from {\tt ZMult}'s definition.

\section{More Monoids}
\subsection{Matrices over some ring}

We all know that multiplication of square matrices is associative and
admits identity matrices as neutral elements. For simplicity's sake let us
restrict our study to $2\times 2$ matrices over some ring.

We first load the {\tt Ring} library, then open a section with some useful
declarations and notations. 
\begin{verbatim}
Require Import Ring.

Section matrices.
  Variables (A:Type)
            (zero one : A)
            (plus mult minus : A -> A -> A)
            (sym : A -> A).
  Notation "0" := zero.
  Notation "1" := one.
  Notation "x + y" := (plus x y).
  Notation "x * y" := (mult x y).

  Variable rt : ring_theory zero one plus mult minus sym (@eq A).

  Add Ring Aring : rt.
\end{verbatim}

We can now define a carrier type for $2\times 2$-matrices, as well as
matrix multiplication and the identiy matrix.
\begin{verbatim}
Structure M2 : Type := {c00 : A; c01 : A; c10 : A; c11 : A}.

Definition Id2 : M2 := Build_M2 1 0 0 1.

Definition M2_mult (m m':M2) : M2 :=
  Build_M2 (c00 m * c00 m' + c01 m * c10 m')
           (c00 m * c01 m' + c01 m * c11 m')
           (c10 m * c00 m' + c11 m * c10 m')
           (c10 m * c01 m' + c11 m * c11 m').
\end{verbatim}

As for multiplication of integers, we can now define an instance of
{\tt Monoid} for the type {\tt M2}.
\begin{verbatim}
Global Instance M2_Monoid : Monoid (M2_mult plus mult) (Id2 0 1).
split.
destruct x; destruct y; destruct z; simpl.
unfold M2_mult. apply M2_eq_intros; simpl; ring.
destruct x; simpl;
unfold M2_mult; apply M2_eq_intros; simpl; ring.
destruct x; simpl;
unfold M2_mult; apply M2_eq_intros; simpl; ring.
Qed.

End matrices.
\end{verbatim}

We want now to play with $2\times 2$ matrices over $\mathbb{Z}$.
We declare an instance {\bf M2Z} for this purpose, and can use directly
the function {\tt power}.
\begin{verbatim}
Instance M2Z : Monoid _ _ := M2_Monoid Zth.

Compute power (Build_M2 1 1 1 0) 40.
 = {|
     c00 := 165580141;
     c01 := 102334155;
     c10 := 102334155;
     c11 := 63245986 |}
   : M2 Z

Definition fibonacci (n:nat) :=
  C00 (power (Build_M2 1 1 1 0) n).

Compute fibonacci 20.
 = 10946
 :Z
\end{verbatim}

\section{Reasoning within a Type Class}

We are now able to consider again the equivalence between two functions
for computing powers. Let use define the binary algorithm for any monoid.

First, we define an auxiliary function. We use the {\tt Program} 
extension to define an efficient version of exponentiation using an
accumulator. The function is defined by well-founded recursion on 
the exponent $n$.
\begin{verbatim}
Function binary_power_mult (A:Type) (dot:A->A->A) (one:A)
     (M: @Monoid A dot one) (acc x:A)(n:nat){measure (fun i=>i) n} : A
  (* acc * (x ** n) *) :=
  match n with 
    | 0%nat => acc
    | _ => if Even.even_odd_dec n
             then binary_power_mult _ acc (dot x x) (div2 n)
             else binary_power_mult _ (dot acc x) (dot x x) (div2 n)
  end.
intros; apply lt_div2; auto with arith.
intros; apply l2_div2; auto with arith.
Defined.

Definition binary_power `{M:Monoid} x n := binary_power_mult M one x n.

Compute binary_power 2 100.
 = 1267650600228229401496703205376 : Z
\end{verbatim}

\subsection{The Equivalence Proof}

The proof of equivalence between {\tt power} and {\tt binary\_power}
is quite long, and can be split in several lemmas. Thus, it is useful
to open a section, in which we fix some arbitrary monoid M. Such a
declaration is made with the {\tt Context} command, which can be 
considered as a version of {\tt Variables} for declaring arbitrary
instances of a given class.
\begin{verbatim}

Section About_power.

Require Import Arith.
 Context `(M:Monoid A dot one ).
\end{verbatim}

It is good practice to define locally some specialized notations and 
tactics.
\begin{verbatim}

Ltac monoid_rw :=
    rewrite (@one_left A dot one M) || 
    rewrite (@one_right A dot one M)|| 
    rewrite (@dot_assoc A dot one M).


Ltac monoid_simpl := repeat monoid_rw.

Local Infix "*" := dot.
Local Infix "**" := power (at level 30, no associativity).
\end{verbatim}

\subsection{Some Useful Lemmas About {\tt power}}

We start by proving some well-known equalities about powers in a monoid.
Some of these equalities are integrated later in simplification tactics.

\begin{verbatim}
Lemma power_x_plus : forall x n p, x ** (n + p) =  x ** n *  x ** p.
Proof.
 induction n as [| p IHp];simpl.
   intros;  monoid_simpl;trivial.
  intro q;rewrite (IHp q);  monoid_simpl;trivial. 
Qed.

Ltac power_simpl := repeat (monoid_rw || rewrite <- power_x_plus).

Lemma power_commute : forall x n p,  
               x ** n * x ** p = x ** p * x ** n. 
Proof.
  intros x n p;power_simpl; rewrite (plus_comm n p);trivial.
Qed.

Lemma power_commute_with_x : forall x n ,  
        x * x ** n = x ** n * x.
Proof.
  induction n;simpl;power_simpl;trivial.
  repeat rewrite <- (@dot_assoc A dot one M); rewrite IHn; trivial.
Qed.

Lemma power_of_power : forall x n p,  (x ** n) ** p = x ** (p * n).
Proof.
   induction p;simpl;[| rewrite power_x_plus; rewrite IHp]; trivial.
Qed.

Lemma power_S : forall x n, x *  x ** n = x ** S n.
Proof. intros;simpl;auto. Qed.

Lemma sqr : forall x, x ** 2 =  x * x.
Proof.
 simpl;intros;monoid_simpl;trivial.
Qed.

Ltac factorize := repeat (
                rewrite <- power_commute_with_x ||
                rewrite  <- power_x_plus  ||
                rewrite <- sqr ||
                rewrite power_S ||
                rewrite power_of_power).

Lemma power_of_square : forall x n, (x * x) ** n = x ** n * x ** n.
  induction n;simpl;monoid_simpl;trivial.
  repeat rewrite dot_assoc;rewrite IHn; repeat rewrite dot_assoc.
 factorize; simpl;trivial.
Qed.
\end{verbatim}

\subsection{Final Steps}

We are now able to prove that the auxilliary function
{\tt binary\_power\_mult} satisfies its intuitive meaning. 
The proof uses well-founded induction and the lemmas proven in the
previous section.

\begin{verbatim}
Lemma binary_power_mult_ok :
  forall n a x,  binary_power_mult a x n = a * x ** n.
Proof.
  intro n; pattern n;apply lt_wf_ind.
  clear n; intros n Hn;   destruct n.
   intros;simpl; monoid_simpl;  trivial.
  intros; rewrite binary_power_mult_equation. 
     destruct (Even.even_odd_dec (S n)).
   rewrite Hn.   rewrite power_of_square;  factorize.
   pattern (S n) at 3;replace (S n) with (div2 (S n) + div2 (S n))%nat;auto.
   generalize (even_double _ e);simpl;auto. 
   apply lt_div2;auto with arith.
   rewrite Hn. 
  rewrite power_of_square ; factorize.
  pattern (S n) at 3;replace (S n) with (S (div2 (S n) + div2 (S n)))%nat;auto.

  rewrite <- dot_assoc; factorize;auto.
  generalize (odd_double _ o);intro H;auto.
  apply lt_div2;auto with arith.
Qed.
\end{verbatim}

Then the main theorem follows immediately:
\begin{verbatim}
Lemma binary_power_ok : forall (x:A) (n:nat), binary_power x n = x ** n.
Proof.
  intros n x;unfold binary_power;rewrite binary_power_mult_ok;
  monoid_simpl;auto.
Qed.
\end{verbatim}

\subsection{Discharging the Context}

It is time to close the section we opened for writing our proof of
equivalence. The theorem {\tt binary\_power\_ok} is now provided
with a universal quantification over all the parameters of any monoid.

\begin{verbatim}
End About_power.

About binary_power_ok.
binary_power_ok :
forall (A : Type) (dot : A -> A -> A) (one : A) (M : Monoid dot one)
  (x : A) (n : nat), binary_power x n = power x n

Arguments A, dot, one M are implicit and maximally inserted
Argument scopes are [type_scope _ _ _ _ nat_scope]
binary_power_ok is opaque
Expands to Constant Top.binary_power_ok

Check binary_power_ok 2 20.
binary_power_ok 2 20
  : binary_power 2 20 = power 2 20

Let Mfib := Build_M2  1 1 1 0.

Check  binary_power_ok Mfib 56.
binary_power_ok Mfib 56
  : binary_power Mfib 56 = power Mfib 56

\end{verbatim}

\subsection{Subclasses}

We could prove many useful equalities in the section {\tt about\_power}.
Nevertheless, we couldn't prove the equality $(xy)^n=x^ny^n$ because it
is false in general -- consider for instance the free monoid of strings,
or simply matrix multiplication. But this equality holds in every 
commutative (a.k.a abelian) monoid.

Thus we say that abelian monoids form a {\sl subclass} of the class of
monoids, and prove this equality in a context declaring an arbitrary 
instance of this subclass.

Structurally, we parameterize the new class {\tt Abelian\_Monoid} by an
arbitrary instance {\sl M} of {\tt Monoid}, and add a new field stating
the commutativity of {\tt dot}. Please keep in mind that we declared
{\tt A}, {\tt dot}, and {\tt one} as {\sl generalizable variables},
hence we can use the backquote symbol here.

\begin{verbatim}
Class Abelian_Monoid `(M:Monoid A dot one) := {
  dot_comm : forall x y, dot x y = dot y x}.
\end{verbatim}

A quick look at the representation of {\sl Abelian\_Monoid} as a record
type helps us understand how this class is implemented.
\begin{verbatim}
Print Abelian_Monoid.
Record Abelian_Monoid (A : Type) (dot : A -> A -> A)
  (one : A) (M : Monoid dot one) : Prop := Build_Abelian_Monoid
  {dot_comm : forall x y : A, dot x y = dot y x }

For Abelian_Monoid: Arguments A, dot, one are implicit and maximally inserted
For Build_Abelian_Monoid: Arguments A, dot, one are implicit
For Abelian_Monoid: Arguemnt scopes are [type_scope _ _ _]
For Build_Abelian_Monoid: Argument scopes are [type_scope _ _ _ _]
\end{verbatim}

For building an instance of {\tt Abelian\_Monoid} we can start from
{\tt ZMult}, the monoid on {\tt Z}, adding a proof that integer
multiplication is commutative.
\begin{verbatim}
Instance ZMult_Abelian : Abelian_Monoid ZMult.
split. 
 exact Zmult_comm.
Qed.
\end{verbatim}

We can now prove our equality by building an appropriate context. Note
that we can specify just the parameters of the monoid here in the binder
of the abelian monoid, an instance of monoid on those same parameters is
automatically generalized. Superclass parameters are automatically
generalized inside quote binders. Again, this is siimply syntactic sugar.
\begin{verbatim}
Section Power_of_dot.
 Context `{M: Monoid A} {AM:Abelian_Monoid M}.
 
Theorem power_of_mult : forall n x y, 
    power (dot x y)  n =  dot (power x  n) (power y n). 
Proof.
 induction n;simpl.
 rewrite one_left;auto.
 intros; rewrite IHn; repeat rewrite dot_assoc.
  rewrite <- (dot_assoc  x y (power x n)); rewrite (dot_comm y (power x n)).
 repeat rewrite dot_assoc;trivial.
Qed.

End Power_of_dot.

Check power_of_mult 3 4 5.
power_of_mult 3 4 5
     : power (4 * 5) 3 = power 4 3 * power 5 3
\end{verbatim}

\chapter{Other Ideas to Explore}
Computerising Mathematical Text\cite{Kama15} explores various ways of 
capturing mathematical reasoning. 

Chlipala\cite{Chli15} gives a pragmatic approach to COQ.

Medina-Bulo et al.\cite{Bulo04} gives a formal verification of
Buchberger's algorithm using ACL2 and Common Lisp.

Th\'ery\cite{Ther01} used COQ to check an implementation of Buchberger's
algorithm.

Pierce\cite{Pier15} has a Software Foundations course in COQ with
downloaded files in Pier15.tgz.

Spitters\cite{Spit11} Type Classes for Mathematics in Coq. Also see
\verb|http://www.eelis.net/research/math-classes/|

Mahboubi\cite{Mahb16} Mathematical Components. This book contains a
proof of the Euclidean algorithm using COQ.

Aczel\cite{Acze13} Homotopy Type Theory

\appendix

\chapter{The Global Environment}

Let $S$ be a set. Let $\circ$ be a binary operation. 
Let $+$ be an additive operation. Let $*$ be a multiplicative operation.

\begin{axiom}[Magma]
A {\bf Magma} is the
set $S$ with a {\bf closed binary operation} 
$S \circ S \rightarrow S$ such that
\[\forall a,b \in S \Rightarrow a \circ b \in S\].
\end{axiom}

\begin{axiom}[Semigroup]
A {\bf Semigroup} is a
{\bf Magma} with the operation $\circ$ that is {\bf associative} such that
\[\forall a,b,c \in S \Rightarrow (a \circ b) \circ c = a \circ (b \circ c)\]
\end{axiom}

\begin{axiom}[Abelian Semigroup]
An {\bf Abelian Semigroup} is a {\bf Semigroup} with the operation 
$\circ$ that is {\bf commutative} such that
\[\forall a,b \in S \Rightarrow a \circ b = b \circ a\]
\end{axiom}

\begin{axiom}[Monoid]
A {\bf Monoid} is a
{\bf Semigroup} with an {\bf identity element} $e \in S$ such that
\[\forall a \in S \Rightarrow e \circ a = a \circ e = a\]
\end{axiom}

\begin{axiom}[Group]
A {\bf Group} is a {\bf Monoid} with an {\bf inverse element} $b \in S$
and an identity element $i \in S$
such that
\[\forall a \in S~\exists b \in S \Rightarrow a \circ b = b \circ a = i\]
\end{axiom}

\begin{axiom}[Group Unique Identity]
A {\bf Group} has a {\bf unique identity element} $e \in S$ such that
\[\exists e \land \forall a,b \in S \land a \ne e \land b \ne e \Rightarrow 
a \circ b \ne e\]
\end{axiom}

\begin{axiom}[Group Unique Inverse]
A {\bf Group} has a {\bf unique inverse element} $i \in S$ such that
\[\exists i \land \forall a,b \in S \land a \ne i \land b \ne i \Rightarrow 
a \circ b \ne i\]
\end{axiom}

\begin{axiom}[Group Right Quotient]
A {\bf Group} has a {\bf Right Quotient} (right division) such that
\[x \circ a = b \Rightarrow x \circ a \circ a^{-1} = b \circ a^{-1} 
\Rightarrow x = b \circ a^{-1}\]
\end{axiom}

\begin{axiom}[Group Left Quotient]
A {\bf Group} has a {\bf Left Quotient} (left division) such that
\[a \circ x = b \Rightarrow a^{-1} \circ a \circ x = a^{-1} \circ b
\Rightarrow x = a^{-1} \circ b\]
\end{axiom}

\begin{axiom}[Abelian Group]
An {\bf Abelian Group} is a {\bf Group} with the operation $\circ$ that is
{\bf commutative} such that
\[\forall a,b \in S \Rightarrow a \circ b = b \circ a\]
\end{axiom}

\begin{axiom}[Abelian Group Quotient]
An {\bf Abelian Group} has a {\bf Quotient} (division) such that
\[a^{-1} \circ a \circ x = a \circ a^{-1} \circ x\]
\end{axiom}



\begin{axiom}[Euclidean Domain]
Let $R$ be an integral domain. Let $f$ be a function from
$R\backslash \{0\}$ to the NonNegativeInteger domain.
If $a$ and $b$ are in $R$ and $b$ is nonzero, then there are
$q$ and $r$ in $R$ such that $a=bq+r$ and either $r=0$ or
$f(r)<f(b)$ 
\end{axiom}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\cleardoublepage
\phantomsection
\addcontentsline{toc}{chapter}{Bibliography}
\bibliographystyle{axiom}
\bibliography{axiom}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\cleardoublepage
\phantomsection
\addcontentsline{toc}{chapter}{Index}
\printindex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}