I think we can always translate everything to the basics operations but it will take a huge number of pages! Fortunately, your question has an elementary answer. We denote by $S_n$ the set of representation of the integer $n$ as sum of squares, we have $(a,b)\in S_n$ if and only if $0\leq b<a$ and $ a^2+b^2=n$;
Given two primes $p$ and $q$ let:
$$ S_p=\{(a_p,b_p)\}\\
S_q=\{(a_q,b_q)\}
$$
Now we want to prove that $S_{pq}$ contains exactly two elements, or in other words we want to prove that the following equations on unknown $(a,b),(c,d)$ has exactly one solution :
- $a>b>0$ and $c>d> 0$ (it's clear that neither $d=0$ or $b=0$ occurs)
- $a^2+b^2=c^2+d^2=pq$
- $a>c$ ( different elements)
The equation in the middle is equivalent to:
$$(a-c)(a+c)=(d-b)(d+b)$$
this implies the existence of integers $x,y,z,t$ pairwise coprime such that:
$$\begin{align} a-c&=&2xy\\ a+c&=&2zt\\d-b&=&2xz\\d+b&=&2yt \end{align}$$
Note that the existence is not hard for example $x=gcd\big(\frac{a-c}{2},\frac{d+b}{2}\big),y=gcd\big(\frac{a-c}{2},\frac{d-b}{2}\big)\cdots$, so this will justify also that they are pairwise corprime.
The result is the fact that:
$$pq=(x^2+t^2)(y^2+z^2)$$
obviously $xyzt\neq 0$ therefore $\{p,q\}=\{x^2+t^2,y^2+z^2\}$ we have two cases:
- $p=x^2+t^2$ and $q=y^2+z^2$, because $a+c>d-b$ we have $t>x$, and $a+c > d+b$ implies $z>y$ so $x=b_p,t=a_p,y=b_q, z=a_q$ so $a=a_pa_q+b_qb_p,b=a_pb_q-b_pa_q,c=a_pa_q-b_qb_p, d=a_pb_q+b_pa_q$
- $q=x^2+t^2$ and $p=y^2+z^2$ because of the symetry it gives different $x,y,z,t$ but the same $a,b,c,d$ as the first case.
So there is only one solution to the equations hence $S_{pq}$ contains exactly two elements.
Note that $a^2 + b^2 = c^2 + d^2$ is equivalent to $a^2 - c^2 = d^2 - b^2$, i.e. $(a-c)(a+c) = (d-b)(d+b)$. If we factor any odd number $m$ as $m = uv$, where $u$ and $v$ are both odd and $u < v$, we can write this as $m = (a-c)(a+c)$ where $a = (u+v)/2$ and $c = (v-u)/2$. So any odd number with more than one factorization of this type gives an example.
Thus from $m = 15 = 1 \cdot 15 = 3 \cdot 5$, we get $8^2 - 7^2 = 4^2 - 1^2$, or $1^2 + 8^2 = 4^2 + 7^2$.
From $m = 21 = 1 \cdot 21 = 3 \cdot 7$ we get $11^2 - 10^2 = 5^2 - 2^2$, or $2^2 + 11^2 = 5^2 + 10^2$.
Best Answer
Consider the equation $x^2+y^2=z^2+w^2=N.$
This is equivalent to: $x^2-z^2=w^2-y^2=D.$
and we need $D$ to have at least two different factorizations with the factors having the same parity. One approach is thus to find the smallest values of $D$ which satisfy this property. Another approach is to find parametric solutions.
Let $(x-z)=ab$ and $(x+z)=cd$ with $(w-y)=ac$ and $(w+y)=bd$,
which gives the parametric solution: $\dfrac{1}{2}(ab+cd,bd-ac,cd-ab,ac+bd)$. We have $N=\dfrac{1}{4}(a^2b^2+c^2d^2+b^2d^2+a^2c^2)$.
Now from the fact that $cd-ab$ and $bd-ac$ are distinct positive integers, we can see that no three of $a,b,c,d$ can be the same.
Thus, the lexicographically smallest value for the tuple $(a,b,c,d)$ is $(1,1,2,3)$, which gives the solution $(7,1,5,5)$ and $N=50$. Once we have this solution, it is easily checked that other tuples of $a,b,c,d$ give larger values of $N$.