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)$?
We will just use that $\mathbb{Z}[i]$ is a factorial ring. As before $n=\prod_{i=1}^{k}p_i^{\alpha_i}$ where $p_1=2$ and $\alpha_1=0,1$ or $p_i\equiv 1 \pmod{4}$ for $i\geq 2$. Then each $p_i$ with $i\geq 2$ factors uniquely as $q_i\overline{q_i}$. Writing $n=a^2+b^2$ corresponds to factoring $n$ as $(a+bi)(a-bi)$ where $(a+bi)$ and $(a-bi)$ are corpime.
Then the statement is clear each factor of $(a+bi)$ is either $q_i^{\alpha_i}$ or $\overline{q_i}^{\alpha_i}$, and if you have $2|n$ you also have a choice for $(1+i)$ and $(1-i)$. Thus you get exactly $2^k$ such writings.
Best Answer
Asymptotics are known for all integer moments of $r_2$: for any $m \ge 1$, there is an explicit constant $a_m$ such that $\sum_{n \le x} r_2(n)^m \sim a_m\cdot x\,(\log x)^{2^{m-1}-1}$. ($m=0$ also holds if we interpret $0^m$ as $0$.)
This result, including the precise constant (and a uniform generalization to other positive binary quadratic forms), can be found in V. Blomer and A. Granville, “Estimates for representation numbers of quadratic forms”, Duke Math. J. 135 (2006). Here's the PDF.