One can show (see below for a sketch of a proof) that every odd prime number $p$
can be written in exactly $(p+1)/2$ different ways as
$$p=a\cdot b+c\cdot d$$
with $a,b,c,d\in\mathbb N$ satisfying $\max(c,d)<\min(a,b)$.
Example for $p=23$:
\begin{matrix}
1\cdot 23+0\cdot 0 & 23\cdot 1+0\cdot 0 \\
2\cdot 11+1\cdot 1 & 11\cdot 2+1\cdot 1 \\
3\cdot 7+1\cdot 2 & 3\cdot 7+2\cdot 1 \\
7\cdot 3+1\cdot 2 & 7\cdot 3+2\cdot 1 \\
4\cdot 5+1\cdot 3 & 4\cdot 5+3\cdot 1 \\
5\cdot 4+1\cdot 3 & 5\cdot 4+3\cdot 1
\end{matrix}
Klein's Vierergruppe $\mathbb V$ acts on all such quadruplets $(a,b,c,d)$
by permuting the first two, permuting the last two, or permuting both the first two and the last two elements. So we get an easy proof that every prime
$p$ congruent to $1$ modulo $4$ must be a sum of squares:
$(p+1)/2$ is then odd and the only fixed points under the action
of $\mathbb V$ are of the form $(a,a,c,c)$.
*Does somebody know a reference for this proof?* It looks a bit like
Zagier's proof which also uses a parity argument for a set acted upon by involutions.
**Motivation:** This is in fact a variation of
https://mathoverflow.net/questions/405035 .
**Added (16th of october 2022):** The final version, called 'A Quixotic Proof of Fermat's Two Squares Theorem for Prime Numbers', has been accepted by the American Math. Monthly. Up to formatting changes, labellings of results and minor differences it coincides with the preprint
https://hal.archives-ouvertes.fr/hal-03813904 .
**Added (first of november 2021):** A detailled proof is contained in
https://hal.archives-ouvertes.fr/hal-03408135/document (also available from the arXiv in a few days).
**Sketch of proof**
Given a solution $(a,b,c,d)$ we consider $u=(a,c),\ v=(-d,b)$.
We associate to $(a,b,c,d)$ the sublattice
$\Lambda=\mathbb Zu+\mathbb Zv$ of index $p$ in $\mathbb Z^2$.
Suppose now $cd>0$ and consider the eight open cones
of $\mathbb R^2$ defined by the complement of the four lines defined by
$xy(x^2-y^2)=0$. We colour these open cones alternatingly black
and white. The four vectors $\pm u,\pm v$
are contained in different black cones (colouring the first cone
above the halfline $(\mathbb R_{>0},0)$ in black).
We say that a sublattice $\Lambda$ of finite index in
$\mathbb Z^2$ has a monochromatic basis
if there exists a basis $b_1,b_2$ of
$\Lambda=\mathbb Z b_1+\mathbb Z b_2$ such that all four elements
$\pm b_1,\pm b_2$ belong to different open cones of the same colour.
(Not every lattice has a monochromatic basis but many do.)
We claim that all monochromatic bases of a lattice (having a monochromatic basis) are
of the same colour: If $b_1,b_2$ is a black monochromatic basis
and $w_1,w_2$ is a white monochromatic basis, then $w_1,w_2$
belong to two open adjacent cones of $\mathbb R^2\setminus(\mathbb Rb_1\cup
\mathbb R b_2)$ which is impossible by the following small Lemma:
**Lemma:** If $f_1,f_2$ and $g_1,g_2$ are two bases of a
lattice $\Lambda=\mathbb Z f_1+\mathbb Z f_2=\mathbb Z g_1+\mathbb Z g_2$
such that $\{\pm f_1,\pm f_2\}$ and $\{\pm g_1,\pm g_2\}$ do not intersect,
then $g_1,g_2$ or $g_1,-g_2$ are contained in a common connected component
of $\mathbb R^2\setminus(\mathbb R f_1\cup \mathbb R f_2)$. (Otherwise
we have up to sign changes and exchanges of indices $f_1=\alpha b_1+\beta b_2$
and $f_2=\gamma b_1-\delta b_2$ with $\alpha,\beta,\gamma,\delta$ strictly positive integers. This implies that $b_1$ belongs to the convex hull of
$(0,0),f_1,f_2$ which is a contradiction.)
Set now $\Lambda_\mu=\{(x,y)\in\mathbb Z,\ \vert\ x+\mu y\equiv 0\pmod p\}$.
If $\mu\in \{2,\ldots,p-2\}$, then $\Lambda_\mu$ contains no
elements of the form $(\pm m,0),(0,\pm m),(\pm m,\pm m)$ with $m$
in $\{1,\ldots,p-1\}$. This implies that every open black or white
cone contains a point with coordinates of absolute value at most $p-1$.
A reduction algorithm implies the existence of a monochromatic basis.
(Start with two arbitrary non-zero elements $e_1,e_2$ of $\Lambda$ having coordinates of
absolute value smaller than $p$ which belong to two different consecutive black cones. If the interior of the convex hull spanned by $\pm e_1,\pm e_2$
contains a non-zero element in a black cone, we can replace $e_1$ or $e_2$
and decrease the area of the convex hull spanned by $\pm e_1,\pm e_2$.
If the interior contains no non-zero elements of $\Lambda$ in black
cones, we get either a monochromatic basis or the convex hull contains
at least four lattice points in four distinct white cones and we switch the working colour to white.)
Moreover, since $\Lambda_\mu$ and $\Lambda_{p-\mu}$ differ by a horizontal
reflection, they have monochromatic bases of different colours.
Retaining only lattices with black monochromatic bases,
We get $(p-2-2+1)/2=(p-3)/2$ such lattices with black monochromatic
bases.
Monochromatic bases of a lattice $\Lambda_\mu$
are not unique but in finite number. It remains
to show that exactly one black monochromatic basis of a lattice $\Lambda_\mu$
has the form $u=(a,c),v=(-d,b)$ as required for a solution of
$p=ab+cd$ (with $\min(a,b)>\max(c,d)$ and $0\leq c,d$). We call
such a basis a reduced monochromatic basis.
First observe that every lattice with a black monochromatic basis
$b_1,b_2$ (where we suppose $b_1\in\mathbb N^2$ and
$b_2$ in $(-\mathbb N)\times\mathbb N$)
has a reduced black monochromatic basis: Replace $b_1$ by $b_1-kb_2$
with $k$ maximal for monochromaticity. Replace then $b_2$ by
$b_2+sb_1$ with $s$ maximal for monochromaticity. The resulting
black monochromatic basis is reduced.
Observing that the two lattices $\mathbb Z(p,0)+\mathbb Z(1,\pm 1)$
contain no vectors $u,v$ (associated to a solution $(a,b,c,d)$ such that ...)
and adding the two trivial solutions $(p,1,0,0),(1,p,0,0)$
(corresponding to the lattices $\mathbb Z(p,0)+\mathbb Z(0,1)$
and $\mathbb Z(1,0)+\mathbb Z(0,p)$) we get a total number of at
least $(p-3)/2+2=(p-1)/2$ solutions and we are done after showing
that every lattice with a black monochromatic basis contains only one
reduced monochromatic basis (also using the fact that sublattices
of prime index $p$ are in one-to-one correspondence with points
of the projective line over $\mathbb F_p$).
Supose now that $u=(a,c),v=(-d,b)$ is a reduced black basis
and let $u'=(a',c'),v'=(-d',b')$ be a second reduced black basis
giving rise to two distinct solutions $(a,b,c,d)$ and $(a',b',c',d')$.
Since $u$ and $v$ determine each other uniquely in a reduced black
basis, we can assume that $u'\not=u$ and $v'\not=v$.
The two vectors $u',v'$ are thus contained in the four open
cones defined by $\mathbb R^2\setminus(\mathbb R u\cup\mathbb R v)$.
The lemma used previously
shows that they can not belong to two adjacent cones of
$\mathbb R^2\setminus(\mathbb R u\cup\mathbb R v)$.
We suppose now that $u',v'$ belong to $\mathcal C\cup (-\mathcal C)$ for
$\mathcal C$ a cone (onnected component) of $\mathbb R^2\setminus(\mathbb R u\cup
\mathbb R v)$. If $u'$ and $v'$ belong to two opposite cones, we exchange
the basis $u,v$ with the basis $u',v'$. We can now assume that both vectors
$u'$ and $v'$ belong to the open cone $(0,+\infty)u+
(0,+\infty)v$ spanned by $u$ and $v$. We have thus $u'=\alpha u+\beta v$ and $v'=\gamma u+\delta v$
with $\alpha,\beta,\gamma,\delta$ strictly positive integers. Reducedness
of the black monochromatic basis $u,v$ implies that $v+u$ is either white
or belongs to the black cone $\mathcal C_u$ containing $u$. If $u+v$ is white, we get a contradiction
by observing that it is contained in the convex hull of $(0,0),u',v'$
(which contains no other points of $\Lambda$). The point $u+v$ is
thus in the black cone $\mathcal C_u$ containing $u$.
Geometric considerations imply now $\delta\geq 3$ and the impossible
inequalities
$$p\geq (3 b+c)a/2>(3ab+ac)/2>ab+a(b+c)/2>ab+cd=p\ .$$
Indeed, since $u$ and $u+v$ belong both to $\mathcal C_u$,
the line $u+\mathbb Rv$ of slope $<-1$ intersects the
white cone separating $\mathcal C_u$ from the black cone $\mathcal C_v$
containing $v$ in a segment containing at least one lattice-point of $\Lambda$. This implies $\delta\geq 3$ and the second coordinate of $v'$ is
thus at least equal to $3b+c$. On the other hand, the first coordinate
of $u'$ has to be at least equal to the first coordinate of the
intersection $u+\mathbb Rv \cap \mathbb R(1,1)$ which is $\geq a/2$
(since $v$ has slope $<-1$).