Promised proof (cobbled together from previous writing such as this, tailored to current situation, but with no effort on doing superscript formatting, or picking best variable names, or any such thing, because I’m lazy and have more important work to do (yet I am clearly here procrastinating)):

A Pythagorean triple is a triple of integers A, B, and C such that A^2 + B^2 = C^2. (It’ll be convenient for now to allow these to be negative; we can cull those later). We’re looking to determine how many ways we can write C^2 as a sum of two squares. In fact, even for values K which aren’t square numbers, we might ask how many ways they can be written as a sum of two squares. (Of course, the answer will be 0 for negative K, and, just as trivially, 1 if K is zero. In the following, K is always presumed positive.)

It turns out there’s another convenient way to think about all this: by a “Gaussian integer”, I mean a complex number of the form A + iB with A and B both integers.

A Pythagorean triple of hypotenuse C is a Gaussian integer of magnitude C; we might also just ask for its squared magnitude to be K, whether or not K is the square of an integer.

The benefit of looking at things this way is that Gaussian integers have a structure very similar to ordinary integers. In particular, they come with a notion of prime factorization very similar to that for ordinary integers.

To see this, let’s review how prime factorization works:

Let’s call two integers “similar” if each is a multiple of the other (so, e.g., x and -x are similar; note that similar values have the same size). From now on, except where otherwise specified, I will not distinguish between similar values; they’re as good as equal to me.

We say a nonzero integer p is “prime” if any way of expressing it as a product of integers involves a factor of p. (Note that 1 (or any value similar to 1; aka, any “unit”) is not considered prime because of the empty product.)

By repeatedly breaking any particular nonzero integer into a product of smaller integers, if possible, and stopping when reaching primes, we find that every nonzero integer has some representation as a product of primes. (Note that this factorization process cannot go on forever, by consideration of how sizes get smaller at each stage)

Not only that, but this representation is unique. (This follows from the fact that if a prime divides a product, it divides one of the product’s factors. This in turn follows from the fact that arithmetic modulo a prime forms a field (that is, if prime p does not divide x, there is some multiple of x which is 1 plus a multiple of p). This follows from the fact that any two (or, indeed, finitely many) values have a greatest common divisor which is a sum of multiples of those values (in particular, the GCD of a prime and any value it does not divide will be 1). This follows from running the extended Euclidean algorithm, which makes use of a “division and remainder (smaller than the divisor)” operator, whose existence just depends on the fact that any ratio of integers has distance less than 1 from some integer. Again, this algorithm cannot go on forever by consideration of how sizes get smaller at each stage)

All of the above was phrased in terms of integers, but, in fact, it works just as well replacing “integer” by “Gaussian integer” throughout! So we have that each nonzero Gaussian integer uniquely factorizes into a product of Gaussian primes.

In particular, each ordinary integer also has a factorization as a product of Gaussian primes, which we can obtain by first factoring it into ordinary primes, and then further factoring those ordinary primes into Gaussian primes. So let’s try to understand how ordinary primes factor into Gaussian primes.

Let p be an ordinary prime and let q be a Gaussian prime factor of p. By symmetry, q’s conjugate must also be a prime factor of p. It may be that q is similar to its conjugate or not.

[ul]

[li]If q and its conjugate aren’t similar, then their product is a non-unit integer factor of p, and thus must be p itself. In this case, let’s say p is a “bifurcating” prime; it’s prime within the integers, but splits into two dissimilar conjugate primes within the Gaussian integers.[/li][li]On the other hand, if q and its conjugate are similar, then q is either similar to an integer or similar to an integer multiple of 1 + i. In the former case, as p is prime within the integers, we must have that q is similar to p itself, and thus p is already a Gaussian prime (let’s call such p “atomic”). In the latter case (with q similar to a multiple of 1 + i), we have that p is divisible by 1 + i, and thus p^2 = |p|^2 is divisible by |1 + i|^2 = 2; this can only occur when p itself is even, which is to say, when p is 2. And, indeed, 2 has the Gaussian prime factorization (1 + i) * (1 - i), a product of two similar conjugates.[/li][/ul]

In summary, every odd ordinary prime is either itself a Gaussian prime (in which case we’ll call it “atomic”), or the product of two dissimilar conjugate Gaussian primes (in which case we’ll call it “bifurcating”). The one other ordinary prime is 2, which has the Gaussian prime factorization (1 + i) * (1 - i), consisting of two similar conjugates.

Thus, the Gaussian prime factorization of K can be determined from the ordinary prime factorization of K: simply take the latter and split each 2 into two copies of (1 + i) [up to similarity], as well as splitting each bifurcating prime into the appropriate two two dissimilar conjugate Gaussian primes, and leaving each atomic prime alone.

Now we can tackle the question of how many ways (up to similarity) there are to write K as A^2 + B^2 = |A + iB|^2 = (A + iB) * (A - iB). This is the number of ways to split the Gaussian prime factors of K into two groups, each group consisting of the conjugates of the other. Our hands are tied for those Gaussian prime factors of K arising from 2: we will have an even number of copies of (1 + i), and must simply split these half and half among our two groups ((1 + i) counting as self-conjugate up to similarity). Similarly, our hands are tied for those Gaussian prime factors of K arising from any particular atomic prime: again, we must split these half and half (if possible; i.e., if it appears an even number of times. Otherwise, we will not be able to write K as a sum of two squares at all). Finally, and most interestingly, for each bifurcating prime factor of K, we get some choice: if it appears e times in the ordinary prime factorization of K, then we get e copies of p and e copies of p’, for dissimilar conjugate Gaussian primes p and p’, in the Gaussian prime factorization of K. Partitioning these into two groups, each consisting of the conjugates of the other, can be done in e + 1 ways in total (we can choose anywhere from 0 to n copies of p to put into the first group, and everything else is uniquely determined by this).

In conclusion, the number of ways (up to similarity) to write K as A^2 + B^2 will be zero if any atomic prime appears with an odd exponent in the prime factorization of K, and otherwise, it will be the product of (e + 1) over each exponent e of a bifurcating prime in the prime factorization of K. (In the case where K = C^2, the prime factorization of K automatically has only even exponents, and the number of ways will become the product of (2e + 1) over each exponent e of a bifurcating prime in the prime factorization of C).

[We’ve counted up to similarity, but luckily, every Gaussian integer is similar to precisely one Gaussian integer of the form A + iB “in the first quadrant” (that is, with neither A nor B negative, and also with A positive if B is). So this is basically the same as counting the number of solutions with positive A and positive B, except it also allows for B to be zero if A is positive. In our case of interest (A^2 + B^2 = C^2 with positive C), there is precisely one solution of this form (A = C) which we can subtract out if we do not wish to count it.

This also counts A^2 + B^2 distinctly from B^2 + A^2; you can divide the count in half to get rid of that if you like (again, luckily, you will never have solutions where A^2 = B^2 as the exponent of 2 on the left-hand side of A^2 + B^2 = C^2 would be odd but on the right-hand side would be even).]

Phew! That’s quite a bit already, but we’re almost done. Now we just need to figure out which odd primes are atomic/bifurcating.

If ordinary prime p equals (A + iB) * (A - iB) = A^2 + B^2, then B won’t be divisible by p, and so A/B will be a well-defined square root of -1 in the field of integers mod p. Conversely, if there is a square root r of -1 in the field of integers mod p, then we have that 1 + r^2 = (1 + ir) * (1 - ir) is a multiple of p, even though p divides neither factor, and thus, p cannot be a Gaussian prime itself. Thus, p is atomic just in case there is no square root of -1 modulo p.

How do we tell if -1 has a square root modulo p? More generally, how do we tell if a value has a square root modulo p? Well, note that two values have the same square just in case they are each other’s negations; restricting attention to the case where p is odd, no nonzero value is its own negation, so the p - 1 nonzero values square to precisely (p - 1)/2 distinct nonzero values. Furthermore, Fermat’s Little Theorem tells us r^(p - 1) = (r^2)^((p - 1)/2) is 1 modulo p for any nonzero r; accordingly, the nonzero values with square roots modulo p are precisely the (p - 1)/2 many solutions to x^((p - 1)/2) = 1.

Does (-1)^((p - 1)/2) = 1? Well, by the nature of -1, this is true just in case (p - 1)/2 is even, which is to say, just in case p is of the form 4k + 1. Thus, we see that the bifurcating primes are precisely those of the form 4k + 1 (with those of the form 4k - 1 being atomic, and 2 being its own odd man out), completing our investigation.