Why is it harder to prove which integers are sums of three squares rather than sums of two squares or four squares? Background: Let $n$ be an integer and let $p$ be a prime. If $p^{e} || n$, we write $v_{p}(n) = e$. A natural number $n$ is a sum of two integer squares if and only if for each prime $p \equiv 3 \pmod 4$, $v_{p}(n)$ is even. Every natural number is a sum of four squares. A natural number $n$ is a sum of three squares if and only if it is not of the form $4^{k}u$ where $u \equiv 7 \pmod 8$. 
I would like to know why it is harder to prove the above result for sums of three squares as opposed to sums of two squares or four squares.
I've heard somewhere that one way to see this involves modular forms... but I don't remember any details. I would also like to know if there is a formula for the number of ways of representing a natural number n as a sum of three squares (or more generally, $m$ squares) that is similar in spirit to the formulas for the number of ways of representing a natural number as the sum of two squares and four squares. 
 A: The modular forms explanation is basically due to the fact
that $3$ is odd and so the generating function for representations
of sums of three squares is a modular form of half-integer weight.
In general if $r_k(n)$ is the number of representations of $n$
as a sum of $k$ squares then
$$\sum_{n=0}^\infty r_k(n)q^n=\theta(z)^k$$
where $q=\exp(\pi i z)$ and 
$$\theta(z)=1+2\sum_{n=1}^\infty q^{n^2}.$$
Then $f_k(z)=\theta(z)^k$ is a modular form of weight $k/2$ for the group
$\Gamma_0(4)$. This means that
$$f_k((az+b)/(cz+d))=(cz+d)^{k/2}f_k(z)$$
whenever the matrix $\begin{pmatrix}a&b\\\\c&d\end{pmatrix}$
lies in $\Gamma_0(4)$, that is $a$, $b$, $c$ and $d$ are integers, $4\mid c$
and $ad-bc=1$.
This definition is easy to understand when $k$ is even, but for odd
$k$ one needs to take the correct branch of $(cz+d)^{k/2}$, and
this is awkward. The space of modular forms of weight $k/2$ is
finite-dimensional for all $k$, and is one-dimensional for small enough $k$.
For these small $k$ the space is spanned by an "Eisenstein series".
Computing the Eisenstein series isn't too hard for even $k$, but is
much nastier for odd $k$ where again square roots need to be
dealt with. See Koblitz's book on modular forms and elliptic functions
for the calculation for $k\ge5$ odd. The calculation for $k=3$
is even nastier as the Eisenstein series does not converge absolutely.
In fact the cases where $k$ is divisible by $4$ are even easier, as
even weight modular forms behave nicer.
For large $k$, Eisenstein series are no longer enough, one needs
also "cusp forms". While fascinating, cusp forms have coefficients
which aren't given by nice formulae unlike Eisenstein series.
Of course there is a formula for $r_3(n)$, due to Gauss in his
Disquisitiones Arithmeticae. It involves class numbers of quadratic
fields (or to Gauss numbers of classes of integral quadratic forms).
A: An explanation from a different direction than you may be expecting: one reason why the sum-of-three-squares problem is so much harder is that there's no normed division algebra of dimension three.  A key element of both the two-squares and four-squares proof is expressing the product of two numbers in the relevant form as another number in that form; these product formulas 'come from' the formulas for multiplying complex and quaternionic numbers respectively.  That is, taking norms of $(a+bi)(c+di) = (ac-bd) + (ad+bc)i$ gives rise to the formula $(a^2+b^2)(c^2+d^2) = (ac-bd)^2+(ad+bc)^2$, and the equivalent formula for quaternions gives Euler's four-square identity.  John Baez wrote a particularly interesting article on why you don't keep getting division algebras as you keep doubling the dimension; you should be able to find it on his homepage if you're curious.
