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 Arithmetic properties of positively reduced $2\times 2$-matrices .
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$).
