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

Second Addition: The "sketch of proof" below is a bit cheesy (it contains a hole): The vector $(u-v)$ can lie outside of the cone $\mathcal C$ (Example: $19=3\cdot 6+1\cdot 1$ gives rise to $u=(3,1)$, $v=(-1,6)$ and $u-v=(4,-5)$ outside $\mathcal C$). I will think it over calmly and hopefully patch it in a few days. (end of second addition).

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.