It is not very difficult to 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:** A sketch of a proof of the claim that for every odd prime $p$ there are $(p+1)/2$ solutions $(a,b,c,d)$
in $\mathbb N^4$
(with $p=ab+cd$ and $\min(a,b)>\max(c,d)$ ) goes as follows:

  Given a solution $(a,b,c,d)$ we consider $u=(a,c),\ v=(-d,b)$.
  We associate now to $(a,b,c,d)$ the sublattice
  $\Lambda=\mathbb Zu+\mathbb Zv$ of index $p$ in $\mathbb Z^2$.
  Suppose now $cd>0$. The line $\mathbb R v$ has negative slope $<-1$
  and the two vectors $u,u-v$ belong respectively to the upper and lower
  halfplane separated by $y=0$. They belong moreover to the interior
  of the cone $\mathcal C$ defined by $\vert y\vert\leq x$. More precisely
  they are the two elements of $\Lambda$ separated by $y=0$ which
  are consecutive on the boundary $\partial H$ of $H$ where $H$
  is the convex hull
  of the set $\mathcal C\cap(\Lambda\setminus\{0\})$.

  Sublattices of index $p$ in $\mathbb Z^2$ are in one-to-one
  correspondence with points of the projective line over $\mathbb F_p$:
  Given a point $[s:t]$ (where we chose $s,t\in \mathbb Z$)
  representing an element of $\mathbb P^1\mathbb F_p$, we consider the lattice
  $$\Lambda_{[s,t]}=\{(x,y)\in\mathbb Z^2\ \vert\ xs+yt\equiv 0\pmod p\}\ .$$
  The previous discussion shows that every sublattice of index $p$ in
  $\mathbb Z^2$ corresponds to at most a unique solution $(a,b,c,d)$.
  The points $[1:0]$ and $[0,1]$ correspond to $(p,1,0,0)$ and $(1,p,0,0)$.
  We leave it to the reader to check that $[1:1]$ and $[1:-1]$ do not correspond to a solution.
  For $[1:\mu]$ and $[1,-\mu]$ with $\mu\in\{2,\ldots,p-2\}$,
  the two associated lattices are symmetrical
  with respect to the vertical axis and exactly one of them
  corresponds to a convex hull $H$ as above whose side crossing $y=0$
  has a negative slope (since this side is compact and contained in the boundary $\partial H$ of the convex $H$ containing
almost all of $\mathcal C$, its slope is $<-1$).

  We get thus a total number of $2+(p-2-2+1)/2=(p+1)/2$ solutions.