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