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 for $a,b,c,d \in \mathbb{R}$ we have
\begin{align}
(a^2 + b^2)(c^2 + d^2) = (ac + bd)^2 + (ad - bc)^2
\end{align}
This means if two numbers can be written as the sum of squares, then their product can also be written as the sum of squares. Consider the $q_j$. Each $b_j$ is even so write $b_j = 2b'_j$. We then have $q_j^{b_j} = (q_j^{b'_j})^2 + 0^2$ so each of them can be written as the sum of sqaures. Thus $q_1^{b_1} \cdots q_l^{b_l}$ can be written as the sum of squares. Similarly be Fermat's theorem on the sum of squares each $p_i$ can be written as the sum of squares, and thus $p_i^{a_i}$ can be written as the sum of squares. Also $2 = 1^2 + 1^2$. Thus the product $m = 2p_1^{a_1} \cdots p_k^{a_k} q_1^{b_1} \cdots q_l^{b_l}$ can be written as the sum of squares.
Best Answer
Since $n$ can be written as a sum of two squares in a number of ways that depends on how many ways there are to split $n$ as: $$ n = z\bar{z} = (a-bi)(a+bi) $$ over $\mathbb{Z}[i]$, and the latter is an Euclidean domain hence a UFD, it follows that for any $n\in\mathbb{N}$:
$$\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2= n\} = 4(\chi_4*1)(n) = 4\left(d_1(n)-d_3(n)\right)$$ where $d_1(n)$ is the number of divisors of $n$ of the form $4k+1$ and $d_1(n)$ is the number of divisors of $n$ of the form $4k+3$. If $p$ and $q$ are primes of the form $4k+1$, it follows that: $$\#\{(a,b)\in\mathbb{N}^2:a^2+b^2=pq\} = \color{red}{4}.$$ These representations can be recovered from Lagrange's identity: $$ (a^2+b^2)(c^2+d^2) = (ac-bd)^2+(ac+bd)^2.$$ For instance, let we assume $p=13$ and $q=41$. Then: $$ p = 2^2 + 3^2,\qquad q=4^2+5^2$$ and $pq=13\cdot 41=533$ can be written as a sum of two squares in the following ways: $$ pq = 2^2 + 23^2 = 7^2 + 22^2$$ since $2\cdot 4+3\cdot 5=23$ and $2\cdot 5+3\cdot 4=22$.