Windmills of the Minds: A Hopping Algorithm for Fermat’s Two Squares Theorem

Abstract

Fermat’s two squares theorem asserts that a prime one more than a multiple of 4 is a sum of two squares. There are many proofs of this gem in number theory, including a remarkable one-sentence proof by Don Zagier based on two involutions on a finite set built from such a prime. Applying the two involutions alternatively leads to an iterative algorithm to find the two squares for the prime. Moreover, a detailed analysis of the computation reveals that it is possible to jump through the iteration nodes, leading to a better hopping algorithm. Here is a formalisation of Zagier’s proof, deriving the involutions using windmill patterns. Theories developed for the formal proof are used to establish the correctness of both algorithms.

1 Introduction

The only even prime ${{2}\,{=}\,{1}^{2}\,{{+}}\,{1}^{{2}}}$, is indeed a sum of two squares. An odd prime, upon division by 4, leaves a remainder of either 1 or 3. We shall invent names for these two types:

• those of the form $4k + 1$ for some k will be called tiks,

• those of the form $4k + 3$ for some k will be called toks.

While a tok prime can never be expressed as a sum of two squares,¹ a tik prime seems always possible, as supported by numerical evidence from Table 1.

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

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

Theorem 1

(Two Squares Theorem) [script]

A tik prime is a sum of two squares of opposite parity.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {prime} \,{n}\,{\wedge }\,\texttt {tik} \,{n}\,{\Rightarrow }\,{\exists \,}({u}{,}{v}).\,\texttt {{odd}} \,{u}\,{\wedge }\,\texttt {{even}} \,{v}\,{\wedge }\,{n}\,{=}\,{u}^{2}\,{{+}}\,{v}^{2} \end{array} \end{aligned}$$

Fermat corresponded with others mathematicians, claiming this assertion without proof. In a rare instance, he outlined his method in a 1659 letter to Pierre de Carcavi, which he called infinite descent. Leonhard Euler was the first to work out the details of Fermat’s method after much effort, and published a proof in 1749.

Although this theorem is proved, other mathematicians offer alternative proofs from time to time. Clearly, the consensus is that the known proofs are logically valid but not satisfactory in some way. It is hard to pinpoint exactly what that is, but Paul Erdős believed that God carries a Book containing the best proof of each mathematical theorem.²

As a dedicated work to honor Paul Erdős, Martin Aigner and Günter M. Ziegler edited “Proofs from THE BOOK”, including Fermat’s Two Squares Theorem as one of the topics [4, Chapter 4]. The chapter explains a proof with illustrations based on involution (a self-inverse function) and fixed points (those values unchanged by the function). Don Zagier condensed such a proof into one sentence [18] in 1990:

The involution on the finite set

$$\begin{aligned} S = \{(x,y,z) \in \mathbb {N}^{3} \mid {{n}\,{=}\,{x}^{2}\,{{+}}\,{4}{y}{z}} \} \end{aligned}$$

defined by

$$\begin{aligned} {({x}{,}{y}{,}{z})} \longmapsto {\left\{ \begin{array}{ll} (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{array}\right. } \end{aligned}$$

(1)

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

The author learned about the windmill patterns in Zagier’s proof from a video by Mathologer [11], became interested in an algorithm found in [2, 9, 14] based on the ingredients in Zagier’s proof, and adapted an improvement of the algorithm by Shiu [15] after an analysis of the computations involved.

This work aims at a coherent presentation of a formalisation of Zagier’s proof, leading to algorithms to find the two squares in Fermat’s Theorem for a tik prime, with formal proofs of the correctness of algorithms.

1.1 Contribution

This paper extends my earlier conference paper [3] which only covers the iterative algorithm. Now the emphasis is about an improved hopping algorithm, how it relates to the iterative algorithm and why hopping is possible. Although the performance boost is slight, the ideas behind the improvement are non-trivial. The effort leads to the first formal proof that the hopping algorithm is correct.

Note that Zagier’s proof has been formalised, in HOL Light [7], in NASA PVS [10] and in Coq [6], without applying the involutions to construct an algorithm.

As mentioned before, ideas in this paper are taken from various sources, with the algorithmic part mostly from Zagier [9] and Shiu [15]. The novel feature of this work is an elegant and pictorial approach in the formalisation, by building the relevant underlying theories.

Formalising an algorithm and showing its correctness involve proving that the algorithm terminates, and proving that when the algorithm terminates, it delivers the desired result. Getting these done means providing formal definitions and developing appropriate theories.

The current work develops general theories about involutions and fixed points, permutations from composition of involutions, and permutation orbits in a finite set. These theories are not only useful for the present work, but will also support further work, e.g., are there still improvements to the hopping algorithm?

1.2 Overview

This work of formalisation features the following:

• introduce geometric ideas to construct two involutions on a finite set in Sect. 2,

• explain an algebraic proof of Fermat’s theorem by involutions and fixed points in Sect. 3,

• compute the two squares by an algorithm, with a brief discussion of the theory behind, in Sect. 4,

• devise an improved algorithm to compute the two squares with correctness proof in Sect. 5.

Finally, some concluding remarks are given in Sect. 6.

1.3 Notation

Statements starting with a turnstile ($\vdash$) are HOL4 theorems, automatically pretty-printed to LaTeX from the relevant theory in HOL4 development. Generally, the notation allows an appealing combination of quantifiers ($\forall , \exists , \exists {!}$), logical connectives ($\wedge$ for “and”, $\vee$ for “or”, $\lnot$ for “not”, also $\Rightarrow$ for “implies” and $\iff$ for “if and only if”), set theory ($\in$ for “an element of”, and comprehensions such as ${\left\{ {x}\,\texttt {|}\,{x}\,{{<}}\,{6}\right\} }$), and functional programming ($\lambda \,$ for abstraction, and juxtaposition for application). Repeated application of a function f is indicated by exponents, e.g., ${{f}\,({f}\,({f}\,{x}))\,{=}\,{f}^{3}({x})}$, and ${{f}\,{\circ }\,{g}}$ indicates the composition of functions f and g from right to left.

For a function f from set S to set T, we write ${{f}\,{:}\,{S}\,{\leftrightarrow }\,{T}}$ to mean a bijection. The empty set is denoted by ${{\emptyset }}$, and a finite set, denoted by ${\texttt {{finite}}\,{S}}$, has cardinality |S|.

The set of natural numbers is denoted by $\mathbb {N}$, counting from 0, and we define , where means ‘equality by definition’. For $n \in \mathbb {N}$, ${\texttt {square}\,{n}}$ means it is a square: ${{\exists \,}{k}.\,{n}\,{=}\,{k}^{2}}$, ${\texttt {prime}\,{n}}$ means it is a prime, and ${\texttt {{even}}\,{n}}$ or ${\texttt {{odd}}\,{n}}$ denotes its parity. When integer m is divided by n, the integer quotient and integer remainder are written as ${{m}\,\texttt {{div}}\,{n}}$ and ${{m}\,\texttt {{mod}}\,{n}}$, respectively. We write ${{n}\,{{\mid }}\,{m}}$ when n divides m, which is equivalent to ${{m}\,{\equiv }\,{0}\,{(}{\bmod }\,{n}{)}}$ when ${{n}\,{\ne }\,{0}}$. Note that ${\texttt {{sqrt}}\,{n}}$ denotes the integer square-root of n, e.g., ${\texttt {{sqrt}}\,{13}\,{=}\,{3}}$.

A list with members a, b, c is denoted by ${[{a};\,{b};\,{c}]}$, and the empty list is []. For a non-empty list ls with length $|ls |\ne 0$, its head is ${\texttt {HD}\,{ls}}$, its last is ${\texttt {LAST}\,{ls}}$, and ${\texttt {FRONT}\,{ls}}$ drops the last, i.e., ${{ls}\,{=}\,\texttt {FRONT}\,{ls}\,{\mathbin {+\hspace{-5.55542pt}+}}\,[\texttt {LAST}\,{ls}]}$ where $\mathbin {+\hspace{-5.55542pt}+}$ denotes list concatenation. We write ${\texttt {ALL\_DISTINCT}\,{ls}}$ for a list with distinct elements, ${\texttt {FILTER}\,{P}\,{ls}}$ is a list with only those elements of ls satisfying the predicate P, and ${\texttt {MAP}\,{f}\,{ls}}$ forms a list by applying the function f to every element of the list ls.

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 f49cbe99d286. In this paper, each theorem has a mark [script], which is hyperlinked to the appropriate line of the corresponding proof script in repository.

2 Geometric Ideas

Let n be an odd number that can be expressed as a sum of two squares: $n = x^{2} + y^{2}$. Note that x and y must be of opposite parity. Let x be odd, then $y = 2z$ for some z. Thus $n = x^{2} + 4z^{2}$. Zagier’s proof considers such representations of n, which have nice geometric pictures.

Fig. 1

Typical windmills, where ${\texttt {windmill}\,({x}{,}{y}{,}{z})\,{=}\,{x}^{2}\,{{+}}\,{4}{y}{z}}$ (Note rightmost one has ${{y}\,{=}\,{z}})$

2.1 Windmills

Suppose a tik n is a sum of two squares, then $n = 4k + 1 = 1^{2} + 4(1)(k) = x^{2} + 4(z)(z)$. Both representations can be unified by a pattern:

Definition 2

A windmill has a central square with four identical rectangular arms.

Figure 1 shows some typical windmills, with ${{x}^{2}}$ as the central square and 4y z as four arms, each a rectangle of width y and height z, arranged clockwise around the square. If a tik n has a windmill with square arms, the four identical squares form a big square, so a picture of n as sum of two squares emerges (rightmost windmill in Fig. 1).

To investigate this possibility, first collect all the windmills representing a number:

Definition 3

The mills of a number is its set of windmills.

Note that ${\texttt {mills}\,{n}}$ is the set S given in Zagier’s one-sentence proof. As long as n is not a square, its windmill, if any, has arms, making ${\texttt {mills}\,{n}}$ finite:

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {{finite}}\,(\texttt {mills}\,{n})\,{\iff }\,{\lnot } \texttt {square}\,{n} \end{array} \end{aligned}$$

Furthermore, the number of windmills is bounded by the non-square n (refer to a proof in [script]):

It is straight-forward to determine all windmill triples (x, y, z) of non-square ${\texttt {tik}\,{n}}$ by taking successive odd x and getting factors for the product ${{y}{z}\,{=}\,({n}\,{{-}}\,{x}^{2})\,\texttt {{div}}\,{4}}$, as shown in Table 2 for $n = 29$. The corresponding windmills are shown in Fig. 2.

Table 2 Determine all the windmill triples (x,y,z) of ${{n}\,{=}\,{29}}$, by odd x and factors of yz

Fig. 2

All the windmills of ${{n}\,{=}\,{29}}$ determined from Table 2

A tik ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$ has ${({1}{,}{1}{,}{k})\,{\in }\,\texttt {mills}\,{n}}$, where ${{k}\,{=}\,{n}\,\texttt {{div}}\,{4}}$. Moreover, for a tik prime, this is the only windmill (x, y, z) with ${{x}\,{=}\,{y}}$:

Theorem 4

[script]

For a prime of the form ${{4}{k}\,{{+}}\,{1}}$, the only windmill with the first and second parameters equal is ${\texttt {windmill} \,({1}{,}{1}{,}{k})}$.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {tik} \,{n}\,{\wedge }\,\texttt {prime} \,{n}\,{\Rightarrow }\,{\forall \,}{x}\,{z}.\,{n}\,{=}\,\texttt {windmill} \,({x}{,}{x}{,}{z})\,{\iff }\,{x}\,{=}\,{1}\,{\wedge }\,{z}\,{=}\,{n}\,\texttt {{div}} \,{4} \end{array} \end{aligned}$$

A proof of this theorem has been given in Chan [3, Theorem 2.5], due to the fact that x divides ${{n}\,{=}\,\texttt {windmill}\,({x}{,}{x}{,}{z})}$.

2.2 Flip Map

In the set of windmills, pairs with arms flipped have their triples swapping width and height, leaving the central square unchanged:

Definition 5

The flip map for a triple.

Obviously the inverse of flip is itself, for any triple t:

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {flip}\,(\texttt {flip}\,{t})\,{=}\,{t} \end{array} \end{aligned}$$

A self-inverse map is called an involution, which we shall define and study later (Sect. 3.1). Meanwhile, there is another intriguing involution relating pairs of windmills.

2.3 Mind of a Windmill

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

Referring to Fig.3., a windmill has a mind (marked in dashes at middle), which is the maximum central square that can be fitted within the outline shape. Through the mind, the leftmost and rightmost windmills are paired up. They have the same mind, although each has a central square of side x and $x'$, respectively. When , going from left to right, x can grow to $x'$, by expanding the central square but shrinking the arms. Conversely, going from right to left, $x'$ can contract to x by shrinking the central square but extending the arms. Indeed, transforming a windmill through the mind is another involution, as we shall see.

Fig. 3

A typical ${\texttt {windmill}\,({x}{,}{y}{,}{z})\,{=}\,{x}^{2}\,{{+}}\,{4}{y}{z}}$ with a mind (in dashes) and transforms to another windmill

To see the mind properly, we need to consider cases of x, y, z in a windmill (x, y, z). There are five cases, depending on the geometric pattern of the windmill. For a detailed treatment of these five cases with pictures, see Chan [3, Sect. 3.3]. Here we only show the summary in Table 3.

Table 3 The five cases of transforming a triple (x, y, z) to $(x',y',z')$ through the mind

Note that the transformation rule:

$$\begin{aligned} (x',y',z') \ \texttt {=}\ (2y - x, y, x + z - y) \end{aligned}$$

happens to be the same for case 2 and case 4. The same rule actually applies to case 3, which has ${{x}\,{=}\,{y}}$. Taking only the distinct cases, the mind of a windmill triple is:

2.4 Zagier Map

Table 3 shows the transformation rule from (x, y, z) to $(x',y',z')$ for windmills of the same mind. Recall that cases 2, 3 and 4 have the same rule, so the map can be succinctly expressed with only three branches:

Definition 6

The zagier map for a triple.

This map is essentially the “magic” one devised by Zagier in Eq. (1). Due to the else-parts, there is a slight difference in the boundary cases, which are irrelevant as they cannot occur for proper windmills, i.e., windmills from non-square tik numbers³.

Referring again to Table 3, the windmills in ${{S}\,{=}\,\texttt {mills}\,{n}}$ can be partitioned into three triple types:

$$\begin{aligned} \begin{array}{cl} S_{x< y} \ \texttt {=}\ \{(x,y,z) \in S \mid x < y\} & \hbox {covering cases} 1\, \hbox {and}\, 2\\ S_{x = y} \ \texttt {=}\ \{(x,y,z) \in S \mid x \ \texttt {=}\ y\} & \hbox {covering case}\ 3\\ S_{x> y} \ \texttt {=}\ \{(x,y,z) \in S \mid x > y\} & \hbox {covering cases} \,4\, \hbox {and}\, 5\\ \end{array} \end{aligned}$$

and the effect of zagier map for each type is:

• 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 windmills with central square and arms:

$$\begin{aligned} \begin{array}{l} \,\,\vdash {x}\,{\ne }\,{0}\,{\wedge }\,{z}\,{\ne }\,{0}\,{\Rightarrow }\,\texttt {zagier}\,(\texttt {zagier}\,({x}{,}{y}{,}{z}))\,{=}\,({x}{,}{y}{,}{z}) \end{array} \end{aligned}$$

(2)

From these definitions, it is simple to verify that the mind is an invariant under the zagier map for any windmill triple:

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {mind}\,(\texttt {zagier}\,({x}{,}{y}{,}{z}))\,{=}\,\texttt {mind}\,({x}{,}{y}{,}{z}) \end{array} \end{aligned}$$

For the windmills of ${{n}\,{=}\,{29}}$, pairing by the zagier map is shown in Fig. 4.

Fig. 4

Partition of windmills of ${{n}\,{=}\,{29}}$ for zagier: those with $x \ \texttt {=}\ y, x < y$, and $x > y$ (Note pairing by minds through arrows)

3 Algebraic Proof

The maps flip and zagier are examples of involution on ${\texttt {mills}\,{n}}$, the set of windmills for n. With these maps defined, the proof of Theorem 1 can proceed after studying involutions.

3.1 Involution

A self-inverse function f is an involution on a set S, denoted by${{f}\,\texttt {involute}\,{S}}$:

In other words, ${{f}\,{:}\,{S}\,{\leftrightarrow }\,{S}}$ is a bijection, pairing up x and ${{f}\,{x}}$, both in S. An element x is fixed by the involution f when ${{x}\,{=}\,{f}\,{x}}$, leading to the following disjoint sets:

Evidently ${\texttt {pairs}\,{f}\,{S}}$ consists of distinct involute pairs, so its cardinality is even:

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {{finite}}\,{S}\,{\wedge }\,{f}\,\texttt {involute}\,{S}\,{\Rightarrow }\,\texttt {{even}}\,{|} \texttt {pairs}\,{f}\,{S}{|} \end{array} \end{aligned}$$

As a result, the parity of |S| depends only on the parity of ${{|} \texttt {fixes}\,{f}\,{S}{|}}$, giving:

Theorem 7

[script]

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

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {{finite}} \,{S}\,{\wedge }\,{f}\,\texttt {involute} \,{S}\,{\wedge }\,{g}\,\texttt {involute} \,{S}\,{\Rightarrow }\\ \,\,\,\,\,\,\,(\texttt {{odd}} \,{|}\texttt {fixes} \,{f}\,{S}{|}\,{\iff }\,\texttt {{odd}} \,{|}\texttt {fixes} \,{g}\,{S}{|}) \end{array} \end{aligned}$$

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

3.2 Two Squares Theorem

As indicated in Table 3, for the zagier map only a triple of case 3 can map to itself:

$$\begin{aligned}\begin{array}{l} \,\,\vdash {x}\,{\ne }\,{0}\,{\Rightarrow }\,(\texttt {zagier}\,({x}{,}{y}{,}{z})\,{=}\,({x}{,}{y}{,}{z})\,{\iff }\,{x}\,{=}\,{y}) \end{array} \end{aligned}$$

Hence $S_{x \ \texttt {=}\ y} \ \texttt {=}\ {\texttt {fixes}\,\texttt {zagier}\,(\texttt {mills}\,{n})}$. Recall Theorem 4 which characterises such triples (x, x, z), we have:

Theorem 8

[script]

A prime of the form ${{4}{k}\,{{+}}\,{1}}$ has only (1, 1, k) fixed by the zagier map.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {prime} \,{n}\,{\wedge }\,\texttt {tik} \,{n}\,{\Rightarrow }\,\texttt {fixes} \,\texttt {zagier} \,(\texttt {mills} \,{n})\,{=}\,\left\{ ({1}{,}{1}{,}{n}\,\texttt {{div}} \,{4})\right\} \end{array} \end{aligned}$$

This is crucial to show the existence of two squares for a tik prime:

Proof of Theorem 1

Let n be a tik prime. A prime is not a square, so the set of windmills ${{S}\,{=}\,\texttt {mills}\,{n}}$ is finite. Both the zagier and flip maps are involutions on S. By Theorem 8, the zagier map has a single fixed point, so ${{|} \texttt {fixes}\,\texttt {zagier}\,{S}{|}\,{=}\,{1}}$. Thus ${{|} \texttt {fixes}\,\texttt {flip}\,{S}{|}}$ is odd by Theorem 7. Hence ${\texttt {fixes}\,\texttt {flip}\,{S}\,{\ne }\,{\emptyset }}$, containing a flip fix triple (x, z, z), so that ${{n}\,{=}\,\texttt {windmill}\,({x}{,}{z}{,}{z})} \ \texttt {=}\ {{x}^{2}\,{{+}}\,{4}{z}^{2}}$. Take ${{u}\,{=}\,{x}}$, and ${{v}\,{=}\,{2}{z}}$, then ${{n}\,{=}\,{u}^{2}\,{{+}}\,{v}^{2}}$. Obviously v is even, and u is odd since tik n is odd. $\square$

This proof of Fermat’s two squares theorem is purely algebraic. Indeed, in his 1990 one-sentence proof [18], Don Zagier did not provide any justification for the zagier map. There had been attempts to derive the map algebraically, e.g., Dijkstra [5] in 1993. A “winged squares” interpretation was later given by Spivak [16] in 2007. Although geometric pictures are not necessary in Zagier’s proof, the windmills do shed some lights on the nature of the zagier involution. Using elementary number theory, the two squares for a tik prime can be shown to be unique (see Chan [3, Theorem 4.3]):

Zagier’s one-sentence proof assumes that readers are familiar with, or at least are able to construct on-the-fly, the basic theory of involution on a finite set (Sect. 3.1). Formalising the proof means building up the involution theory in a theorem-prover. This is done by the author in HOL4, but any theorem-prover supporting sets and functions will do. Current formalisations of Zagier’s proof include HOL Light [7], NASA PVS [10] and Coq [6]. Zagier’s proof is an abridged form of Heath-Brown’s proof [8], which has been formalised in Mizar [13] and ProofPower [1]. There are other formalisations of Fermat’s two squares theorem, not based on Zagier’s idea. Please refer to Chan [3, Sect. 7.5, Table 6] for more information.

4 Basic Algorithm

The triple fixed by flip map in the set of windmills for a tik prime n gives a sum of two squares for n. Zagier’s proof asserts its existence, we shall construct an algorithm to compute the two squares.

Fig. 5

The iteration chain of ${{n}\,{=}\,{29}}$ by the composition ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ from zagier fix to flip fix

4.1 Iterative Algorithm

Let ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$ be a tik prime. By Theorem 8, the only zagier fixed point is ${{u}\,{=}\,({1}{,}{1}{,}{k})}$, meaning ${\texttt {zagier}\,{u}\,{=}\,{u}}$. To alter the triple u, applying flip is the obvious choice. Next application should be zagier to keep altering the triple, and so on alternatively. Thus by applying the composition ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ repeatedly from the known zagier fixed point, there is hope that the chain will lead to the only flip fixed point. Figure 5 shows that this is indeed the case for $\, {{n}\,{=}\,{29}}$.

For a tik n, once the flip fixed triple (x, z, z) is found, ${{n}\,{=}\,{x}^{2}\,{{+}}\,({2}{z})^{2}}$ as shown in the proof in Sect. 3.2, giving the two squares in Fermat’s theorem. This suggests the following iterative algorithm to locate the flip fix:

Since HOL4 supports while-loops, this algorithm can be implemented closely following the pseudo-code⁴:

Definition 9

Computing the flip fixed point of ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$ using a WHILE loop.

Details of a formal proof of its correctness are given in Chan [3, Sect. 7.2]. The proof is based on a theory of permutation orbits.

4.2 Permutation Orbits

In general, the composition of two involutions is just a permutation, not an involution. Let $\varphi :S \rightarrow S$ be a permutation, a bijection on the set S. For an element ${{x}\,{\in }\,{S}}$ the iteration nodes $\varphi (x)$, ${{\varphi }^{2}({x})}$, ${{\varphi }^{3}({x})}$, etc., form its orbit. The period of x under $\varphi$ is the smallest positive index n such that ${{\varphi }^{n}({x})\,{=}\,{x}}$. If such a positive index does not exist, the period is defined to be 0. For a finite set S, the period of any orbit is always nonzero.

When the permutation ${{\varphi }\,{=}\,{f}\,{\circ }\,{g}}$ is a composition of two involutions f and g, the theory developed in Chan [3, Sect. 6] shows that fixed points of either f or g are connected by a chain of composition iteration nodes. The result depends on the parity of the period for the orbit from a fixed point:

Theorem 10

[script]

For the orbit of $({{f}\,{\circ }\,{g}})$ from a fixed point of f, the node at half of the period p is another fixed point, either of f for even p, or of g for odd p.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {{finite}} \,{S}\,{\wedge }\,{f}\,\texttt {involute} \,{S}\,{\wedge }\,{g}\,\texttt {involute} \,{S}\,{\wedge }\,{x}\,{\in }\,\texttt {fixes} \,{f}\,{S}\,{\wedge }\,{p}\,{=}\,\texttt {{period}} \,({f}\,{\circ }\,{g})\,{x}\,{\wedge }\,\texttt {{even}} \,{p}\,{\Rightarrow }\\ \,\,\,\,\,\,\,{\forall \,}{j}.\,{0}\,{{<}}\,{j}\,{\wedge }\,{j}\,{{<}}\,{p}\,{\Rightarrow }\,(({f}\,{\circ }\,{g})^{j}({x})\,{\in }\,\texttt {fixes} \,{f}\,{S}\,{\iff }\,{j}\,{=}\,{p}\,\texttt {{div}} \,{2})\\ \,\,\vdash \texttt {{finite}} \,{S}\,{\wedge }\,{f}\,\texttt {involute} \,{S}\,{\wedge }\,{g}\,\texttt {involute} \,{S}\,{\wedge }\,{x}\,{\in }\,\texttt {fixes} \,{f}\,{S}\,{\wedge }\,{p}\,{=}\,\texttt {{period}} \,({f}\,{\circ }\,{g})\,{x}\,{\wedge }\,\texttt {{odd}} \,{p}\,{\Rightarrow }\\ \,\,\,\,\,\,\,{\forall \,}{j}.\,{j}\,{{<}}\,{p}\,{\Rightarrow }\,(({f}\,{\circ }\,{g})^{j}({x})\,{\in }\,\texttt {fixes} \,{g}\,{S}\,{\iff }\,{j}\,{=}\,{p}\,\texttt {{div}} \,{2}) \end{array} \end{aligned}$$

For a tik prime ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$, the only zagier fix (1, 1, k) starts an orbit in set of windmills ${\texttt {mills}\,{n}}$. The algorithm in Sect. 4.1 goes through a node in the orbit during each iteration of the while-loop. Since the zagier fix is unique by Theorem 8, the first part of preceding Theorem 10 ensures that the period is odd. The second part of Theorem 10 guarantees that the while-loop will exit when the iteration reaches the flip fix at the halfway node of the orbit.

It is easy to convert ${\texttt {two\_sq}\,{n}}$ to the following algorithm:

Definition 11

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

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

Theorem 12

[script]

The algorithm by Definition 11 gives Fermat’s two squares.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {prime} \,{n}\,{\wedge }\,\texttt {tik} \,{n}\,{\Rightarrow }\,({{{\textbf {{\textsf {let}}}}}}\,({u}{,}{v})\,=\,\texttt {two\_squares} \,{n}\,{{{\textbf {{\textsf {in}}}}}}\,{n}\,{=}\,{u}^{2}\,{{+}}\,{v}^{2}) \end{array} \end{aligned}$$

5 Improving the Algorithm

There are ways to improve the basic Iterative Algorithm in Sect. 4.1, by analysing details in the ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ composition.

5.1 Iteration Cases

Since the zagier map splits into 3 cases, the ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ composition also has 3 cases when applied to a windmill (x, y, z):

(3)

We shall name them as ping, pong and pung, respectively. The conditions in Eq. (3) give three predicates, with the corresponding transforms,

Definition 13

The predicates and transforms for three cases of ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$.

Note that each transform is given by a matrix, and the matrix for flip map is evident:

$$\begin{aligned} \begin{array}{llll} {\texttt {ping}} = \begin{bmatrix}1& \quad 2& \quad 0\\ 0& \quad 1& \quad 0\\ -1& \quad -1& \quad 1\end{bmatrix}, & {\texttt {pong}} = \begin{bmatrix}-1& \quad 0& \quad 2\\ 0& \quad 0& \quad 1\\ 1& \quad 1& \quad -1\end{bmatrix}, & \\ & {\texttt {pung}} = \begin{bmatrix}1& \quad 0& \quad -2\\ 1& \quad 1& \quad -1\\ 0& \quad 0& \quad 1\end{bmatrix}, & {\texttt {flip}} = \begin{bmatrix}1& \quad 0& \quad 0\\ 0& \quad 0& \quad 1\\ 0& \quad 1& \quad 0\end{bmatrix} \end{array} \end{aligned}$$

With these definitions Eq. (3) becomes, for a triple ${{t}\,{=}\,({x}{,}{y}{,}{z})}$,

5.2 Iteration Hops

The iterative algorithm in Sect. 4.1 can be depicted by a series ping, pong and pung. Some examples are given in Fig. 6.

Fig. 6

Triples of iteration by ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$, from zagier fix to flip fix, for various tik primes (in circle, note the initial chain of ping’s)

For a tik $n = 4 k + 1$, the starting zagier fix (1, 1, k) is a ping for ${{2}\,{{<}}\,{k}}$, i.e., ${{9}\,{{<}}\,{n}}$. Note that the iteration triples begin with a chain of ping’s, with triples of the form (x, 1, z) having x increasing and z decreasing. It is easy to verify,⁵:

$$\begin{aligned} {\texttt {ping}}^{m}\ (1,1,k)= & \begin{bmatrix}1& \quad 2& \quad 0\\ 0& \quad 1& \quad 0\\ -1& \quad -1& \quad 1\end{bmatrix}^{m} \begin{bmatrix}1\\ 1\\ k\end{bmatrix} = \begin{bmatrix}1& \quad 2m& \quad 0\\ 0& \quad 1& \quad 0\\ -m& \quad -m^{2}& \quad 1\end{bmatrix} \begin{bmatrix}1\\ 1\\ k\end{bmatrix} \\= & \begin{bmatrix}1 + 2m\\ 1\\ k - m - m^{2}\end{bmatrix} = \begin{bmatrix}x\\ 1\\ z\end{bmatrix} \end{aligned}$$

The resulting triple remains a proper windmill with ${{0}\,{{<}}\,{z}}$ as long as ${{m}\,{{<}}\,{\frac{1}{2}({{1}\,{{+}}\,\texttt {{sqrt}}\,{n}})}}$, as illustrated in Table 4. Thus if the iteration is just a chain of ping’s, the number of ping’s can be predicted, and several iterations can be hopped over in one go.

From definitions, it is not hard to prove that a ping on a proper windmill will result in an increase of its mind:

$$\begin{aligned}\begin{array}{l} \,\,\vdash {0}\,{{<}}\,{x}\,{\wedge }\,{0}\,{{<}}\,{y}\,{\Rightarrow }\,\texttt {mind}\,({x}{,}{y}{,}{z})\,{{<}}\,\texttt {mind}\,(\texttt {ping}\,({x}{,}{y}{,}{z})) \end{array} \end{aligned}$$

Therefore a chain of ping’s keeps expanding the mind of a windmill, until it becomes improper, i.e., x or y is zero or negative.

Table 4 Comparing counts of initial ping’s for various tik primes n with formula in last column

Once a non-ping triple (x, y, z) is reached, this “hopping” gets stuck, because the ping-condition ${{x}\,{{<}}\,{z}\,{{-}}\,{y}}$ is false. But then ${{{-}}{x}\,{{<}}\,{{-}}({z}\,{{-}}\,{y})\,{\iff }\,{{-}}{x}\,{{<}}\,{y}\,{{-}}\,{z}}$ is true.⁶ This is the ping-condition for $(-x,z,y)$, an inverted form of (x, y, z) by a matrix J:

$$\begin{aligned} J = \begin{bmatrix}-1& \quad 0& \quad 0\\ 0& \quad 0& \quad 1\\ 0& \quad 1& \quad 0\end{bmatrix}\quad \text {with}\quad J^{2} = I,\text { the identity; }\quad \text {and}\quad J\ \begin{bmatrix}x\\ y\\ z\end{bmatrix} = \begin{bmatrix}-x\\ z\\ y\end{bmatrix} \end{aligned}$$

The inverted triple $(-x,z,y)$ with a negative mind is a virtual windmill, in contrast to a proper windmill with a positive mind. It can start a new chain of ping’s, covering virtual windmills at first, but may reach proper windmills eventually. This is shown in the continuation of the previous examples in Table 5. For both $n = 61$ and $n = 101$, the first inversion starts a chain of ping’s and ends at a non-ping, which happens to be the flip fix. For $n = 149$, the ping of the first inversion gives (3, 7, 5), which is a pong but not a flip fix, so another inversion is applied, followed by ping to reach the flip fix.

Table 5 Continuation of chain of ping’s for various tik primes n for non-ping with virtual windmills

Further investigation reveals (see Theorem 17 later) the following transition from virtual to proper windmill,

• the last virtual windmill upon inversion gives a proper pong.

• virtual windmills before the pong, if any, are pung’s upon inversion, and

• proper windmills after the pong and before the next non-ping, if any, are ping’s.

This suggests an improved algorithm to find the two squares of a tik prime n: hop by ping for the maximum allowed to arrive at a non-ping. If a flip fix is not reached, invert the non-ping and start hopping by ping again to the next non-ping. Repeat inversion and ping-hopping until a flip fix is reached.

Such an algorithm will involve negative numbers⁷ for triples due to inversion by matrix J. In HOL4 the type num is used for non-negative numbers, while the type int handles both positive and negative numbers. The theories so far, including those for involutions, composition permutations, orbits and periods, up to the basic algorithm, are developed in HOL4 using the type num, as no negative numbers arise. Rather than recasting all definitions and theorems from num to int, we take an approach to avoid using type int altogether.

5.3 Hopping Algorithm

Recall from Sect. 5.1 that ping, pong, pung have matrices, and Sect. 5.2 that inversion is given by the matrix J. Using matrix multiplication, we can verify directly,

$$\begin{aligned} \begin{array}{rl} {\texttt {ping}}\ J & \begin{bmatrix}1& \quad 2& \quad 0\\ 0& \quad 1& \quad 0\\ -1& \quad -1& \quad 1\end{bmatrix} \begin{bmatrix}-1& \quad 0& \quad 0\\ 0& \quad 0& \quad 1\\ 0& \quad 1& \quad 0\end{bmatrix} = \begin{bmatrix}-1& \quad 0& \quad 2\\ 0& \quad 0& \quad 1\\ 1& \quad 1& \quad -1\end{bmatrix} = {\texttt {pong}}\\[15pt] J\ {\texttt {ping}}\ J & \begin{bmatrix}-1& \quad 0& \quad 0\\ 0& \quad 0& \quad 1\\ 0& \quad 1& \quad 0\end{bmatrix} \begin{bmatrix}1& \quad 2& \quad 0\\ 0& \quad 1& \quad 0\\ -1& \quad -1& \quad 1\end{bmatrix} \begin{bmatrix}-1& \quad 0& \quad 0\\ 0& \quad 0& \quad 1\\ 0& \quad 1& \quad 0\end{bmatrix} = \begin{bmatrix}1& \quad 0& \quad -2\\ 1& \quad 1& \quad -1\\ 0& \quad 0& \quad 1\end{bmatrix} = {\texttt {pung}}\\ \end{array} \end{aligned}$$

Therefore, a chain of ping’s, pong’s and pung’s can be simplified due to $J^{2} = I$, e.g.,

$$\begin{aligned} {\texttt {ping}\,{\circ }\,\texttt {ping}\,{\circ }\,\texttt {pong}\,{\circ }\,\texttt {pung}\,{\circ }\,\texttt {pung}} & \ \texttt {=}\ ({\texttt {ping}})^{2}({\texttt {ping}}\ J)(J\ {\texttt {ping}}\ J)(J\ {\texttt {ping}}\ J) \nonumber \\ & \ \texttt {=}\ {\texttt {ping}}^{5}\ J \end{aligned}$$

(4)

Indeed, when a number of ping’s is applied to an inverted triple,

We shall define hop based on this expression when $x < 2\,m z$, which is a proper windmill. When $2\,m z \le x$, the virtual windmill is inverted back into a proper windmill. Thus,

Definition 14

A hop from a triple (x, y, z) with m steps giving a proper windmill.

To predict the hop step m, we first determine the possible values for m:

Theorem 15

[script]

The range of step m giving proper windmills of n for hopping.

This means ${\texttt {hop}\,{m}\,({x}{,}{y}{,}{z})}$ will remain a proper windmill when ${{m}\,{=}\,({x}\,{{+}}\,{c})\,\texttt {{div}}\,{2}{z}}$ is within range $0 < c \le {\texttt {{sqrt}}\,{n}}$. Such considerations inspire the next definition,

Definition 16

The number of steps based on a triple (x, y, z) and a parameter c.

Hopping through nodes will reach different node types:

Theorem 17

[script]

The node types of hop m from a non-${\texttt {ping} }$ triple t.

These results match the transition from virtual to proper windmills mentioned after Table 5. Next we define how to hop with steps:

Definition 18

Hopping from a triple t with parameter c for step.

Note that hop takes a non-ping triple, to be inverted by J and then hop by ping’s. When hopping for tik $n = 4 k + 1$, instead of starting from the zagier fix (1, 1, k) which is usually a ping, we can make use of its flip (1, k, 1) which is a pong, not a ping. This simplifies the overall hopping algorithm:

It remains to prove that this hopping algorithm indeed gives the desired two squares.

5.4 Hopping Nodes

First, we define two types of skipping:

Definition 19

Skipping all ping or pung triples by WHILE loop.

Applying the idea presented in Eq. (4), it is possible to establish:

Theorem 20

[script]

Hopping with the parameter $(\texttt {sqrt} \,{n})$ is equivalent to a well-defined sequence of skipping triples from a non-${\texttt {ping} }$ triple t.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {tik} \,{n}\,{\wedge }\,{\lnot }\texttt {square} \,{n}\,{\wedge }\,{n}\,{=}\,\texttt {windmill} \,{t}\,{\wedge }\,{\lnot }\texttt {is\_ping} \,{t}\,{\Rightarrow }\\ \,\,\,\,\,\,\,\texttt {hopping} \,(\texttt {{sqrt}} \,{n})\,{t}\,{=}\,(\texttt {skip\_ping} \,{\circ }\,\texttt {pong} \,{\circ }\,\texttt {skip\_pung} )\,{t} \end{array} \end{aligned}$$

Note that skip_ping just returns the starting triple if it is not a ping; likewise for skip_pung. Hence ${\texttt {hopping}\,(\texttt {{sqrt}}\,{n})\,{t}}$ may have no pung or ping to skip, but it always skips over a single pong.

This suggests forming a list of ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ iteration windmill triples, called nodes. For a tik $n = 4k + 1$, we start with the non-ping (1, k, 1), then iterate for h times, to form a list ${\texttt {path}\,{n}\,{h}}$:

Definition 21

A path from the initial flip of zagier fix of a number n, follow by successive iterations.

From a node list ls we can select only the pong nodes:

Definition 22

Picking up only pong’s in a list of node iterations.

For each selected node, we apply pong to skip it, then skip any following ping’s to arrive at the end node of ${\texttt {hopping}\,(\texttt {{sqrt}}\,{n})}$, forming another list:

Definition 23

Manipulation of pong_nodes to hop_nodes.

We take ${\texttt {FRONT}\,{ls}}$ because when ${{ls}\,{=}\,\texttt {path}\,{n}\,{h}}$ we shall discard the last flip fix, which can be a pong. Then we prove the following:

Theorem 24

[script]

If iteration $({\texttt {path} \,{n}\,{h}})$ ends at a ${\texttt {flip} }$ fix triple, this is also the triple for the last of ${\texttt {hop\_nodes} }$.

Here $|ns |$ denotes the length of list ns. This is the critical result to show the correctness of the Hopping Algorithm.

5.5 Two Squares by Hopping

The implementation of the Hopping algorithm of Sect. 5.3 in HOL4 is,

Definition 25

The improved hopping algorithm intended for a tik prime n.

We can now state and prove the main result about this improved algorithm.

Theorem 26

[script]

The Hopping Algorithm of Definition 25 gives a ${\texttt {flip} }$ fix for a tik prime n.

$$\begin{aligned}\begin{array}{l} \,\,\vdash \texttt {tik} \,{n}\,{\wedge }\,\texttt {prime} \,{n}\,{\Rightarrow }\,\texttt {two\_sq\_hop} \,{n}\,{\in }\,\texttt {fixes} \,\texttt {flip} \,(\texttt {mills} \,{n}) \end{array} \end{aligned}$$

Proof

Let tik prime ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$, ${{p}\,{=}\,\texttt {{period}}\,(\texttt {zagier}\,{\circ }\,\texttt {flip})\,({1}{,}{1}{,}{k})}$, $h = 1 + \lfloor \frac{p}{2}\rfloor$ and ${{ls}\,{=}\,\texttt {path}\,{n}\,{h}}$. The path starts with (1, k, 1), goes to the zagier fix (1, 1, k), and continues for h iterations of ${\texttt {zagier}\,{\circ }\,\texttt {flip}}$. Due to ${{\texttt {odd}}\,{p}}$ for the iteration period of a tik prime, the path ends with a flip fix at ${\texttt {LAST}\,{ls}}$ by Theorem 10. That is, ${\texttt {flip}\,(\texttt {LAST}\,{ls})\,{=}\,\texttt {LAST}\,{ls}}$, or ${\texttt {LAST}\,{ls}\,{\in }\,\texttt {fixes}\,\texttt {flip}\,(\texttt {mills}\,{n})}$.

When ${{n}\,{=}\,{5}}, {{k}\,{=}\,{1}}$, and (1, k, 1) is already a flip fix. The algorithm immediately returns this triple. For ${{n}\,{\ne }\,{5}}$, ${\texttt {path}\,{n}\,{h}}$ with $0 < h$ has all iteration nodes distinct, i.e., ${\texttt {ALL\_DISTINCT}\,{ls}}$. Then ${\texttt {pong\_nodes}\,(\texttt {FRONT}\,{ls})}$ has all pong nodes distinct, and also ${{ns}\,{=}\,\texttt {hop\_nodes}\,{ls}}$ has all hopping nodes distinct. By Theorem 24, ${\texttt {LAST}\,{ns}\,{=}\,\texttt {LAST}\,{ls}}$ and ${\texttt {LAST}\,{ls}\,{=}\,(\texttt {hopping}\,(\texttt {{sqrt}}\,{n}))^{|{ns}|}(\texttt {HD}\,{ls})}$.

In general, for $j < |ns |$, the j-th element of $ns \,\texttt {=}\,({\texttt {hopping}\,(\texttt {{sqrt}}\,{n})})^{(j+1)} ({\texttt {HD}\,{ls}})$, which is ${\texttt {LAST}\,{ns}}$ when $j + 1 = |ns |$. Since ${\texttt {zagier}}\, \textrm{fix} {({1}{,}{1}{,}{k})}$ is in ls, it leads to a flip fix only at half of the iteration period by Theorem 10. Thus none of the hopping nodes in ns is a flip fix except ${\texttt {LAST}\,{ns}}$.

Therefore, by the property of WHILE in HOL4,

$$\begin{aligned} \texttt {WHILE}\,(\lambda \,{t}.\,\texttt {flip}\,{t}\,{\ne }\,{t})\,(\texttt {hopping}\,(\texttt {{sqrt}}\,{n}))\,(\texttt {HD}\,{ls})\,{=} \,\,\texttt {LAST}\,{ls} \end{aligned}$$

showing that the Hopping Algorithm correctly gives a flip fix for any tik prime n. $\square$

The iterations for tik primes in Fig. 6 are improved by hopping, as shown in Fig. 7. These examples show that the improvement by hopping is most significant at the start.

Fig. 7

Triples of hopping from flip of zagier fix to flip fix, for various tik primes (in circle, note that each hop contains a pong)

Implementation

Recall that hop in Definition 14 has two cases, depending on whether the triple (x, y, z) satisfies the condition ${{x}\,{{<}}\,{2}{m}{z}}$. For ${\texttt {hopping}\,(\texttt {{sqrt}}\,{n})}$ in the algorithm of Sect. 5.3, this condition is always true, so only the first case applies. Naming the first case of hop as pop, this allows direct computations in the implementation:

In HOL4, an interactive session provides EVAL for evaluation and prefix time to obtain timing statistics:

The sample HOL4 interactive dialog above shows that the 12-digit tik prime ${{123456789061}\,{=}}{\,{240081}^{2}\,{{+}}\,{256550}^{2}}$. The whole computation takes nearly 3 minutes by iteration using two_squares, but only 11 seconds by hopping using two_squares_pop, which is two_squares_hop of Definition 25 but replacing hop by pop, as shown in the pseudo-code. Note that EVAL executions are based on optimised symbolic rewriting in HOL4, thus orders of magnitude slower than running native code.

A comparison between the hopping algorithm in Sect. 5.3 and the iterative algorithm in Sect. 4.1 is shown in Table 6:

• First column shows the tik prime ${{n}\,{=}\,{4}{k}\,{{+}}\,{1}}$,

• Second column shows hopping from start to end, with hopping steps in between,

• Third column shows the number of hops,

• Fourth column shows the total number of nodes, which is $1 + \left\lfloor \frac{p}{2}\right\rfloor$ where p is the iteration period of${\texttt {zagier}\,{\circ }\,\texttt {flip}}$ for (1, 1, k), and $\left\lfloor \frac{p}{2}\right\rfloor$ is the number of steps for the iterative algorithm.

• Last column shows the timing, in milliseconds (ms), when running in HOL4 interactive mode on a typical laptop.

For this algorithm, the initial hopping step is $\left\lfloor \frac{{{1}\,{{+}}\,\texttt {{sqrt}}\,{n}}}{2}\right\rfloor$, which is large for large n. When the hopping step drops to 1, this is just the same as iteration to the next node. A typical example is the tik prime $n = 2221$: only the initial hop step 24 is large, other hop steps are small, with many one-step hops.

Table 6 Summary of sample runs for tik primes using the hopping algorithm (refer to text for explanation of columns)

Other algorithms

There are many methods to compute the two squares for a tik prime. For example, after his one-sentence proof, Don Zagier referred to an effective algorithm [17], based on ideas from the Hermite-Serret Algorithm [12]. When comparing performance, these alternative methods are better. Our algorithms are simple, deterministic, and straightforward to implement:

• For the iterative algorithm in Sect. 4.1, only addition and subtraction are performed, since doubling is just adding by itself.

• For the hopping algorithm in Sect. 5.3, the square-root computation is done before the while-loop, and only the four arithmetic operations are performed, since squaring is just multiplying by itself.

These elementary arithmetic operations are usually supported by arbitrary-precision libraries, allowing our algorithms to deal with tik primes with multi-digits.

6 Conclusion

In 1832 Gauss introduced the so-called Gaussian integers:

$$\begin{aligned} \mathbb {Z}[i] = \{x + iy \mid x \in \mathbb {Z}, y \in \mathbb {Z}\} \qquad \text {where }\mathbb {Z}\, \text {is the set of integers and }i^{2} = -1. \end{aligned}$$

A prime in $\mathbb {Z}$ may or may not remain as a prime in $\mathbb {Z}[i]$. By Fermat’s theorem, a tik prime $n = u^{2} + v^{2} = (u + iv)(u - iv)$ has a factorisation in Gaussian integers, so it is not a Gaussian prime. Thus Fermat’s two square theorem leads to the study of primes in an algebraic number field; in particular, the question whether the structure possesses unique factorisation. In view of its impact, a formal proof of Fermat’s claim is a worthy exercise in theorem-proving.

Our visual proof of Fermat’s Theorem 1 involves 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, taking up elegant ideas due to Heath-Brown [8], Zagier [18] and Spivak [16]. By developing a theory of involution-composition iteration, an algorithm to compute the two squares of the theorem can be formally shown to be correct. Details of the iterations reveal features to improve the algorithm by hopping over nodes, which uses only simple arithmetic and is straight-forward to implement.