To calculate the fixed points of $g$ we'll have to solve the equation $$ g(c)=c $$ Since $g(S_1)=S_3$,$g(S_3)=S_1$ and $g(S_2)=S_2$ the only fixed points of $g$ must lie in $S_2$. $$ g(x,y,z) = (2y-x,y,x-y+z)= (x,y,z) $$ holds if $y=x$, but at the same time $p=x^2+4yz=x(x+4z)$ and so if $p$ is prime the only option is $y=x=1$ and $p=4z+1$. Note that if $x\neq 1$, $p$ would be a product of 2 numbers thus not a prime. There's also a really well done video by the Numberphile team about this paper [](https://www.youtube.com/watch?v=SyJlRUBoVp0) \(c\) is called a "fixed point" of \(f\) because if you start a recursive sequence: \[a_{n+1}=f(a_{n})\] \[a_{n+2}=f(a_{n+1})\] and it reaches \(c\) the sequence will stay fixed at \[a_n = c = f(c)\] To sum up: Zagier starts by creating a function $g$ which is an involution in $S$. He finds that $g$ has a fixed point which implies that $p=4z+1$, $z\in \mathbb{N}$. He then uses the lemma about the cardinalities to prove that the other involution $f(x,y,z)=(x,z,y)$ has at least one fixed point. Calculating the fixed point of $f$ we finally get to the result $p=x^2+(2z)^2$, $z,x\in \mathbb{N}$. In other words $p$ is a sum of 2 squares! I think you mean $g(S_1) = S_3, g(S_3) =S_1$, right? Lemma: The cardinalities of a finite set and of its fixed-point set under any involution have the same parity. Proof: Let $S=\{s_1,...,s_n\}$ be the finite set in question. If $s_i$ is not a fixed point then there exists an $s_j$ such that $f(s_i)=s_j$. Note that this forces $f(s_j)=s_i$, so if we delete $s_i$ and $s_j$ from $S$ we obtain a new set $S'$ such that f is an involution in this set, and we have $|S|\equiv |S'|\mod{2}$, in such a manner delete all the non-fixed points, so the result follows. Using the lemma above and the fact that g has a unique fixed point (so in particular has odd parity), implies that $|S|$ has odd parity, so the number of fixed points of $f$ has odd parity, which means that there is at least one fixed point of $f$. For $f$, the fixed point equation $$ f(x,y,z) = (x,z,y)= (x,y,z) $$ holds if $y=z$ and so $p=x^2+(2z)^2$. We finally proved that $p$ is a sum of 2 squares! In this paper Don Zagier constructs a really short proof of the very famous Fermat's theorem on sums of two squares. Fermat's theorem on sums of two squares says that an odd prime $p$ is expressible as $$p = x^2 + y^2$$ with $x$ and $y$ integers, if and only if $$ p \equiv 1 \pmod{4} $$ which is the same thing as saying if and only if $p=4k+1$, where $k \in \mathbb{N}$. The prime numbers for which this is true are called Pythagorean primes. For example, the primes 5, 13, 17, 29, 37 and 41 are all congruent to 1 modulo 4, and they can be expressed as sums of two squares in the following ways: $$ 5 = 1^2 + 2^2 \\ \quad 13 = 2^2 + 3^2\\ \quad 17 = 1^2 + 4^2 $$ Fermat himself claimed to have proven the theorem in a letter to Marin Mersenne in 1640, but didn't provide any proof. The theorem was later solved by several mathematicians including Euler, Gauss, Heath-Brown and Don Zagier. First we should recall what is an involution. An involution is a function $f$ that is its own inverse: $$ f(f(x))=x $$ for all x in the domain of $f$. A fixed point of a function is an element of the function's domain that is mapped to itself by the function $$ f(c)=c $$ From a visual perspective a fixed point is the point where the function intersects the line $y=x$. Here is an example of a function with 3 fixed points  Don Zagier, born in 1951, is an American mathematician whose main area of work is number theory. He was a child prodigy and got a bachelor's and master's degree from MIT at the age of 16 and a PhD from the University of Bonn at the age of 20.  Zagier starts by defining two involutions $$f(x,y,z) = (x,z,y)$$ $$ g(x,y,z)=\begin{cases} (x+2z,~z,~y-x-z),\quad \textrm{if}\,\,\, x < y-z \\ (2y-x,~y,~x-y+z),\quad \textrm{if}\,\,\, y-z < x < 2y\\ (x-2y,~x-y+z,~y),\quad \textrm{if}\,\,\, x > 2y \end{cases} $$ It's easy to see that $f$ is an involution because $$ f(f(x,y,z))=f(x,z,y)=(x,y,z) $$ It's a little more cumbersome to prove that $g(x,y,z)$ is an involution. To do that we first should note that $g$ divides $S$ in 3 regions: $$S_1 = \{(x,y,z) ∈ S : x < y−z \}$$ $$S_2 = \{ (x,y,z) ∈ S : y−z < x < 2y \}$$ $$S_3 = \{ (x,y,z) ∈ S : x > 2y \}$$ such that $S=S_1 \cup S_2 \cup S_3 $. An interesting property about $g$ is that $g(S_1)=S_3$, $g(S_1)=S_3$, $g(S_2)=S_2$. In other words the solutions of $S_1$ are mapped to $S_3$ and vice-versa and the solutions of $S_2$ are mapped to $S_2$  Let's now prove this for every region: For $S_1$ $$ g(x,y,z)=(x+2z,z,y-x-z) $$ To check that $g(S_1)=S_3$ we need to see if the condition $x>2y$ is satisfied. In this case $x>2y \rightarrow x+2z>2z \equiv x>0 $ which is always true since $x \in \mathbb{N}$. We do the same for $S_3$ $$ g(x,y,z)=(x-2y,x-y+z,y) $$ To check that $g(S_3)=S_1$ we need to see if the condition $x<y-z$ is satisfied. In this case $x<y-z \rightarrow x-2y<x-y+z-y \equiv z>0 $ which is always true since $z \in \mathbb{N}$. Finally for $S_2$ $$ g(x,y,z)=(2y-x,y,x-y+z) $$ To check that $g(S_2)=S_2$ we need to see if the conditions $y-z<x$ and $x<2y$ are satisfied. In this case $y-z<x \rightarrow y-x+y-z<2y-x \equiv z>0 $ and $x<2y \rightarrow 2y-x<2y \equiv x>0 $ which are both true since $z,x \in \mathbb{N}$. It's now easy to see that $g$ is an involution since: $$ g(g(S_1))=g(S_3)=S_1 $$ $$ g(g(S_3))=g(S_1)=S_3 $$ $$ g(g(S_2))=g(S_2)=S_2 $$ We should also verify that $g$ maps $S$ into itself, in other words that the image will also satisfy $x^2+4xy = p$. I'm just going to verify that for $S_1$ (for $S_2$ and $S_3$ is exactly the same procedure) $$ (x+2z)^2+4y(y-x-z)=x^2+4z^2+4xz+4zy-4xz-4z^2=\\ =x^2+4zy \in S $$