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.
Please help me! Thank you in advance.
 A: 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)$?
A: 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

