An impressively terse presentation, but that makes it not particularly clear or likely to spark further understanding in someone who doesn't already understand

An impressively terse presentation, but that makes it not particularly clear or likely to spark further understanding in someone who doesn't already understand the phenomenon intuitively. Here's the approach I prefer for clarity and understanding: First, recognize that prime p can be written as a sum of two squares just in case A^2 + B^2 = 0 (mod p) has a nonzero solution. In one direction (left to right), this is obvious; in the other direction (right to left), given A^2 + B^2 = kp with A and B indivisible by p, we can compute A' + B'i = gcd(A + Bi, p) in the Gaussian integers using the Euclidean algorithm. The result is a proper factor of p; however, it is not 1, as A + Bi is not co-prime to p (if it were, so would be its conjugate, and therefore so would be their product A^2 + B^2 = kp, which clearly isn't). Thus, |A' + B'i|^2 = A'^2 + B'^2 is a nontrivial proper factor of |p|^2 = p^2, which means A'^2 + B'^2 = p on the nose. Next, recognize that A^2 + B^2 = 0 (mod p) just in case (A/B)^2 = -1 (mod p). [Recall that we can do division in arithmetic modulo a prime; in jargon, such arithmetic comprises a "field"]. Thus, the question of which primes can be written as a sum of two squares reduces to the question of which primes admit a square root of -1 in their modular arithmetic. Finally, we can apply Euler's criterion to this, which tells us x is a square modulo odd prime p just in case x^((p + 1)/2) = x (mod p) [I can expound on why this is true in more detail if readers are interested; essentially, there are (p + 1)/2 many squares modulo p (one is 0, and other (p - 1)/2 come from the other values mod p paired off under negation), and Fermat's Little Theorem tells us they're all solutions to this equation]. Applying this criterion with x = -1, we see that -1 has a square root modulo odd prime p just in case p is 1 (mod 4), which, by all the above, means an odd prime is a sum of two squares just in case it is 1 (mod 4). Pushing further on the ideas above by thinking about factorization of Gaussian integers more generally gives us more: it gives us a convenient formula for exactly how many ways there are to write an arbitrary integer N as a sum of two squares. (And by carrying out the analogue with quaternions rather than complex numbers, we get a formula for four square representations as well). But, again, all this I will leave to followup discussion if desired.