Explanation:
first of all a rather trivial point: a square is a sum of two squares (of which usually one will need to be 0²)
now you need a lemma: if a is a sum of two squares and b is a sum of two squares, then the product ab is also a sum of two squares.
this is easy to see in terms of the multiplication of two Gaussian integers:
(a²+b²)(c²+d²) = |a+ib|² |c+id|² =|(a+ib)(c+id)|² = |(ac−bd) + i(ad+bc)|² = (ac−bd)² + (ad+bc)²
since 2=1²+1² the lemma implies any power of two is a sum of two squares.
a well-known theorem tells us that any prime of the form 4n+1 is the sum of two squares, so the lemma implies that any product of such primes p1^(k1) p2^(k2)⋯pm^ (km) is a sum of two squares.
finally to multiply a sum of two squares by a square (of the form q2j11q2j22⋯q2jnn) is again a sum of two squares.
the proof of the theorem about primes of the form 4n+1 is a fairly easy consequence of two basic results: Wilson's theorem and the fact that the Gaussian integers are a Euclidean ring and therefore are a Unique Factorization Domain. Wilson's theorem is a simple exercise in modular arithmetic. To demonstrate a division algorithm for the Gaussian integers is a little more technical but of an elementary nature.