Number of integer solutions of $x^2 + y^2 = k$

This is worked out in detail in many an introductory textbook on Number Theory. Here is a summary of the treatment in An Introduction to the Theory of Numbers (5th Edition) by Ivan Niven, Herbert S. Zuckerman, and Hugh L. Montgomery (Wiley, 1991, section 3.6).

Theorem (Fermat) Let $n$ be a positive integer, and write $n$ in the form $$n = 2^{\alpha}\prod_{p\equiv 1\pmod{4}} p^{\beta} \prod_{q\equiv 3 \pmod{4}} q^{\gamma}.$$ with $p$ and $q$ primes. Then $n$ can be expressed as a sum of two squares if and only if all the exponents $\gamma$ are even.

Now, define the following:

• $R(n)$: the number of ordered pairs $(x,y)$ of integers such that $x^2+y^2=n$.
• $r(n)$: the number of ordered pairs $(x,y)$ such that $\gcd(x,y)=1$ and $x^2+y^2=n$ (called "proper representations of $n$");
• $P(n)$: the number of representations of $n$ by the form $x^2+y^2$ with $\gcd(x,y)=1$, $x\gt 0$ and $y\geq 0$;
• $N(n)$: the number of solutions to the congruence $s^2\equiv -1\pmod{n}$.

Theorem. A positive integer $n$ is properly represented by the form $x^2+y^2$ if and only if $-4$ is a square modulo $4n$.

In particular, since $-4$ is a square modulo $8$ but not modulo $16$, $n$ may be divisible by $2$ but not by $4$. If $p$ is an odd prime of the form $4k+1$, then $-4$ is a square modulo $p$, and by Hensel's Lemma it lifts to a solution modulo $p^k$ for an y $k$; so $n$ may be divisible by arbitrary powers of primes of the form $4k+1$. On the other hand, if $p$ is a prime that divides $n$ and is congruent to $3$ modulo $4$, then $-4$ is not a square modulo $p$. So we get:

Theorem. A positive integer $n$ is properly representable as a sum of two squares if and only if the prime factors of $n$ are all of the form $4k+1$, except for the prime $2$ which may occur to at most the first power.

Now suppose that $n$ is positive and $n=x^2+y^2$ is an arbitrary representation. Set $g=\gcd(x,y)$; then $g^2|n$, so $n = g^2m$, $\gcd(\frac{x}{g},\frac{y}{g}) = 1$, so $m$ is properly represented as the sum of the squares of $\frac{x}{g}$ and $\frac{y}{g}$. Note that $g$ may have prime factors congruent to $3$ modulo $4$, which divide $n$ to an even power, and the power of $2$ that divides $n$ may be arbitrary as well.

Theorem. Suppose that $n\gt 0$. Then $P(n)=N(n)$, $r(n) = 4N(n)$ and $R(n)=\sum r\left(\frac{n}{d^2}\right)$, where the sum is taken over all positive $d$ such that $d^2|n$.

Theorem. Let $n$ be a positive integer and write $$ n = 2^{\alpha}\prod_{p}p^{\beta}\prod_q q^{\gamma}$$ where $p$ runs over prime divisors of $n$ congruent to $1$ modulo $4$ in the first product, and $q$ runs over prime divisors of $n$ of the form $4k+3$ in the second. If $\alpha=0$ or $1$, and all $\gamma$ are $0$, then $r(n) = 2^{t+2}$, where $t$ is the number of primes $p$ of the form $4k+1$ that divide $n$ (counted with multiplicity); otherwise, $r(n)=0$. If all the $\gamma$ are even, then $R(n) = 4\prod_p(\beta+1)$. Otherwise, $R(n)=0$.

$25$ has two prime divisors: $5$ and $5$, and you have to count them both because they are both congruent to $1$ modulo $4$. Now multiply by $4$. You have $8$.

This counts the decompositions $$\begin{array}{cc} 3^2+4^2,& 3^2+(-4)^2, \\ (-3)^2+4^2, & (-3)^2+(-4)^2, \\ 4^2+3^2, & 4^2+(-3)^2, \\ (-4)^2+3^2, & (-4)^2+(-3)^2 \end{array} $$

As you can see, multiplying by $4$ at the end is done to cater for signs: if you are only interested in "positive" decompositions, then don't do that. Also, this counts $a^2+b^2$ as different from $b^2+a^2$. If you don't want that, well, you have to substract the number $N$ of decompositions of the form $a^2+a^2$ (which is easy to compute...) divide by $2$, and then add $N$ back.

Finally, if you want to know why this works, well, there are very few better things than to consult an introduction to number theory!

A simpler way to look at this would be to consider the factorization of $\displaystyle n$ in Gaussian Integers: $\displaystyle \mathbb{Z}[i] = \{a + bi \ \ | \ a \in \mathbb{Z}, \ b \in \mathbb{Z}\}$.

It is well known that that

• Gaussian integers have unique factorization property (upto units $\displaystyle \pm 1, \pm i$).
• Primes (of $\displaystyle \mathbb{Z}$) of the form $\displaystyle 4k+3$ are also primes in $\displaystyle \mathbb{Z}[i]$.
• Primes (of $\displaystyle \mathbb{Z}$) of the form $\displaystyle 4k+1$ factorize into $\displaystyle w w'$ for a prime $\displaystyle w$ (in $\displaystyle \mathbb{Z}[i]$) and $\displaystyle w'$ is the conjugate of $\displaystyle w$. In fact, if $\displaystyle w = a+ib$, then the prime is of the form $\displaystyle a^2 + b^2$.
• $\displaystyle 2 = (1+i)(1-i)$ and $\displaystyle 1+i$ is prime.

Now if $\displaystyle n = x^2 + y^2$ then $\displaystyle n = (x+iy)(x-iy)$, which corresponds to factorization of $\displaystyle n$ in $\displaystyle \mathbb{Z}[i]$, which you can get from the prime factorization of $\displaystyle n$ in $\mathbb{Z}[i]$ by choosing a subset of primes to multiply, to get $\displaystyle x + iy$. The conjugates of those primes go toward $\displaystyle x-iy$. The remaining need to be a perfect square to distribute once each into $\displaystyle x+iy$ and $\displaystyle x-iy$. Note this implies that primes of the form $\displaystyle 4k+3$ (which are also primes in $\displaystyle \mathbb{Z}[i]$) need to have an even power in the factorization.

So it basically amounts to finding out the different factorizations of the perfect squares which are factors of $\displaystyle n$ (after ignoring primes of the form $\displaystyle 4k+3$).

Note: It is the units which gives the multiple of $4$, if you count $(x,y), (y,x), (-x,y),(y,-x)$.

(Of course, this is not a formal proof...)

Hope the helps.