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 (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$).