Number of integer solutions of $x^2 + y^2 = k$ I'm looking for some help disproving an answer provided on a StackOverflow question I posted about computing the number of double square combinations for a given integer.  
The original question is from the Facebook Hacker Cup  

Source: Facebook Hacker Cup
  Qualification Round 2011
A double-square number is an integer $X$
  which can be expressed as the sum of
  two perfect squares. For example, 10
  is a double-square because $10 = 3^2 + 1^2$.
  Given $X$, how can we determine the
  number of ways in which it can be
  written as the sum of two squares? For
  example, $10$ can only be written as $3^2 + 1^2$
  (we don't count $1^2 + 3^2$ as being different). On the other hand,
  $25$ can be written as $5^2 + 0^2$ or as
  $4^2 + 3^2$.
You need to solve this problem for $0 \leq X \leq 2,147,483,647$.
Examples:
$10 \Rightarrow 1$
$25 \Rightarrow 2$
$3 \Rightarrow 0$
$0 \Rightarrow 1$
$1 \Rightarrow 1$

In response to my original question about optimizing this for F#, I got the following response which I'm unable to confirm solves the given problem correctly.  

Source: StackOverflow answer by Alexandre C.
Again, the number of integer solutions
  of $x^2 + y^2 = k$ is four times the
  number of prime divisors of $k$ which
  are equal to $1 \bmod 4$.
Knowing this, writing a program which
  gives the number of solutions is easy:
  compute prime numbers up to $46341$ once
  and for all.
Given $k$, compute the prime divisors of
  $k$ by using the above list (test up to
  $\sqrt{k}$). Count the ones which are
  equal to $1 \bmod 4$, and sum. Multiply
  answer by $4$.

When I go through this algorithm for $25$, I get $8$ which is not correct.  

  
*
  
*For each prime factor (pf) of $25$ ($5$, $5$ are prime factors of $25$)
  
  
*
  
*if pf % $4$ = $1$ (true for both $5$'s), add $1$ to count
  
  
*return $4$ * count (count would be $2$ here). 
  
  
  So for $25$, this would be $8$

 A: $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: 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.
A: 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$.
