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.
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$.
Best Answer
This is not the smallest solution, but it gives you a simple method to find an integer that is the sum of three positive cubes in an arbitrary number of ways. All you have to do is use Lehmer's Identity,
$$(-9 p^3q + q^4)^3 + (-9p^4 + 3p q^3)^3 + (3p^2)^6 = q^{12}$$
If $q = 1$, then you get near-misses to the eqn $a^3+b^3 = c^3$.
However, you can also choose a large enough constant $q$ such that the first two addends are positive for $p = 1,2,3,\dots,n$. For $n=13$, I find $q=28$ is enough, so,
$$\begin{aligned} 614404^3+65847^3+3^6 &= 28^{12}\\ 612640^3+131568^3+12^6 &=28^{12}\\ 607852^3+196839^3+27^6 &=28^{12}\\ 598528^3+261120^3+48^6 &=28^{12}\\ \vdots\end{aligned}$$
and so on for 13 $p$ before the first addend becomes negative. Of course, if you choose larger $q$, then the addends will be positive for an even longer range of $p$.