# Need help understanding Sum of Squares function $r_2(n)$

I am trying to understand how to quickly find the number of two squares that can be added to form a number 'n'. This is my reference.

I have written a function that I believe gives me my proper answers, but I need it to run faster. Pretty much what it does is it goes through possible numbers below Square Root of n, and sees if those to numbers squared equals n. I then add 4 to the sum, since it includes when any number is negative (positive when squared). If you understand Java, here is my code:

static int SquaresR2(int n) {
int sum = 0;
outer:
for(int a=0; a<Math.sqrt(n)-1; a++) {
for(int b=0; b<Math.sqrt(n)-1; b++) {
if( a*a + b*b == n ) {
if(a>b) break outer;
sum+=4;
System.out.println(n+" = "+a+"^2 + "+b+"^2");
}
}
}
sum*=2;

if(Math.sqrt(n)==(int)Math.sqrt(n)) sum+=4;
return sum;
}


On Wolfram MathWorld it says that finding the Sum of Squares k=2 relates to factoring n to $n=2^{a_{0}} p_{1}^{2_{a_{1}}} ... p_{r}^{2_{a_{r}}} q_{1}^{b_{1}} ... q_{s}^{b_{s}}$ where the $p_{i}$s are primes of the form 4k+3 and the $q_{i}$s are primes of the form 4k+1. I have (almost) no idea what this means.

It also talks about $B=(b_{1}+1)(b_{2}+1)...(b_{r}+1)$, which I also have no understanding of.

I am thinking it has something to do with factoring n using primes. I am a high school senior taking Calculus 1.

-
You add four? Shouldn't you be multiplying by four? – anon Feb 11 '12 at 0:36
I have n't tried to logically verify your code but at first glance if it is speed that u seek u might wanna call Math.sqrt(n) only once and cache the result in a local stack variable u seem to call it over and over for each iteration of the two loops and afterwards u should good speed up from this alone. – Hardy Feb 11 '12 at 0:56
@AWk34 also ++a (pre-increment) should all ways a faster operation than a++ (post-increment). Although modern compilers might be able to optimize it but do n't count on it. – Hardy Feb 11 '12 at 1:08
The references given here might be of use in constructing your own implementation. – J. M. Feb 11 '12 at 1:41

I'll leave alone the questions of (a) if actually computing the prime factorization of numbers is realistic or practical for this algorithm or (b) why what MW has down is true. Instead, I will just focus on showing what MW is saying. I assume you know about primes and prime factorizations.

Every positive integer has a unique factorization into prime powers. For example,

$$60=2^2\cdot3^1\cdot5^1,\qquad 1225=5^2\cdot7^2,\qquad 7118280=2^3\cdot3^4\cdot5\cdot13^3.$$

All odd primes are either of the form $4k+1$ or $4k+3$. Here's a table for the first few primes $p$:

$$\begin{array}{|r|l|} \hline p & 4(\circ)+\square \\ \hline 3 & 4(0)+\color{Blue}3 \\ 5 & 4(1)+\color{Red}1 \\ 7 & 4(1)+\color{Blue}3 \\ 11 & 4(2)+\color{Blue}3 \\ 13 & 4(3)+\color{Red}1 \\ 17 & 4(4)+\color{Red}1 \\ 19 & 4(4)+\color{Blue}3 \\ 23 & 4(5)+\color{Blue}3 \\ \hline \end{array}$$

So we can take the primes in any prime power factorization and categorize them as either $\color{Green}2$ or of one of the forms $\color{Red}{4k+1}$ or $\color{Blue}{4k+3}$. Let's look at the previous examples again:

$$60=\color{Green}{2^2} \cdot \color{Blue}{3^1} \cdot \color{Red}{5^1},\qquad 1225=\color{Red}{5^2} \cdot \color{Blue}{7^2},\qquad 7118280=\color{Green}{2^3} \cdot \color{Blue}{3^4} \cdot \color{Red}{5} \cdot \color{Red}{13^3}.$$

What MW says is now

• The power of $2$ in $n$'s prime factorization does not affect the value of $r_2(n)$.
• If any power of a blue prime ($4k+3$) in $n$'s factorization is odd e.g. $3^1$ in $60$, then $r_2(n)=0$.
• Otherwise, in no particular order, let $b_1, b_2, b_3, \dots$ be the powers of the red primes ($4k+1$) in the factorization of $n$. We then have the equality $r_2(n)=4(b_1+1)(b_2+1)\cdots$.

Can you use these facts and the given examples to compute $r_2(60)$, $r_2(1225)$ and $r_2(7118280)$?

-
You may be interested in not allowing negative numbers, and in counting for example the representations $13=2^2+3^2$ and $13=3^2+2^2$ as the same. Let $n>0$. (i) If $r_2(n)$ is divisible by $8$, the number of essentially distinct representations is $r_2(n)/8$ (ii) Otherwise, it turns out that $r_2(n)$ is divisible by $4$ but not by $8$. Let $r_2(n)=8k+4$. Then the number of essentially distinct representations is $k+1$. – André Nicolas Feb 11 '12 at 1:29
@André: the sum of two squares function defined in MathWorld indeed accounts for order, since $r_2(1)$ is $4$. The function you are considering is a different, but still interesting, beastie... – J. M. Feb 11 '12 at 2:32
@J.M.: I was not sure which one the OP was really after. If the motivation is geometric, then it is probably $r_2(n)$. – André Nicolas Feb 11 '12 at 2:38
The motivation is geometric. Specifically it is for finding lattice points of a circle. – Awk34 Feb 11 '12 at 4:11
Hooray! I think that I have gotten my $r_{2}(n)$ function works for small numbers, but still not as fast as I would like. I tried these tests: f(1328125) = 180?, f(84246500) = 420?, f(248431625) = 420?, f(N) = SquaresR2($n^{2}$). The first two work, but the second one outputs 0. – Awk34 Feb 11 '12 at 5:40

I'm not sure about quickness (you'll have to do your own tests), but here's a method you might want to consider, based on this one line in the OEIS:

Euler transform of period $4$ sequence $4,-6,4,-2,\dots$ - Michael Somos, Jul 19 2004

More concretely, consider the function

$$a_n=\begin{cases}4&\text{if }n\text{ odd}\\-2&\text{if }n\bmod 4=0\\-6&\text{if }n\bmod 4=2\\\end{cases}$$

and form the function

$$c_n=\sum_{d\mid n} d a_d$$

where $\sum\limits_{d\mid n}$ means that one should sum over all integers $d$ that are divisors of $n$.

We can then compute $r_2(n)$ via the recursion relation

$$r_2(n)=\frac1{n}\left(c_n+\sum_{j=1}^{n-1} c_j\; r_2(n-j)\right)$$

with the initial condition $r_2(1)=4$.

Implementing $a_n$ in code should be a snap; implementing $c_n$ would only require checking for the divisors of $n$ and adding up terms corresponding to those; implementing the actual sum-of-two squares function $r_2(n)$ will require caching and whatever other tricks are needed for recursion.

Here is a sample Mathematica implementation:

c[k_Integer] := DivisorSum[k, # Piecewise[{{4, OddQ[#]},
{-6, Mod[#, 4] == 2},
{-2, Mod[#, 4] == 0}}] &];

r2[1] = 4;
r2[k_Integer] := r2[k] = (c[k] + Sum[c[j]*r2[k - j], {j, 1, k - 1}])/k


Check for first fifty positive integers:

Apply[And, Thread[(b /@ Range[50]) - SquaresR[2, Range[50]] == 0]]
True

-