Windmills of the Minds: An Algorithm for Fermat's Two Squares Theorem]Windmills of the Minds: An Algorithm for
Fermat's Two Squares Theorem

Hing Lun Chan

Australian National University

/Author (Hing Lun Chan)
/Title  (Windmills of the minds: an algorithm for Fermat's Two Squares Theorem)

/Subject (Theorem Proving)
/Keywords (number theory;algorithm;prime)

The two squares theorem of Fermat is a gem in number theory,
with a spectacular one-sentence “proof from the Book”.
Here is a formalisation of this proof, with an interpretation using windmill patterns.
The theory behind involves involutions on a finite set,
especially the parity of the number of fixed points in the involutions.
Starting as an existence proof that is non-constructive, there is an ingenious way to turn it into a constructive one. This gives an algorithm to compute the two squares by iterating the two involutions alternatively from a known fixed point.

<concept_desc>Theory of computation</concept_desc>
<concept_desc>Theory of computation Automated reasoning</concept_desc>

[100]Theory of computation
[500]Theory of computation Automated reasoning

§ INTRODUCTION

Fermat's two squares theorem, dated back to 1640, states that a prime $n$ that is one more than a multiple of $4$ can be uniquely expressed as a sum of odd and even squares (Section <ref>, Theorem <ref>).
Of the many proofs of this classical number theory result, this one-sentence proof by Zagier <cit.> caused a sensation in 1990:

The involution on the finite set

\[
S = \{(x,y,z) \in \mathbb{N}^{3} \mid \HOLinline{\HOLFreeVar{n}\;\HOLSymConst{=}\;\HOLFreeVar{x}\HOLSymConst{\ensuremath{\sp{2}}}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{4}\HOLSymConst{\ensuremath{}}\HOLFreeVar{y}\HOLSymConst{\ensuremath{}}\HOLFreeVar{z}} \}
\]

defined by

\begin{equation}
\label{eqn:zagier-map}
\HOLinline{(\HOLFreeVar{x}\HOLSymConst{,}\HOLFreeVar{y}\HOLSymConst{,}\HOLFreeVar{z})} \longmapsto
\begin{cases}
(x + 2 z,z,y - z - x) & \text{if }x < y - z \\
(2 y - x,y,x + z - y) & \text{if }y - z < x < 2y\\
(x - 2 y,x + z - y,y) & \text{if }x > 2y
\end{cases}
\end{equation}

has exactly one fixed point, so S is odd, and the involution defined by
$\HOLinline{(\HOLFreeVar{x}\HOLSymConst{,}\HOLFreeVar{y}\HOLSymConst{,}\HOLFreeVar{z})} \longmapsto \HOLinline{(\HOLFreeVar{x}\HOLSymConst{,}\HOLFreeVar{z}\HOLSymConst{,}\HOLFreeVar{y})}$
also has a fixed point.

Those who are perplexed by this multi-line sentence are not alone.

Even knowing involution, a self-inverse function, and fixed points, those values unchanged by a function, the proof is not obvious at a glance!

Listed as number 20 in Formalizing 100 Theorems <cit.>,
there are many formal proofs of this theorem.
Some are based on textbook proofs, others follow the ideas in Zagier's proof.
All show the existence of the two squares, only a few (Coq <cit.> and Lean <cit.>) include the uniqueness part.
Therefore, a formalisation of this one-sentence proof, in a constructive way, is an interesting exercise in theorem-proving.
As a bonus, the exercise is a path of discovery due to recent progress in understanding this proof.

As Don Zagier remarked after the one sentence, his proof was a condensed version of a 1984 proof by Roger Heath-Brown <cit.>, who in turn acknowledged prior work in number theory taken up by Joseph Liouville <cit.>.
This one-sentence proof invokes two involutions: the second one is obvious, but the first one in Equation (<ref>) has been called “black magic” <cit.>. The algebraic formulation of this involution has been given a geometric interpretation by Alexander Spivak <cit.> in 2007. These are the windmills (Section <ref>). They explain why the magic works, and suggest an interplay of the involutions to identify fixed points of each other. Moreover, this provides an algorithm to find the two squares in Fermat's theorem.
Thus the one-sentence proof can be made constructive, as elucidated by Zagier <cit.> in 2013.

§.§ Contribution

This paper gives the first formal proof of an algorithm to compute the two squares in Fermat's two squares theorem, by following a constructive version of Zagier's proof in HOL4.

As noted before, Zagier's proof has been formalised, in HOL Light <cit.>, in NASA PVS <cit.> and in Coq <cit.>, although not in this constructive form.

All the ideas used in this paper can be found in Shiu <cit.> and Zagier <cit.>.

The novel feature of this work is an elegant and pictorial approach for our formalisation.
The emphasis is in providing formal definitions and developing appropriate theories, not only for the present work, but also for supporting further work.

§.§ Overview

Major features in this formalisation are:

* the groundwork for Zagier's proof in Section <ref>,

* the two involutions for windmills in Section <ref>,

* the existence and uniqueness of two squares in Section <ref>,

* an algorithm to compute the two squares in Section <ref>,

* theories of involutions and iterations in Section <ref>, and

* a correctness proof of our algorithm in Section <ref>.

After a review of the work done, we conclude in Section <ref>.

§.§ Notation

Statements starting with a turnstile ($\vdash$) are HOL4 theorems,
automatically pretty-printed to  from the relevant theory in the HOL4 development.
Generally, our notation allows an appealing combination of quantifiers ($\forall, \exists, \exists{!}$),
logical connectives ( for “and”,  for “or”,  for “not”, also  for “implies” and  for “if and only if”),
set theory ($\in$ for “element of”, $\times$ for Cartesian product, and comprehensions such as x  x  6),
and functional programming ( for abstraction, and juxtaposition for application).
Repeated application of a function $f$ is indicated by exponents,
, f (f (f x)) = f.

For a function $f$ from set $S$ to set $T$,
we write f  S  T to mean a bijection.

The empty set is denoted by ,
and a finite set, denoted by finite S,

has cardinality S.

The set of natural numbers is denoted by $\mathbb{N}$, counting from $0$,
and count n  x  x  n,
where  means `equality by definition'.
For a natural number $n \in \mathbb{N}$,
square n means it is a square: k. n = k,

prime n means it is a prime, and
even n or odd n denotes its parity.

The integer quotient and remainder of $m$ divided by $n$ are written as m div n and m mod n, respectively.

We write n  m when $n$ divides $m$,
which is equivalent to m  0  n when n  0.

These are basic notations. Others will be introduced as they first appear.

*HOL4 Sources
Proof scripts are located in a repository at <https://bitbucket.org/jhlchan/project/src/master/fermat/twosq/>.

The scripts are compiled using HOL4, version .
In this paper, each theorem has windmill197,
which is hyperlinked to the appropriate line of the corresponding proof script in repository.
Proof scripts of all theorems are located in a repository (omitted for anonymous review).
The scripts are compiled using HOL4, version .
In this paper, each theorem has windmill197,
which is hyperlinked to the appropriate line of the corresponding proof script in repository.
(For anonymous review, this feature has been removed. See Appendix <ref> for a cross-reference of theorems in this paper and the proof scripts in supplementary material.)

§ SUM OF TWO SQUARES

The only even prime is 2 = 1  1, a sum of two squares.
An odd prime, upon division by 4, leaves a remainder of either $1$ or $3$.
Only an odd prime of the first type can be expressed as a sum of two squares,
as supported by numerical evidence from Table <ref>.

Examples of odd primes that can be expressed as a sum of two squares.

This table provides examples of odd primes that can be expressed as a sum of two squares.

\[
\begin{array}{r@{\quad\ee\quad}l@{\quad\ee\quad}l}
5 & 4(1) + 1  & 1^{2} + 2^{2}\\
13 & 4(3) + 1  & 3^{2} + 2^{2}\\
17 & 4(4) + 1  & 1^{2} + 4^{2}\\
29 & 4(7) + 1  & 5^{2} + 2^{2}\\
37 & 4(9) + 1  & 1^{2} + 6^{2}\\
41 & 4(10) + 1  & 5^{2} + 4^{2}\\
53 & 4(13) + 1  & 7^{2} + 2^{2}\\
61 & 4(15) + 1  & 5^{2} + 6^{2}\\
\end{array}
\]

Pierre de Fermat, in a letter to Marin Mersenne on Christmas day 1640,
claimed that he had an “irrefutable” proof of this:

A prime $n$ can be expressed uniquely as a sum of odd and even squares
if and only if n = 4k  1 for some $k$.

prime n

(n  1  4

(u,v). odd u  even v  n = u  v)

This paper concentrates on formalising an elementary proof of this result by Roger Heath-Brown, later simplified by Don Zagier.
As shown in his one-sentence proof in Section <ref>,
the idea is this: look at the representations of $n$ not by squares, but in another form.

Consider the following set $S_{n}$ of triples $(x,y,z)$:

\begin{equation}
\label{eqn:mills-set}
S_{n} = \{(x,y,z) \in \mathbb{N}\times{\mathbb{N}}\times{\mathbb{N}}
\mid \HOLinline{\HOLFreeVar{n}\;\HOLSymConst{=}\;\HOLFreeVar{x}\HOLSymConst{\ensuremath{\sp{2}}}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{4}\HOLSymConst{\ensuremath{}}\HOLFreeVar{y}\HOLSymConst{\ensuremath{}}\HOLFreeVar{z}} \}.
\end{equation}

For a prime $n$ of the form $4k + 1$, we have $(1,1,k) \in S_{n}$.
Thus the set $S_{n}$ is non-empty, and there are only finitely many triples in $S_{n}$.

A triple with y = z will give n = x  4y, , a sum of two squares.
If we can show that $S_{n}$ has only one such triple,
we have a proof of Fermat's Theorem <ref>, with both existence and uniqueness.

Meanwhile, some general theories will be developed as an exercise in formal proofs, so that they can be applied to similar problems.
In addition, we extend the theories to establish not only an algorithm, but also a proof of its correctness, to compute the two unique squares for primes of the form 4k  1.

at (1.5,2.2) $x$;
at (4.5,4.5) $y$;
at (8.5,3.5) $z$;

at (15.5,2.2) $x$;
at (16.0,5.5) $y$;
at (18.5,4.2) $z$;

at (29.5,3.4) $x$;
at (28.0,6.6) $y$;
at (29.5,5.0) $z$;

Typical windmills, where windmill x y z = x  4yz. The rightmost one has y = z.

This figure shows typical windmills. The central square is x by x, the four rectangles arranged clockwise around the square are all y by z. The rightmost one has y and z equal.

§.§ Windmills

The following expression will be our main focus:

A windmill consists of a central square with four identical rectangular arms.

windmill x y z  x  4yz

Some typical windmills are shown in Figure <ref>.
The first term x is given by a central square of side $x$,
and the second term 4yz is given by four arms, each a rectangle of width $y$ and height $z$, arranged clockwise around the square.

Therefore each triple in the set $S_{n}$ of Equation (<ref>) can be represented by a windmill, that is,
each triple $(x,y,z)$ satisfies n = windmill x y z.
Given a prime n = 4k  1,
we shall look for a windmill with four square arms (the one on the far right in Figure <ref>),
, y = z, so that n = x  4yy x  (2y).
First, we collect all triples $(x,y,z)$ which are solutions of n = x  4yz:

The mills of a number is its set of windmills.

mills n  (x,y,z)  n = windmill x y z

This is the formal definition of the set $S_{n}$ of Equation (<ref>).

The conditions for a proper windmill, with all lengths nonzero, are:

square n  n  0  4

x y z. (x,y,z)  mills n  x  0  y  0  z  0

\begin{equation}
\label{eqn:mills-triple-nonzero}
\end{equation}

When $n$ is a square, $\HOLinline{\HOLFreeVar{n}\;\HOLSymConst{=}\;\HOLFreeVar{x}\HOLSymConst{\ensuremath{\sp{2}}}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{4}\HOLSymConst{\ensuremath{}}\HOLFreeVar{y}}(0)$ for any value of $y$.
This would make mills n infinite. Otherwise:

The number of windmills for a number $n$ is finite if and only if $n$ is not a square.

finite (mills n)  square n

Given an odd $n$ that is not a square, we can determine all its windmill triples $(x,y,z)$ by noting that, since 4yz is even, $x$ must be odd, and $y$ and $z$ form the product yz = (n  x) div 4.
In Table <ref> this is worked out for n = 29, using successive odd $x$ and factors for the product $yz$.
The corresponding windmills are shown in Figure <ref>.

Determine all the windmill triples of n = 29, by odd $x$ and factors of $yz$.

This table shows how to detemine all the windmill triples for n = 29, using odd x and factors of the product yz.

odd $x$     n  x     4yz     triple $(x,y,z)$     comment

$1$     $29 - 1^{2}$     $28 \ee 4(7)$     $(1,1,7), (1,7,1)$     factors of $7$ are $1, 7$.

$3$     $29 - 3^{2}$     $20 \ee 4(5)$     $(3,1,5), (3,5,1)$     factors of $5$ are $1, 5$.

$5$     $29 - 5^{2}$      $4 \ee 4(1)$     $(5,1,1)$              factor of $1$ is $1$.

[step=1, color=white!60!black] (0,0) grid (65,17);

at (4,1)  $(1,1,7)$;
at (22,1) $(3,1,5)$;
at (36,1) $(5,1,1)$;
at (45,1) $(3,5,1)$;
at (55,1) $(1,7,1)$;

All the windmills of n = 29, determined from Table <ref>.

This figure shows all the windmills of n = 29, determined from the previous table.

When a number $n$ has the form 4k  1,

\[
n \ee 1^{2} + 4(1)k \ee \HOLinline{\HOLConst{windmill}\;\HOLNumLit{1}\;\HOLNumLit{1}\;\HOLFreeVar{k}},
\]

showing that its mills n  :

n  1  4  (1,1,n div 4)  mills n

Moreover, when this form corresponds to a prime, this is the only triple $(x,y,z)$ with x = y:

For a prime of the form 4k  1, the only windmill with the first and second parameters equal is windmill 1 1 k.

prime n  n  1  4

x z. n = windmill x x z  x = 1  z = n div 4

Note that k = n div 4 for prime n = 4k  1.
Consider (x,y,z)  mills n with x = y.
This implies,

\[
\HOLinline{\HOLFreeVar{n}\;\HOLSymConst{=}\;\HOLConst{windmill}\;\HOLFreeVar{x}\;\HOLFreeVar{x}\;\HOLFreeVar{z}} \ee \HOLinline{\HOLFreeVar{x}\HOLSymConst{\ensuremath{\sp{2}}}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{4}\HOLSymConst{\ensuremath{}}\HOLFreeVar{x}\HOLSymConst{\ensuremath{}}\HOLFreeVar{z}\;\HOLSymConst{=}\;\HOLFreeVar{x}\HOLSymConst{\ensuremath{}}(\HOLFreeVar{x}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{4}\HOLSymConst{\ensuremath{}}\HOLFreeVar{z})}.
\]

Therefore x  n.
As prime $n$ is not a square, x  n.
Hence x = 1, so y = 1, and z = k.

§.§ Involution

We are going to study involutions on mills n, the set of windmills for $n$.
A function $f$ is an involution on a set $S$, denoted by f involute S, when it is its own inverse:

f involute S  x. x  S  f x  S  f (f x) = x

That is, $f$ is a bijection f  S  S, pairing up $x$ and f x, both in $S$. When x = f x, the element $x$ is fixed by the involution $f$.

We define the following sets:

The pairs and fixes of an involution $f$ on a set $S$.

pairs f S  x  x  S  f x  x

fixes f S  x  x  S  f x = x

Clearly they are disjoint.
The subset $\HOLinline{\HOLConst{pairs}\;\HOLFreeVar{f}\;\HOLFreeVar{S}}$ consists of distinct involute pairs, so its cardinality is even:

finite S  f involute S  even pairs f S

So both S and fixes f S have the same parity. This leads to:

If two involutions act on the same finite set $S$, their fixes have the same parity.

finite S  f involute S  g involute S

(odd fixes f S  odd fixes g S)

We shall meet the two involutions on mills n, a set which is finite for non-square $n$ (by Theorem <ref>).

§ WINDMILL INVOLUTIONS

Zagier's one-sentence proof is the interplay of two involutions on the set of windmills (mills n) for a prime n = 4k  1.

§.§ Flip Map

The first involution just swaps the $y$ and $z$ in the triple $(x,y,z)$:

The flip map for a triple.

flip (x,y,z)  (x,z,y)

The set S = mills n of windmill triples of a number $n$ can be partitioned by $y, z$ into:

\[
\begin{array}{c}
S_{y < z} \ee \{(x,y,z) \in S \mid y < z\}\\
S_{y \ee z} \ee \{(x,y,z) \in S \mid y \ee z\}\\
S_{y > z} \ee \{(x,y,z) \in S, \mid y > z\}\\
\end{array}
\]

An example for n = 29 is shown in Figure <ref>.
Clearly there is a bijection: $\HOLConst{flip}\colon S_{y < z} \leftrightarrow S_{y > z}$,

and $S_{y \ee z} \ee \HOLinline{\HOLConst{fixes}\;\HOLConst{flip}\;\HOLFreeVar{S}}$.
Thus the inverse of flip is itself:

flip (flip (x,y,z)) = (x,y,z)

showing that:

The flip map is an involution on the set of windmills.

flip involute mills n

[step=1, color=white!60!black] (0,0) grid (65,17);

at (4,1)  $(1,1,7)$;
at (22,1) $(3,1,5)$;
at (36,1) $(5,1,1)$;
at (45,1) $(3,5,1)$;
at (55,1) $(1,7,1)$;

(a) at (0.5,0.5);
(b) at (31,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;
(a) at (32.5,0.5);
(b) at (40,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;
(a) at (41.5,0.5);
(b) at (64.5,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;

(a) at (1.0,1.1);
(b) at (6.0,1.1);
[<->] (a) to [bend left,looseness=0.8] node[midway,below] flip (b);
(a) at (2.5,1.2);
(b) at (4.6,1.2);
[<->] (a) to [bend left,looseness=0.8] node[midway,below] flip (b);
(a) at (3.6,0.5);
(b) at (3.6,0.51);
[<->] (a) to [out=-30, in=-150, looseness=200]
node[midway,above] flip (b);

Partition of windmills of n = 29 for flip: those with $y < z, y \ee z$, and $y > z$. Note left and right pairing.

This figure shows a partition of the windmills of n = 29 for the flip map: those with y < z, y = z, and y > z. Note the left and right pairing between y < z and y > z.

§.§ Zagier Map

The other involution is the one devised by Don Zagier, as shown in Equation (<ref>):

The Zagier map for a triple.

zagier (x,y,z)

if x  y  z then (x  2z,z,y  z  x)

else if x  2y then (2y  x,y,x  z  y)

else (x  2y,x  z  y,y)

Algebraically, this is indeed an involution, as HOL4 can verify without a blink:

x  0  z  0  zagier (zagier (x,y,z)) = (x,y,z)

\begin{equation}
\label{eqn:zagier-involute}
\end{equation}

That HOL4 can verify this directly from definition is a showcase of its excellent algebraic simplifier, especially for natural numbers.
However, we would like to see the magic behind, in terms of the geometry of windmills.
Note that this definition differs slightly from Equation (<ref>) since the else-parts include boundary cases.
They actually correspond to improper windmills, and they are irrelevant for the values of $n$ satisfying Equation (<ref>).

§.§ Mind of a Windmill

The main purpose of introducing windmills is to read their minds.

at (1.5,2.2) $x$;
at (5.0,4.5) $y$;
at (8.5,3.5) $z$;

(14,2.8) – (17,2.8) node[midway,below,yshift=-7pt] $x$;
(13,4.2) – (18,4.2) node[midway,above,yshift=7pt] $x'$;
at (29.5,3.4) $x'$;
at (27.5,8.6) $y'$;
at (28.5,6.2) $z'$;

A typical windmill x y z = x  4yz, with a mind (in dashes) and transforms to another windmill.

This figure shows a typical windmill, with a mind in dashes, on the left. The figure also illustrates how the left windmill transforms to another windmill on the right, with the same mind.

Referring to Figure <ref>, a windmill has a mind (marked in dashes at middle), which is the maximum central square, with side $x'$, that can be fitted with the four arms.
When x  , the original square x can grow to the mind , forming another windmill but keeping the overall shape (on the right).
Conversely, going from right to left, we can use the mind as a reference to shrink the square term from  to x by trimming four sides,
thereby restoring the arms to original.
Transforming a windmill's square term through the mind is the geometric interpretation of Equation (<ref>).

The five cases of Zagier map, transforming a triple $(x,y,z)$ to $(x',y',z')$.

This figure shows the five cases of Zagier map, transforming a triple (x,y,z) to (x',y',z').

Case     Type     condition     Mind     Picture     $x'$     $y'$     $z'$     condition     Type

$1$      2*$x < y$     $x < y - z$       $x + 2z$     Figure <ref> (a)
$x + 2z$     $z$     $y - x - z$     $2y' < x'$      2*$y' < x'$

$2$          $y - z < x$     $2y - x$     Figure <ref> (b)
$2y - x$     $y$     $x + z - y$     $x' < 2y'$

$3$     $x = y$                 $x$     Figure <ref> (c)
$x$          $y$     $z$                  $x' = y'$

$4$      2*$y < x$     $x < 2y$      $x$     Figure <ref> (d)
$2y - x$     $y$     $x + z - y$     $y' - z' < x'$      2*$x' < y'$

$5$         $2y < x$        $x$     Figure <ref> (e)
$x - 2y$     $x + z - y$     $y$     $x' < y' - z'$

at (1.5,2.2) $x$;
at (4.0,5.6) $y$;
at (8.5,4.0) $z$;
[->,thick] (10,3) – (12,3) node[midway,above] zagier;
at (18.5,4.0) $x'$;
at (26.0,4.0) $y'$;
at (23.5,1.5) $z'$;

[color=red,ultra thick] (5,3) – (8,3);
[color=red,ultra thick] (22,3) – (25,3);
(5,3) – (8,3);
(22,3) – (25,3);

at (0.8,10) (a) Case $1$: x  y  z.;
at (33,3) $\large{
\begin{array}{r@{\;=\;}l}
x' & x + 2z\\
y' & z\\
z' & y - x - z\\
\end{array}

at (1.5,2.2) $x$;
at (3.5,7.5) $y$;
at (6.5,4.5) $z$;
[->,thick] (10,3) – (12,3) node[midway,above] zagier;
at (19.5,5.0) $x'$;
at (21.0,7.6) $y'$;
at (17.0,6.8) $z'$;

[color=red,ultra thick] (0,6) – (0,7);
[color=red,ultra thick] (18,6) – (18,7);
(0,6) – (0,7);
(18,6) – (18,7);

at (6,10) ;
at (6,9) (b) Case $2$: y  z  x and x  y, so x  2y.;
at (33,3) $\large{
\begin{array}{r@{\;=\;}l}
x' & 2y - x\\
y' & y\\
z' & x + z - y\\
\end{array}
at (0,-4) ;

at (2.2,3.5) $x$;
at (2.2,6.5) $y$;
at (4.5,5.2) $z$;
[->,thick] (10,3) – (12,3) node[midway,above] zagier;
at (18.0,3.5) $x'$;
at (18.0,6.5) $y'$;
at (14.8,5.2) $z'$;

[color=red,ultra thick] (0,4) – (0,6);
[color=red,ultra thick] (16,4) – (16,6);
(0,4) – (0,6);
(16,4) – (16,6);

at (7,9) ;
at (6.5,8) (c) Case $3$: y = x, so y  z  x and x  2y.;
at (33,3) $\large{
\begin{array}{r@{\;=\;}l}
x' & 2y - x = x\\
y' & y\\
z' & x + z - y = z\\
\end{array}

at (3.0,5.5) $x$;
at (2.5,8.5) $y$;
at (4.5,6.8) $z$;
[->,thick] (10,3) – (12,3) node[midway,above] zagier;

at (19.0,3.6) $x'$;
at (18.0,8.6) $y'$;
at (14.5,6.0) $z'$;

[color=red,ultra thick] (0,4) – (0,8);
[color=red,ultra thick] (16,4) – (16,8);
(0,4) – (0,8);
(16,4) – (16,8);

at (3,11) ;
at (3,10) (d) Case $4$: y  x, but x  2y.;
at (34,3) $\large{
\begin{array}{r@{\;=\;}l}
x' & 2y - x\\
y' & y\\
z' & x + z - y\\
\end{array}

at (3.5,6.2) $x$;
at (1.5,9.5) $y$;
at (3.5,8.0) $z$;
[->,thick] (10,3) – (12,3) node[midway,above] zagier;
at (19.5,3.5) $x'$;
at (22.5,7.8) $y'$;
at (25.6,5.8) $z'$;

[color=yellow,ultra thick] (0,9) – (3,9);
[color=yellow,ultra thick] (4,-2) – (7,-2);
[color=yellow,ultra thick] (16,3) – (19,3);
[color=yellow,ultra thick] (20,4) – (23,4);
(0,9) – (3,9);
(4,-2) – (7,-2);
(16,3) – (19,3);
(20,4) – (23,4);

at (3,12) ;
at (3,11) (e) Case $5$: 2y  x, so y  x.;
at (34,3) $\large{
\begin{array}{r@{\;=\;}l}
x' & x - 2y\\
y' & x + z - y\\
z' & y\\
\end{array}

All five cases of the Zagier map, from $(x,y,z)$ to $(x',y',z')$ through the mind of a windmill.

This figure shows all the 5 cases of the Zagier map, transform through the mind of a windmill.

The Zagier map transforms $(x,y,z)$ to $(x',y',z')$ via the mind of the windmill, keeping its overall shape.
There are three types, depending on whether $x < y$, $x = y$, or $y < x$.
Both the first and last types are divided into two cases, as the geometry for the mind is different.
Altogether there are five cases, as analysed in Table <ref>, and illustrated in Figure <ref>.[Dubach and Muehlboeck <cit.> also identified five types for windmills.]

Although five cases of Zagier map have been identified, note that the transformation rule:

\[
(x',y',z') \ee (2y - x, y, x + z - y)
\]

happens to be the same for case $2$ and case $4$.
The same rule actually applies to case $3$, which has x = y.
Thus the Zagier map can be succinctly expressed as in Definition <ref> with only three branches.

Moreover, we can define the mind of a windmill triple as (see Table <ref>):

mind (x,y,z)

if x  y  z then x  2z

else if x  y then 2y  x

else x

and verify that the mind is an invariant under the Zagier map for any triple:

mind (zagier (x,y,z)) = mind (x,y,z)

Referring again to Table <ref>,
the windmills in S = mills n can be partitioned into three triple types:

\[
\begin{array}{c@{\quad}l}
S_{x < y} \ee \{(x,y,z) \in S \mid x < y\}
& \text{covering cases $1$ and $2$}\\
S_{x \ee y} \ee \{(x,y,z) \in S \mid x \ee y\}
& \text{covering case $3$}\\
S_{x > y} \ee \{(x,y,z) \in S \mid x > y\}
& \text{covering cases $4$ and $5$}\\
\end{array}
\]

Such a partition for n = 29 is shown in Figure <ref>.
Table <ref> also shows that, for triples with proper windmills:

* a triple of case $1$ maps to case $5$ and vice versa,

* a triple of case $2$ maps to case $4$ and vice versa, and

* a triple of case $3$ maps to itself.

Therefore the Zagier map is its own inverse for proper triples.
Combining Equation (<ref>) and Equation (<ref>)
for the windmills of a prime, we have:

The Zagier map is an involution on mills n for a prime $n$.

prime n  zagier involute mills n

[step=1, color=white!60!black] (0,0) grid (65,17);

at (4,1)  $(1,1,7)$;
at (22,1) $(1,7,1)$;
at (36,1) $(3,5,1)$;
at (46,1) $(3,1,5)$;
at (60,1) $(5,1,1)$;

(a) at (0.5,0.5);
(b) at (16.5,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;
(a) at (18,0.5);
(b) at (40,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;
(a) at (41.5,0.5);
(b) at (64.5,16.5);
[draw,thick,color=purple,fit= (a) (b),rounded corners=.55cm,inner sep=2pt] ;

[ultra thick]881

(a) at (2.5,1.2);
(b) at (4.6,1.2);
[<->] (a) to [bend left] node[midway,below] zagier (b);
(a) at (3.4,0.5);
(b) at (6.0,0.5);
[<->] (a) to [bend right, looseness=0.8] node[midway,above] zagier (b);
(a) at (1.2,0.6);
(b) at (1.2,0.61);
[<->] (a) to [out=-30, in=-150, looseness=250]
node[midway,above] zagier (b);

Partition of windmills of n = 29 for zagier: those with $x \ee y, x < y$, and $x > y$. Note pairing by minds.

This figure shows a partition of windmills of n = 29 for the Zagier map, those with x = y, x < y, and x > y. Note the pairing of minds between those x < y and x > y.

§ TWO SQUARES THEOREM

Now we have enough tools to formalise Fermat's two squares theorem.

§.§ Existence of Two Squares

For the Zagier map,
it is straightforward to verify, as indicated in Table <ref>, that only a triple of case $3$ can map to itself:

x  0  (zagier (x,y,z) = (x,y,z)  x = y)

Hence $S_{x \ee y} \ee \HOLinline{\HOLConst{fixes}\;\HOLConst{zagier}\;(\HOLConst{mills}\;\HOLFreeVar{n})}$.
Applying Theorem <ref> which characterises such triples, for certain primes $S_{x \ee y}$ is a singleton:

A prime of the form 4k  1 has only $(1,1,k)$ fixed by the Zagier map.

prime n  n  1  4

fixes zagier (mills n) = (1,1,n div 4)

The fixed points of two involutions play crucial roles in the existence of two squares for Theorem <ref>:

A prime of the form 4k  1 is a sum of two squares of different parity.

prime n  n  1  4

(u,v). odd u  even v  n = u  v

A prime is not a square, so mills n is finite by Theorem <ref>.
, and both Zagier and flip maps are involutions on mills n,
by Theorem <ref> and Theorem <ref>.
Note that Zagier map has a single fixed point by Theorem <ref>.
Thus fixes zagier (mills n) = 1,
so fixes flip (mills n) is odd by Theorem <ref>.
Hence fixes flip (mills n)  , containing a triple $(x,y,y)$.
Thus n = windmill x y y $\ee$ x  4y.
Take u = x, and v = 2y, then n = u  v.
Evidently v is even, and $u$ is odd since $n$ is odd.

Current formalisations of Zagier's proof (HOL Light <cit.>, NASA PVS <cit.> and Coq <cit.>), or its close relative Heath-Brown's proof (Mizar <cit.> and ProofPower <cit.>), stop at just showing the existence of two squares for the primes in Fermat's Theorem <ref>,
most likely because this already meets the Formalizing 100 Theorems challenge <cit.>.

See also related work in Section <ref>.

§.§ Uniqueness of Two Squares

The uniqueness of the two squares in Fermat's Theorem <ref>
is a consequence of the following property of a prime:

If a prime $n$ can be expressed as a sum of two squares, the expression is unique up to commutativity.

prime n  n = a  b  n = c  d

a; b = c; d

The proof is purely number-theoretic, which has also been formalised by Laurent Théry in Coq <cit.>.
Moreover, we have:

A number of the form 4k  3 cannot be expressed as a sum of two squares.

n  3  4  u v. n  u  v

This is an elementary result from possible remainders after division by 4:
while a number, such as $u$ or $v$, may have a remainer $0, 1, 2$, or $3$,
a square, such as $u^{2}$ or $v^{2}$, can only have a remainder $0$ or $1$.
Thus the sum of such remainders can never be $3$.

Now we can complete the proof of Fermat's two squares Theorem <ref>:

prime n

(n  1  4

(u,v). odd u  even v  n = u  v)

For the if part $(\Rightarrow)$,
existence is given by Theorem <ref>, and
uniqueness is provided by Theorem <ref>.
For the only-if part $(\Leftarrow)$,
an odd prime with n  1  4 cannot be a sum of two squares by Theorem <ref>.

§ TWO SQUARES ALGORITHM

To make Zagier's proof constructive, we need to compute that single triple fixed by flip map.

Let $n$ be a prime of the form 4k  1.
By Theorem <ref>, the only Zagier fixed point is u = (1,1,k), meaning zagier u = u.
To change the triple $u$, applying flip is the obvious choice.
To keep changing the triple, zagier should be applied.
Thus by applying the composition zagier  flip repeatedly from the known Zagier fixed point, there is hope that the chain will lead to the only flip fixed point.
Figure <ref> shows that this is indeed the case for n = 29.

[step=1, color=white!60!black] (0,0) grid (65,17);

at (4,1)  $(1,1,7)$;
at (22,1) $(1,7,1)$;
at (38,1) $(3,1,5)$;
at (51,1) $(3,5,1)$;
at (60,1) $(5,1,1)$;

[ultra thick]881

(a) at (1.0,1.1);
(b) at (2.2,1.1);
[->] (a) to [bend left] node[midway,above] flip (b);
(a) at (2.8,0.6);
(b) at (3.8,0.6);
[->] (a) to [bend right] node[midway,below] zagier (b);
(a) at (4.2,1.2);
(b) at (5.2,1.2);
[->] (a) to [bend left] node[midway,above] flip (b);
(a) at (5.2,0.5);
(b) at (6.0,0.5);
[->] (a) to [bend right] node[midway,below] zagier (b);

The iteration chain of n = 29 by the composition zagier  flip, from Zagier fix to flip fix.

This figure show the iteration chain of n = 29, by the composition of first flip then Zagier. The chain starts from Zagier fix, ends in flip fix.

In terms of windmills, the flip map keeps the central square, but flips the arms of rectangles from $y$-by-$z$ to $z$-by-$y$. This generally changes the mind of the windmill. The Zagier map keeps the mind, but changes the central square.
Similar to the mind being an invariant of the Zagier map,
the absolute difference $\left|y - z\right|$ is an invariant of the flip map.
If the Zagier map can reduce this difference, the successive iterations of zagier  flip will be able to locate the flip fixed point.

§.§ Flip Fix Search

To find the fixed point of the flip map, we can experiment with this pseudo-code:

* Input: a number n = 4k  1.

* Output: a triple fixed by the flip map.

* Method:

* start with u = (1,1,k), the Zagier fix.

* while ($u$ is not a flip fix) :

*     $u \leftarrow \HOLinline{(\HOLConst{zagier}\;\HOLSymConst{\HOLTokenCompose}\;\HOLConst{flip})\;\HOLFreeVar{u}}$

* end while.

In an HOL4 interactive session, this pseudo-code can be implemented directly as:[This pseudo-code can be implemented directly in any programming language that supports while-loops and tuples.]

Computing the flip fixed point of n = 4k  1 using a WHILE loop.

two_sq n

WHILE (()  found) (zagier  flip) (1,1,n div 4),

found (x,y,z)  y = z

This simple while-loop may or may not terminate. We shall take up this issue in Section <ref>.

For primes of the form 4k  1, it terminates and seems to work.
To prove its correctness, we shall develop a theory of permutation iteration, then apply the theory to this algorithm.

§ PERMUTATION ORBITS

In general, the composition of two involutions is no longer an involution, but just a permutation.
Let $\varphi \colon S \rightarrow S$ be a permutation,
a bijection on the set $S$, denoted by φ permutes S.
For an element x  S
the iteration sequence $\varphi(x)$, φ, φ, , form its orbit.
The smallest positive index $n$ such that φ = x is called the period of $x$ under $\varphi$.
If such a positive index does not exist, the period is defined to be $0$.

In HOL4, the definition makes use of OLEAST, the optional LEAST operator:

The period of function iteration of an element is the least nonzero index for the element iterate to wrap around, otherwise zero.

period φ x

case OLEAST k. 0  k  φ = x of

none  0

some k  k

When the set $S$ is finite, the iterates cannot be always distinct.
Thus the permutation orbit of any $x \in S$ is finite, with a nonzero period,
denoted by p = period φ x:

finite S  φ permutes S  x  S

p. 0  p  p = period φ x

and by definition the period is minimal, which means that there is no wrap around for element iterates when the index is less than the period:

0  j  j  period φ x  φ  x

This implies a criterion for an exponent index to be divisible by period:

For a nonzero period $p$ of $x$, $x$ is fixed by the $k$-th iterate of $\varphi$ if and only if $k$ is a multiple of period $p$.

0  p  p = period φ x

(φ = x  k  0  p)

Moreover, the period is the same for all iterates in the same orbit:

finite S  φ permutes S  x  S  y = φ

period φ y = period φ x

§.§ Involution Composition

When the permutation φ = f  g, a composition of two involutions $f$ and g, we shall investigate whether their fixed points are connected by a chain of composition iterations.
Note the following pattern of function application:

\begin{equation*}
\label{eqn:function-assoc}
\begin{split}
f\ \circ\ (\HOLinline{\HOLFreeVar{g}}\ \circ\ f)\ \circ\ (\HOLinline{\HOLFreeVar{g}}\ \circ\ f)\ \circ\ (\HOLinline{\HOLFreeVar{g}}\ \circ\ f)\\
\ee (f\ \circ\ \HOLinline{\HOLFreeVar{g}})\ \circ\ (f\ \circ\ \HOLinline{\HOLFreeVar{g}})\ \circ\ (f\ \circ\ \HOLinline{\HOLFreeVar{g}})\ \circ\ f
\end{split}
\end{equation*}

by associativity.
Also, $(\HOLinline{\HOLFreeVar{f}\;\HOLSymConst{\HOLTokenCompose}\;\HOLFreeVar{g}})^{-1} \ee \HOLinline{\HOLFreeVar{g}}^{-1}\ \circ\ f^{-1} \ee \HOLinline{\HOLFreeVar{g}\;\HOLSymConst{\HOLTokenCompose}\;\HOLFreeVar{f}}$ for involutions,
so inverse is just reversal of application order in this case.
Let p = period (f  g) x for x  S.
With these notations, we can establish some basic results:

When $f$ fixes $x$, the period for $x$ is $1$ if and only if g also fixes $x$.

f involute S  g involute S  x  fixes f S

p = period (f  g) x

(p = 1  x  fixes g S)

Pick an element $x$ in the set $S$.
For involutions, an iterate of (f  g) can be equal to another iterate of (g  f):

The $i$-th iterate of (f  g) equals the $j$-th iterate of (g  f) if and only if (i  j) is a multiple of period $p$.

finite S  f involute S  g involute S  x  S

p = period (f  g) x

((f  g) = (g  f)  i  j  0  p)

When $f$ fixes point $x$,
the iterates (f  g) and (f  g) are related when the sum (i  j) is special:

When $f$ fixes $x$, the $i$-th and $j$-th iterate of (f  g) differ by one $f$ application if and only if (i  j) is a multiple of period $p$.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x

((f  g) = f ((f  g))  i  j  0  p)

There is a related result, with a similar proof:

When $f$ fixes $x$, the $i$-th and $j$-th iterate of (f  g) differ by one g application if and only if (i  j  1) is a multiple of period $p$.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x

((f  g) = g ((f  g))

i  j  1  0  p)

These theorems are useful in the study of iteration orbits starting from fixed points.

§.§ Period Parity

Given a finite set $S$, and an element $x \in S$, the iterates (f  g) form an orbit, with length equal to the period p = period (f  g) x.
Figure <ref> shows two orbits, one with an even period, the other with an odd period.

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$x$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ x$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ x$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ x$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ x$] (e) at (2,1) ;
[ele,label=above:$\varphi^{5}\ x$] (f) at (1,1) ;
[ele,label=above:$\varphi^{6}\ x\qquad$] (g) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);

[ele,draw,fill=white] (ab) at ($(a)!0.5!(b)$) ;
[ele,draw,fill=white] (bc) at ($(b)!0.5!(c)$) ;
[ele,draw,fill=white] (cd) at ($(c)!0.5!(d)$) ;
[ele,draw,fill=white] (de) at ($(d)!0.5!(e)$) ;
[ele,draw,fill=white] (ef) at ($(e)!0.5!(f)$) ;
[ele,draw,fill=white] (fg) at ($(f)!0.5!(g)$) ;

[thick,color=purple] (a) to [bend right] node[midway,below] g (ab);
[thick,color=blue] (ab) to [bend right] node[midway,below] $f$ (b);
[thick,color=purple] (b) to [bend right] node[midway,below] g (bc);
[thick,color=blue] (bc) to [bend right] node[midway,below] $f$ (c);
[thick,color=purple] (c) to [bend right] node[midway,below] g (cd);
[thick,color=blue] (cd) to [bend right] node[midway,below] $f$ (d);
[thick,color=purple] (d) to [bend right] node[midway,above] g (de);
[thick,color=blue] (de) to [bend right] node[midway,above] $f$ (e);
[thick,color=purple] (e) to [bend right] node[midway,above] g (ef);
[thick,color=blue] (ef) to [bend right] node[midway,above] $f$ (f);
[thick,color=purple] (f) to [bend right] node[midway,above] g (fg);
[thick,color=blue] (fg) to [bend right] node[midway,above] $f$ (g);

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$x$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ x$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ x$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ x$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ x$] (e) at (2.5,1) ;
[ele,label=above:$\varphi^{5}\ x$] (f) at (1.5,1) ;
[ele,label=above:$\varphi^{6}\ x$] (g) at (0.5,1) ;
[ele,label=above:$\varphi^{7}\ x\qquad$] (h) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);
[->,dashed,shorten <=2pt,shorten >=2] (g) – (h);

[ele,draw,fill=white] (ab) at ($(a)!0.5!(b)$) ;
[ele,draw,fill=white] (bc) at ($(b)!0.5!(c)$) ;
[ele,draw,fill=white] (cd) at ($(c)!0.5!(d)$) ;
[ele,draw,fill=white] (de) at ($(d)!0.5!(e)$) ;
[ele,draw,fill=white] (ef) at ($(e)!0.5!(f)$) ;
[ele,draw,fill=white] (fg) at ($(f)!0.5!(g)$) ;
[ele,draw,fill=white] (gh) at ($(g)!0.5!(h)$) ;

[thick,color=purple] (a) to [bend right] node[midway,below] g (ab);
[thick,color=blue] (ab) to [bend right] node[midway,below] $f$ (b);
[thick,color=purple] (b) to [bend right] node[midway,below] g (bc);
[thick,color=blue] (bc) to [bend right] node[midway,below] $f$ (c);
[thick,color=purple] (c) to [bend right] node[midway,below] g (cd);
[thick,color=blue] (cd) to [bend right] node[midway,below] $f$ (d);
[thick,color=purple] (d) to [bend right] node[midway,above] g (de);
[thick,color=blue] (de) to [bend right] node[midway,above] $f$ (e);
[thick,color=purple] (e) to [bend right] node[midway,above] g (ef);
[thick,color=blue] (ef) to [bend right] node[midway,above] $f$ (f);
[thick,color=purple] (f) to [bend right] node[midway,above] g (fg);
[thick,color=blue] (fg) to [bend right] node[midway,above] $f$ (g);
[thick,color=purple] (g) to [bend right] node[midway,above] g (gh);
[thick,color=blue] (gh) to [bend right] node[midway,above] $f$ (h);

Orbits of φ = f  g for point $x$. Left one has even period $6$, right one has odd period $7$.

This figure shows the orbits for point x of the composition: first g than f. Left one has even period 6, right one has odd period 7.

In the figure, black dots indicate iterates of φ = f  g, in dashes, and white dots indicate the intermediates, with g first, then $f$, through the arcs.
Since $f$ and g are involutions, the arcs can go both ways: forward or backward.

Let $\alpha$ denote a fixed point of $f$, and $\beta$ denote a fixed point of g,
, f α = α, and g β = β.

We shall look at how these fixed points are related, which is crucial in the correctness proof of our algorithm (see Definition <ref>).

§.§ Fixed Point Period Even

Consider an orbit with even period starting with $\alpha$, a fixed point of $f$.
Figure <ref> shows one on the left, and its real picture on the right.

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$\alpha$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ \alpha$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ \alpha$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ \alpha$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ \alpha$] (e) at (2,1) ;
[ele,label=above:$\varphi^{5}\ \alpha$] (f) at (1,1) ;
[ele,label=above:$\varphi^{6}\ \alpha\qquad$] (g) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);

[ele,draw,fill=white] (ab) at ($(a)!0.5!(b)$) ;
[ele,draw,fill=white] (bc) at ($(b)!0.5!(c)$) ;
[ele,draw,fill=white] (cd) at ($(c)!0.5!(d)$) ;
[ele,draw,fill=white] (de) at ($(d)!0.5!(e)$) ;
[ele,draw,fill=white] (ef) at ($(e)!0.5!(f)$) ;

[thick,color=blue] (a) to [out=130, in=-130,looseness=50]
node[midway,left] $f$ (a);
[thick,color=purple] (a) to [bend right] node[midway,below] g (ab);
[thick,color=blue] (ab) to [bend right] node[midway,below] $f$ (b);
[thick,color=purple] (b) to [bend right] node[midway,below] g (bc);
[thick,color=blue] (bc) to [bend right] node[midway,below] $f$ (c);
[thick,color=purple] (c) to [bend right] node[midway,below] g (cd);
[thick,color=blue] (cd) to [bend right] node[midway,below] $f$ (d);
[thick,color=purple] (d) to [bend right] node[midway,above] g (de);
[thick,color=blue] (de) to [bend right] node[midway,above] $f$ (e);
[thick,color=purple] (e) to [bend right] node[midway,above] g (ef);
[thick,color=blue] (ef) to [bend right] node[midway,above] $f$ (f);
[thick,color=purple] (f) to [bend right] node[midway,above] g (g);

[ultra thick, glow=green] (cd) – (d);

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$\alpha$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ \alpha$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ \alpha$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ \alpha$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ \alpha$] (e) at (2,1) ;
[ele,label=above:$\varphi^{5}\ \alpha$] (f) at (1,1) ;
[ele,label=above:$\varphi^{6}\ \alpha\qquad$] (g) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);

[thick,color=blue] (a) to [out=130, in=-130,looseness=50]
node[midway,left] $f$ (a);
[thick,color=blue] (f) to [bend right] node[midway,left] $f$ (b);
[thick,color=blue] (e) to [bend left] node[midway,right] $f$ (c);
[thick,color=blue] (d) to [out=50, in=-50,looseness=50]
node[midway,right] $f$ (d);

[thick,color=purple] (c) to [bend right] node[midway,below] g (d);
[thick,color=purple] (e) to [bend left] node[midway,above] g (b);
[thick,color=purple] (f) to [bend right] node[midway,above] g (g);

Orbit from an $f$ fixed point $\alpha$ with even period $6$. Identical points on the left (marked by two parallel lines) are merged on the right (move white dot to black dot). In particular, on the left the two vertices of the shaded line are the same, forming a fixed point of $f$ on the right.

This figure shows an orbit from an f fixed point alpha, with even period 6. Identical points on the left, marked by two parallel lines, are merged on the right, by moving white dot to black dot. In particular, on the left the two vertices of the shaded line are the same, forming a fixed point of f on the right.

This orbit is formed by taking the left diagram of Figure <ref>,
but identifying the black dot on $\alpha$ (the leftmost one) with its preceding white dot from $f$, since f α = α, giving the left $f$-loop.
This node $\alpha$ is now preceded by two g-arcs, one from a black dot and one from a white dot. However, g is an involution, which is injective, so the two dots are identical. The same reasoning shows that all the dots linked by double lines are identical, so that the orbit on the left can be simplified to the one on the right, taking only black dots.

Moreover, the rightmost black dot and a preceding white dot from $f$ must be the same, due to g-arcs from identical dots.
This means the half-period iterate, the rightmost black dot, is another fixed point of $f$, say $\alpha'$.
Note that $\alpha' \ne \alpha$, for otherwise the period will be affected.
This example motivates the following:

When $f$ fixes $x$, and (f  g) has an even period $p$ for $x$,
then $f$ also fixes (f  g), which is not $x$ itself.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x

y = (f  g)  even p

y  fixes f S  y  x

First we show that $f$ fixes $y$.
Let h = p div 2.
Since period $p$ is even, $p \ee \HOLinline{\HOLNumLit{2}\HOLSymConst{\ensuremath{}}\HOLFreeVar{h}\;\HOLSymConst{=}\;\HOLFreeVar{h}\;\HOLSymConst{\ensuremath{+}}\;\HOLFreeVar{h}}$.
This implies that h  h  0  p,
so y = f y by Theorem <ref>.
Since (f  g) is a permutation, y  S, so y  fixes f S.
Next we show that y  x.
Suppose y = x.
Since for finite $S$ the period p  0,
Theorem <ref> shows that
$p$ divides h = p div 2.
Hence p = 1, which is not even.

Therefore if a fixed point of $f$ has an even period under φ = f  g, it is not alone. This leads directly to:

If $f$ fixes only a single $x$, then f  g has an odd period for $x$.

finite S  f involute S  g involute S

fixes f S = x  p = period (f  g) x

odd p

§.§ Fixed Point Period Odd

Now consider an orbit with odd period starting with $\alpha$, a fixed point of $f$.
Figure <ref> shows one on the left, and its real picture on the right.
This orbit is formed by taking the right diagram of Figure <ref>,
but identifying the black dot on $\alpha$ (the leftmost one) with its preceding white dot from $f$, since f α = α, giving the left $f$-loop.
The same reasoning as the even period orbit of Section <ref> shows that all the dots linked by double lines are identical, so that the orbit on the left can be simplified to the one on the right, again taking only black dots.

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$\alpha$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ \alpha$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ \alpha$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ \alpha$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ \alpha$] (e) at (2.5,1) ;
[ele,label=above:$\varphi^{5}\ \alpha$] (f) at (1.5,1) ;
[ele,label=above:$\varphi^{6}\ \alpha$] (g) at (0.5,1) ;
[ele,label=above:$\varphi^{7}\ \alpha\qquad$] (h) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);
[->,dashed,shorten <=2pt,shorten >=2] (g) – (h);

[ele,draw,fill=white] (ab) at ($(a)!0.5!(b)$) ;
[ele,draw,fill=white] (bc) at ($(b)!0.5!(c)$) ;
[ele,draw,fill=white] (cd) at ($(c)!0.5!(d)$) ;
[ele,draw,fill=white] (de) at ($(d)!0.5!(e)$) ;
[ele,draw,fill=white] (ef) at ($(e)!0.5!(f)$) ;
[ele,draw,fill=white] (fg) at ($(f)!0.5!(g)$) ;

[thick,color=blue] (a) to [out=130, in=-130,looseness=50]
node[midway,left] $f$ (a);
[thick,color=purple] (a) to [bend right] node[midway,below] g (ab);
[thick,color=blue] (ab) to [bend right] node[midway,below] $f$ (b);
[thick,color=purple] (b) to [bend right] node[midway,below] g (bc);
[thick,color=blue] (bc) to [bend right] node[midway,below] $f$ (c);
[thick,color=purple] (c) to [bend right] node[midway,below] g (cd);
[thick,color=blue] (cd) to [bend right] node[midway,below] $f$ (d);
[thick,color=purple] (d) to [bend right] node[midway,above] g (de);
[thick,color=blue] (de) to [bend right] node[midway,above] $f$ (e);
[thick,color=purple] (e) to [bend right] node[midway,above] g (ef);
[thick,color=blue] (ef) to [bend right] node[midway,above] $f$ (f);
[thick,color=purple] (f) to [bend right] node[midway,above] g (fg);
[thick,color=blue] (fg) to [bend right] node[midway,above] $f$ (g);
[thick,color=purple] (g) to [bend right] node[midway,above] g (h);

[ultra thick, glow=green] (de) – (d);

ele/.style=fill=black,circle,minimum width=.8pt,inner sep=1pt,
every fit/.style=ellipse,draw,inner sep=5pt]
[ele,label=left:$\alpha$] (a) at (0,0.5) ;
[ele,label=below:$\varphi\ \alpha$] (b) at (1,0) ;
[ele,label=below:$\varphi^{2}\ \alpha$] (c) at (2,0) ;
[ele,label=right:$\varphi^{3}\ \alpha$] (d) at (3,0.5) ;
[ele,label=above:$\varphi^{4}\ \alpha$] (e) at (2.5,1) ;
[ele,label=above:$\varphi^{5}\ \alpha$] (f) at (1.5,1) ;
[ele,label=above:$\varphi^{6}\ \alpha$] (g) at (0.5,1) ;
[ele,label=above:$\varphi^{7}\ \alpha\qquad$] (h) at (0,0.5) ;

[->,dashed,shorten <=2pt,shorten >=2] (a) – (b);
[->,dashed,shorten <=2pt,shorten >=2] (b) – (c);
[->,dashed,shorten <=2pt,shorten >=2] (c) – (d);
[->,dashed,shorten <=2pt,shorten >=2] (d) – (e);
[->,dashed,shorten <=2pt,shorten >=2] (e) – (f);
[->,dashed,shorten <=2pt,shorten >=2] (f) – (g);
[->,dashed,shorten <=2pt,shorten >=2] (g) – (h);

[thick,color=purple] (f) to [bend left] node[midway,right] g (b);
[thick,color=purple] (g) to [bend left] node[midway,below] g (a);
[thick,color=purple] (e) to [bend left] node[midway,right] g (c);
[thick,color=purple] (d) to [out=50, in=-50,looseness=50]
node[midway,right] g (d);

[thick,color=blue] (a) to [out=130, in=-130,looseness=50]
node[midway,left] $f$ (a);
[thick,color=blue] (d) to [bend right] node[midway,above] $f$ (e);
[thick,color=blue] (c) to [bend right] node[midway,right] $f$ (f);
[thick,color=blue] (b) to [bend right] node[midway,left]  $f$ (g);

Orbit from an $f$ fixed point $\alpha$ with odd period $7$. Identical points on the left (marked by two parallel lines) are merged on the right (move white dot to black dot). In particular, on the left the two vertices of the shaded line are the same, forming a fixed point of $g$ on the right.

This figure shows an orbit from an f fixed point alpha with odd period 7. Identical points on the left, marked by two parallel lines, are merged on the right, by moving white dot to black dot. In particular, on the left the two vertices of the shaded line are the same, forming a fixed point of g on the right.

Moreover, the rightmost black dot and a preceding white dot from g must be the same, due to $f$-arcs from identical dots.
This means the half-period iterate, the rightmost black dot, must be a fixed point of g, say $\beta$.
If β = α, then period p = 1, in accordance with Theorem <ref>.
This example motivates the following:

When $f$ fixes $x$, and (f  g) has an odd period $p$ for $x$,
then g fixes (f  g), which is not $x$ itself if and only if p  1.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x

y = (f  g)  odd p

y  fixes g S  (y = x  p = 1)

First we show that g fixes $y$.
Let h = p div 2.
Since period $p$ is odd, $p \ee \HOLinline{\HOLNumLit{2}\HOLSymConst{\ensuremath{}}\HOLFreeVar{h}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{1}\;\HOLSymConst{=}\;\HOLFreeVar{h}\;\HOLSymConst{\ensuremath{+}}\;\HOLFreeVar{h}\;\HOLSymConst{\ensuremath{+}}\;\HOLNumLit{1}}$.
Thus h  h  1  0  p,
so y = g y by Theorem <ref>.
As (f  g) is a permutation, y  S, so y  fixes g S.
Theorem <ref> ensures that: y = x  p = 1.

§.§ Fixed Point Orbits

Let φ = f  g,
and $\alpha, \beta$ be fixed points of $f, \HOLinline{\HOLFreeVar{g}}$, respectively.
Theorem <ref> and Theorem <ref> show that:

* if the period $p$ of $\alpha$ is even, its orbit has another fixed point of $f$ at the h = p div 2 iterate: φ.

* if the period $p$ of $\alpha$ is odd, its orbit has another fixed point of g at the h = p div 2 iterate: β = φ.

Figure <ref> and Figure <ref> show that these orbits have no more fixed points. The only fixed point, of either $f$ or $g$, occurs at halfway point of the orbit.

Thus, fixed point orbits lead directly from one fixed point to another. This is because, assuming one of the intermediate iterate is a fixed point, the iteration path will turn back, due to either $f$ or $g$, both being involutions. This will produce an orbit with a shorter period, but period for an orbit is minimal.

Such considerations lead to the following stronger forms of Theorem <ref> and Theorem <ref>:

When $f$ fixes $x$, the $j$-th iterate of (f  g) from $x$ is a fixed point of either $f$ or g if and only if $j$ is half of the period $p$.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x  even p

j. 0  j  j  p

((f  g)  fixes f S  j = p div 2)

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x  odd p

j. 0  j  j  p

((f  g)  fixes g S  j = p div 2)

This completes our tour of the theory of permutation orbits and fixed points.
The results provide the key to formally prove that our two-squares algorithm by iterations is correct.

§ CORRECTNESS OF ALGORITHM

The algorithm to compute the flip fixed point from the known Zagier fixed point, given in Definition <ref>, makes use of a while-loop.
A while-loop consists of a guard $G$ and a body $B$, starting with an element $x$. The body is a function on $x$, producing iterates
$B(x)$, B, B, .
The guard is a predicate on each iterate: the loop continues only if the test result by the guard stays true.

In HOL4, the WHILE loop with guard $G$ and body $B$ starting with $x$ is defined as:

WHILE G B x  if G x then WHILE G B (B x) else x

from which one can easily show by induction that:

(j. j  k  G (B))

WHILE G B x =

if G (B) then WHILE G B (B) else B

giving this expected result:

The WHILE loop delivers the first body iterate that fails the guard test.

(j. j  k  G (B))  G (B)

WHILE G B x = B

§.§ Iterate with WHILE

From Section <ref>, we learn that for two involutions $f$ and g,
a fixed point $\alpha$ of $f$ is paired up with a fixed point $\beta$ of g
whenever the period of $\alpha$ under the composition φ = f  g is odd.
In fact, $\beta$ lies in the orbit of $\alpha$ at halfway point, the iterate at half period.
Since a while-loop also gives an iterate, we have:

For two involutions $f$ and g, if $f$ fixes $x$ with an odd period,
a WHILE loop with f  g from $x$ can reach a fixed point of g.

finite S  f involute S  g involute S

x  fixes f S  p = period (f  g) x  odd p

WHILE (t. g t  t) (f  g) x  fixes g S

Let guard G = (t. g t  t), and body B = f  g.
If period p = 1,
then x  fixes g S by Theorem <ref>.
So G x,
and WHILE G B x = x since the condition is not met at the start.
Therefore WHILE G B x  fixes g s.

If period p  1, let h = p div 2,
and z = B.
Since 1  p, $0 < h < p$.
Also $f$ and g are involutions, so B permutes S.
Hence z  fixes g S by Theorem <ref>.
so G z.

We claim j. j  h  G (B).
To see this, let y = B, which is an element of $S$.
If j = 0, then y = x.
Since period p  1,
y  fixes g S by Theorem <ref>, so G y.
If j  0, then $0 < j < h < p$, and j  h.
Hence y  fixes g S by Theorem <ref>,
so G y again. The claim is proved.

By the claim and G z,
apply Theorem <ref> to conclude
$\HOLinline{\HOLConst{WHILE}\;\HOLFreeVar{G}\;\HOLFreeVar{B}\;\HOLFreeVar{x}\;\HOLSymConst{=}\;\HOLFreeVar{z}} \in \HOLinline{\HOLConst{fixes}\;\HOLFreeVar{g}\;\HOLFreeVar{S}}$.

§.§ Two Squares by WHILE

We have developed the theory to show that the algorithm in Section <ref> is correct:

For a prime of the form 4k  1,
the two squares algorithm of Definiton <ref> gives a flip fixed point.

prime n  n  1  4

two_sq n  fixes flip (mills n)

Let S = mills n, φ = zagier  flip,
u = (1,1,n div 4), and period p = period φ u.
By Definition <ref>,
and noting that ()  found = (t. flip t  t),
this is to show: WHILE (t. flip t  t) φ u  fixes flip S.

Since a prime is not a square, we have finite S.
Now φ permutes S as Zagier map and flip map are both involutions,
by Theorem <ref> and Theorem <ref>,
and fixes zagier S = u by Theorem <ref>.
Thus u  fixes zagier S,
and period $p$ is odd by Corollary <ref>.
So WHILE (t. flip t  t) φ u  fixes flip S by Theorem <ref>.

It is almost trivial to convert two_sq n to following algorithm:

Compute the two squares for Fermat's two squares theorem.

two_squares n  (let (x,y,z) = two_sq n in (x,y  z))

giving the two squares in a pair,
and its correctness is readily demonstrated:

The algorithm by Definition <ref> gives indeed Fermat's two squares.

prime n  n  1  4

(let (u,v) = two_squares n in n = u  v)

Running Fermat's two squares algorithm in a HOL4 session, with timing.

This table shows a sample run of Fermat's two squares algorithm in a HOL4 session, with timing information.

Table <ref> shows a sample run in HOL4 session on a typical laptop,
using EVAL for evaluation and prefix time to obtain timing statistics.
Note that these EVAL executions are based on optimised symbolic rewriting in HOL4, thus orders of magnitude slower than running native code.

*Other algorithms
A prime has a finite set of windmill triples, by Theorem <ref>.
Fermat's two squares for a prime $n$ with n  1  4, which must exist by Theorem <ref>,
can be found by a brute-force search: subtract $n$ by successive odd squares, and check whether the difference is a square. Although there are better ways to test a square than the  square-root test, they are not simple to implement.

Don Zagier, after his one-sentence proof, referred to an effective algorithm by Wagon <cit.> to compute the two squares. The algorithm requires finding a quadratic non-residue of the given prime $n$.

The advantage of our algorithm in Definition <ref> over such alternative methods is that only addition and subtraction are performed.
The implementation is rather straightforward. The issue of termination is discussed next.

§.§ Terminating Condition

As mentioned in Section <ref>, for our algorithm the WHILE loop may or may not terminate. To gaurantee termination, convert the WHILE loop to a countdown loop, as follows. First, ensure that the input number $n$ is not a square, so that mills n is finite (Theorem <ref>), and check n  1  4, so that (1,1,n div 4)  mills n, , mills n  .

Obviously for any triple (x,y,z)  mills n, each $x, y$ or $y$ is less than $n$, hence mills n  n.

Now, use a countdown loop from n to $0$, start with the triple (1,1,n div 4) for the zagier  flip iteration.

The iterations trace an orbit. At half-way point, the orbit hits either a flip fixed point, detected by y = z, when the period is odd (Theorem <ref>), or another Zagier fixed point, detected by x = y, when the period is even (Theorem <ref>). They provide actual exits from the countdown loop, much earlier than the count drops to zero.

§.§ Lessons Learnt

This formalisation work can be a self-contained project in a theorem-proving workshop.
The ideas are simple, but formulating the theorems properly is not simple.
For example, at first the author would like to prove:

\begin{equation*}
\label{eqn:zagier-inv}
\HOLinline{\HOLSymConst{\HOLTokenForall{}}\HOLBoundVar{x}\;\HOLBoundVar{y}\;\HOLBoundVar{z}.\;\HOLConst{zagier}\;(\HOLConst{zagier}\;(\HOLBoundVar{x}\HOLSymConst{,}\HOLBoundVar{y}\HOLSymConst{,}\HOLBoundVar{z}))\;\HOLSymConst{=}\;(\HOLBoundVar{x}\HOLSymConst{,}\HOLBoundVar{y}\HOLSymConst{,}\HOLBoundVar{z})}.
\end{equation*}

The interactive session produces several subgoals which he cannot resolve immediately.
A comparison of Definition <ref> with Equation (<ref>)
shows differences in boundary cases.
Finally, some insight from windmills resolves why the boundaries are ignored, and provides the pre-condition x  0  z  0, see Equation (<ref>). The result is Theorem <ref>.

Explaining the Zagier map is an involution through the mind of a windmill poses some challenges. Definition <ref> of the Zagier map has $3$ branches, so the initial effort is to treat just $3$ cases. The first case is immediate, but the second case runs into a mess. It is only after drawing a lot of windmills that the author realises these finer points:

* there are $3$ types: $x < y, x \ee y$, and $x > y$ for the windmill triple $(x,y,z)$,

* the $3$ types are further subdivided due to geometry of the mind, giving $5$ cases in total,

* the $5$ cases can be condensed into $3$ branches, as the definition shows.

The result is Table <ref> in Section <ref>.

For the permutation orbits in Section <ref>,
the proofs about relations between iterates start as long-winded arguments treating if-part and only-if part separately.
Putting them in this paper prompts the author to rethink the logic.
The polished proofs simply employ a chain of logical equivalences.

Fixed point orbits have either even or odd period, as treated in Section <ref> and Section <ref>.
The drawing of the diagrams helps to refine the proofs to be short and sweet, making good use of theorems already proved.

About the correctness proof of the algorithm using a while-loop in Section <ref>, the author initially applied Hoare logic assertions to derive the desired iterate upon loop exit.
This is awkward, as pointed out by Michael Norrish who knows the HOL4 theorem-prover inside out.
This is awkward, as pointed out by (omitted for anonymous review).
The reason is that WHILE is defined as iteration of the body in HOL4.
The section had since been rewritten.

*Development Effort

The proofs have been streamlined after several revisions.
Such refinements result in the script line counts for various theories developed, shown in Table <ref>.

Statistics of various theories in this work.

This table gives the statistics of various theories in this formalisation work.

\[
\begin{array}{llr}
\text{HOL4 Theory}    & \text{Description}                        & \text{\#Lines}\\
\hline
\text{involute}       & \text{basic involution}                   & 231\\
\text{iteration}      & \text{function iteration and period}      & 917\\
\text{iterateCompose} & \text{iteration of involute composition}  & 1648\\
\text{iterateCompute} & \text{iteration period computation}       & 939\\
\text{windmill}       & \text{windmills and their involutions}    & 1844\\
\text{twoSquares}     & \text{two-squares by windmills}           & 1317\\
\end{array}
\]

The scripts are fully documented, including the traditional proofs as comment before each theorem. Although comments almost double the script size, the line counts are still indicative of the effort to convert ideas into formal proofs.

§.§ Related Work

As noted in Section <ref>, Fermat's two squares theorem has been formalised. However, none of these formal proofs is constructive, in the sense that there is no formal proof of an algorithm to compute the two squares for a prime satisfying the theorem.

Fermat's two squares theorem has two parts: existence and uniqueness.
All formal proofs include the existence part (see Theorem <ref>) , using classic and modern existence proofs: the method of infinite descent is used in one system, Gaussian integers are employed in three systems, both Heath-Brown's proof and Zagier's proof are treated in two systems.

Only two formal proofs include the uniqueness part (see Theorem <ref>): Théry <cit.> proved by algebraic identities and divisibility, and Hughes <cit.> proved by unique factorisation of Gaussian integers.

Recently, Dubach and Muehlboeck <cit.> formalised Zagier's proof using involutions in Coq's Mathematical Components Library. They illustrated their proof using the windmills as per this paper, and extended the use of involutions on the same set to formalise also an integer-partition proof of Fermat's two squares theorem by Christopher <cit.>.

A summary of these formal proofs, in chronological order, is given in Table <ref>.

Chronology of formalisation of Fermat's two squares theorem.

This table lists, in chronological order, the formal proofs of Fermat's two squares theorem, by various authors in different theorem provers.

Year     Author(s)[reference]     Theorem Prover     Comment

2004     Laurent Théry <cit.>     Coq     Gaussian integers, with uniqueness

2007     Roelof Oosterhuis <cit.>     Isabelle     Euler's proof with infinite descent

2009     Marco Riccardi <cit.>     Mizar     Heath-Brown's proof with involutions

2010     John Harrison <cit.>     HOL Light     Zagier's proof with involutions

2012     Anthony Narkawicz <cit.>     NASA PVS     Zagier's proof with involutions

2015     Mario Carneiro <cit.>     MetaMath     Gaussian integers

2016     Rob Arthan <cit.>     ProofPower     Heath-Brown's proof with involutions

2019     Chris Hughes <cit.>     Lean     Principal Ideal Ring of Gaussian integers, with uniqueness

2021     Dubach and Muehlboeck <cit.>     Coq     Zagier's and Christopher's proofs with involutions

§ CONCLUSION

About Fermat's two squares theorem,
G. H. Hardy wrote in his 1940 essay A Mathematician's Apology <cit.>:

This is Fermat's theorem, which is ranked, very justly, as one of the finest of arithmetic. Unfortunately, there is no proof within the comprehension of anybody but a fairly expert mathematician.

This work has been a rewarding exercise in formalisation,
delivering a proof of Fermat's Theorem <ref> using only natural numbers, involutions, and counting.
There is a certain sense of mathematical beauty when a non-trivial result can be shown by elementary means, borrowing elegant ideas by Zagier and Spivak.
Moreover, by developing a theory of involution iteration, an algorithm to compute the two squares of the theorem can be formally shown to be correct.

*Future Work
The theory in Section <ref>, about orbits and fixed points, can be developed using group actions, since the iteration indices form an addition cyclic group under mod $p$, where $p$ is the orbit period.
One can exploit the symmetry in permutation orbits, especially for permutations arising from two involutions, to improve the algorithm, as shown in the analysis by Shiu <cit.>.
In HOL4, this direction can start from the algebra of group theory in Chan and Norrish <cit.>.
A formal analysis of the performance of the algorithm for two squares described in Definition <ref> can be modelled using an approach in Chan <cit.>.

§ ACKNOWLEDGEMENTS

Many thanks to Michael Norrish for his careful review of the draft, providing useful advice and helpful recommendations to improve this paper.
The author is also grateful to the anonymous reviewers who pointed out typographical errors and suggested clarifications.
This paper has been revised to incorporate their comments.

§ APPENDICES

§ CROSS-REFERENCE WITH PROOF SCRIPT

This is supplementary material for review.
Refer to Table <ref> for a
mapping of theorems in this paper with names in scripts.

Cross-reference of theorems and definitions. HOL4 proof scripts uploaded as script.zip.

This is a table giving the cross-reference of therems and definitions in the scripts. HOL4 proof scripts have been uploaded as script.zip file.

Location of $\dots$     Script: (name)Script.sml
Name of theorem or definition     Line #

Theorem <ref>
twoSquares     fermat_two_squares_iff     619

Definition <ref>
windmill     windmill_def     348

Definition <ref>
windmill     mills_def     481

Theorem <ref>
windmill     mills_finite_non_square     714

Theorem <ref>
windmill     windmill_trivial_prime     428

Definition <ref>
involuteFix     pairs_def     269

Definition <ref>
involuteFix     fixes_def     264

Theorem <ref>
involuteFix     involute_two_fixes_both_odd     1182

Definition <ref>
windmill     flip_def     884

Theorem <ref>
windmill     flip_involute_mills     933

Definition <ref>
windmill     zagier_def     959

Theorem <ref>
windmill     zagier_involute_mills_prime     1475

Theorem <ref>
twoSquares     zagier_fixes_prime     162

Theorem <ref>
twoSquares     fermat_two_squares_exists_odd_even     441

Theorem <ref>
twoSquares     fermat_two_squares_unique_thm     205

Theorem <ref>
helperTwosq     mod_4_not_squares     419

Definition <ref>
twoSquares     two_sq_def     816

Definition <ref>
iteration     iterate_period_def     348

Theorem <ref>
iteration     iterate_period_mod     653

Theorem <ref>
iterateCompose     involute_involute_fixes_both     558

Theorem <ref>
iterateCompose     iterate_involute_mod_period     401

Theorem <ref>
iterateCompose     involute_involute_fix_orbit_1     583

Theorem <ref>
iterateCompose     involute_involute_fix_orbit_2     689

Theorem <ref>
iterateCompose     involute_involute_fix_orbit_fix_even_distinct     884

Corollary <ref>
iterateCompose     involute_involute_fix_sing_period_odd     1009

Theorem <ref>
iterateCompose     involute_involute_fix_orbit_fix_odd_distinct_iff     980

Theorem <ref>
iterateCompose     involute_involute_fix_even_period_fix     1100

Theorem <ref>
iterateCompose     involute_involute_fix_odd_period_fix     1146

Theorem <ref>
iterateCompute     iterate_while_thm     922

Theorem <ref>
iterateCompose     involute_involute_fix_odd_period_fix_while     1536

Theorem <ref>
twoSquares     two_sq_thm     840

Definition <ref>
twoSquares     two_squares_def     1025

Theorem <ref>
twoSquares     two_squares_thm     1041

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

authorpersonRob Arthan.
titleMathematical Case Studies: Some Number Theory.
howpublishedSupplementary ProofPower Examples.
Available from
<http://www.lemma-one.com/ProofPower/examples/wrk074.pdf>, accessed 22
August, 2021.

authorpersonMario Carneiro.
titleTheorem 2sq.
howpublishedMetamath 100 proof.
Available from <http://us.metamath.org/mpeuni/2sq.html>, accessed
22 August, 2021.

authorpersonHing Lun Chan.
titlePrimality Testing is Polynomial-time: a
Mechanised Verification of the AKS Algorithm.
PhD. schoolAustralian National University,
addressCanberra, Australia.

[Chan and NorrishChan and Norrish2012]
authorpersonHing Lun Chan and
personMichael Norrish. year2012.
A String of Pearls: Proofs of Fermat's Little
Theorem. In booktitleProceedings of Certified Programs and
Proofs (seriesLNCS, number7679),
editorpersonChris Hawblitzel and
personDale Miller (Eds.). publisherSpringer,

authorpersonA. David Christopher.
A partition-theoretic proof of Fermat's Two
Squares Theorem.
journalDiscrete Mathematics volume339,
number4 (date06 April year2016),

[Dubach and MuehlboeckDubach and
authorpersonGuillaume Dubach and
personFabian Muehlboeck. year2021.
Formal verification of Zagier's one-sentence
(date24 April year2021).
In arXiv: <https://arxiv.org/abs/2103.11389>, accessed 22 August,

authorpersonGodfrey H. Hardy.
booktitleA Mathematician's Apology.
publisherCambridge University Press.

authorpersonJohn Harrison.
titleRepresentation of primes $\equiv 1\pmod{4}$ as sum
of $2$ squares.
howpublishedFormalizing 100 Theorems in HOL Light.
Available from
accessed 22 August, 2021.

authorpersonRoger Heath-Brown.
Fermat’s two squares theorem.
journalInvariant  volume11
(year1984), pages3–5.
Available from <https://core.ac.uk/reader/97080>, accessed 22
August, 2021.

authorpersonChris Hughes.
titleSums of two squares.
howpublishedLean Community, 100 theorems.
Available from <https://leanprover-community.github.io/100.html>,
accessed 22 August, 2021.

[Malter, Schleicher, and ZagierMalter
et al2013]
authorpersonAimeric Malter, personDierk
Schleicher, and personDon Zagier.
New Looks at Old Number Theory.
journalThe American Mathematical Monthly
volume120, number3 (year2013),

authorpersonAnthony Narkawicz.
titleprimes_sum_squares: THEORY.
howpublishedNASA PVS Library 6.0.9.
accessed 22 August, 2021.

authorpersonRoelof Oosterhuis.
titleTheory TwoSquares.
howpublishedIsabelle, Archive of Formal Proofs.
Available from <https://www.isa-afp.org/entries/SumSquares.html>,
accessed 22 August, 2021.

authorpersonMarco Riccardi.
titletheorem :: NAT_5:23.
howpublishedFormalizing 100 Theorems in Mizar.
Available from <http://www.mizar.org/100/index.html>, accessed 22
August, 2021.

authorpersonPeter Shiu.
Involutions associated with Sums of Two Squares.
journalPublications de l'Institut
Mathématique  volume59 (year1996),

authorpersonAlexander Spivak.
Winged Squares [in Russian].
journalPopular lectures in mathematics, 2006-2007
academic year, Lecture 15 (165) (date10 March
Available from <http://mmmf.msu.ru/lect/spivak/zagir_!.pdf>,
accessed 22 August, 2021.

authorpersonLaurent Théry.
titleNumbers Equal to the Sum of Two Square Numbers.
howpublishedFormalizing 100 theorems in Coq.
Available from
<https://github.com/coq-contribs/sum-of-two-square>, accessed 22 August,

[Trimble and LamaTrimble and Lama2008]
authorpersonTodd Trimble and personVishal
Lama. year2008.
titleA proof from The Book?
howpublishedTodd and Vishal's blog.
Available from
accessed 22 August, 2021.

authorpersonStan Wagon.
Editor's Corner: The Euclidean Algorithm Strikes
journalThe American Mathematical Monthly
volume97, number2 (year1990),

authorpersonFreek Wiedijk. yearlast update:
titleFormalizing 100 Theorems.
Available from <https://www.cs.ru.nl/ freek/100/>, accessed 22
August, 2021.

authorpersonKenneth S. Williams.
booktitleNumber Theory in the Spirit of
publisherCambridge University Press, Cambridge, UK.

ISBN: 9780511751684.

authorpersonDon Zagier.
A One-Sentence Proof That Every Prime $p\equiv
1\pmod{4}$ Is a Sum of Two Squares.
journalThe American Mathematical Monthly
volume97, number2 (year1990),