Is there a simple way to compute the number of ways to write a positive integer as the sum of three squares? It's a standard theorem that the number of ways to write a positive integer N as the sum of two squares is given by four times the difference between its number of divisors which are congruent to 1 mod 4 and its number of divisors which are congruent to 3 mod 4.  Alternatively, there are no such representations if the prime factorization of N contains any prime of form 4k+3 an odd number of times.   If the prime factorization of N contains all such primes an even number of times, then we have
r2(N) = 4(b1 + 1)(b2 + 1)...(br+1)
where b1, ..., br are the exponents of the primes congruent to 1 mod 4 in the factorization of N.
For example, 325 = 52 × 13 can be written in 4(2+1)(1+1) = 24 ways as a sum of squares.  These are 182 + 12, 172 + 62, 152 + 102, and the representations obtained from these by changing signs and/or permuting.
Is there an analogous formula in the three-square case?  I know that an integer can be written as the sum of three squares if and only if it is not of the form 4m(8n+7).  There is a simple argument that shows that the number of ways to write all integers up to N as a sum of three squares is asymptotically 4πN3/2/3 -- representations of an integer less than N as a sum of three squares can be identified with points in the ball in R3 centered at the origin with radius N1/2.   Differentiating, a "typical" integer near N should have about 2πN1/2 representations as a sum of three squares.   From playing around with some data it looks like 
limn → ∞ #{ k ≤ n and r3(k)/k1/2 ≤ x} / n 
might be a nonzero constant.  That is, for each positive real x, the probability that a random integer k can be written in no more than x k1/2 ways approaches some constant in the open interval (0, 1) as k → ∞. 
One way to prove this (if it is in fact true) would be if there were some formula for r3(k), in terms of the prime factorization, which is why I'm curious.
(I apologize if this is something that is well-known to number theorists, although I'd appreciate a pointer if it is.  I am not a number theorist, I just play around with this sort of thing every so often and generate amusing conjectures.)
 A: See Granville-Soundararajan, ``The distribution of values of L(1, \chi),''
for (whatever is known about) the distribution of r_3(n)/sqrt(n). 
A: I just can't stop myself from putting up the following, from the MOTD on the Berkeley server:

Oct  2: Warning: Due to a known bug, the default Linux document viewer
        evince prints N*N copies of a PDF file when N copies requested.
        As a workaround, use Adobe Reader acroread for printing multiple
        copies of PDF documents, or use the fact that every natural number
        is a sum of at most four squares.

A: I put a fair amount of effort into this, just about as the recent duplicate question was being closed. So I am moving it.
I wanted to include the viewpoint of Burton Wadsworth Jones, given in his little book "The Arithmetic Theory of Quadratic Forms." The theorem, with many cases, is that the number of "primitive" or "proper" representations $R_{0}(n)$ of a number by $x^2 + y^2 + z^2,$ (meaning $\gcd (x,y,z) = 1$) is a multiple of the class number of binary quadratic forms of discriminant $-4n,$ but the multiple changes depending on congruence properties of $n.$ Also there are "ground" cases, here $n=1,$ which are done separately anyway.To get the actual number of representations for a number that is not squarefree it is  necessary to take a sum.
Let's see, if $n$ is a multiple of 4 there are no primitive representations, as $x^2 + y^2 + z^2 \equiv 0 \pmod 4$ means that $x,y,z$ are all even. But that is fine, because this also means that the number of representations of $4n$ is exactly the same as the number of representations of $n.$ Also, if
 $ n \equiv 7 \pmod 8$ there are no representations at all.
For $n > 1$ and $ n \equiv 1 \pmod 8,$  $\; \; R_{0}(n) = 12 h(-4n).$
For $ n \equiv 3 \pmod 8,$  $ \; \; R_{0}(n) = 8 h(-4n).$
For $ n \equiv 5 \pmod 8,$  $ \; \; R_{0}(n) = 12 h(-4n).$
For $ n \equiv 2 \pmod 8,$  $ \; \; R_{0}(n) = 12 h(-4n).$
For $ n \equiv 6 \pmod 8,$  $ \; \;R_{0}(n) = 12 h(-4n).$
Just to include something that is not entirely about proper representations, from the Hecke eigenform method one gets, with p an odd prime,
$$ R(p^2 n) = (p + 1 - (-n|p) ) \; \; R(n) - \; \; p \; R( n / p^2) $$
where $R(n)$ is the number of representations including both proper and improper, the Jacobi symbol $(-n|p)$ is taken to be 0 if $p | n,$ while $R(n/p^2)$ is taken to be 0 if $p^2$ does not divide $n.$ This appears in an article by Hirschhorn and Sellers called, and I think this is clever, "On representations of a number as a sum of three squares" which appeared about 1999 in a journal with the word "Discrete" in the title. I just have a preprint here.
A: Short answer: no.
Medium answer: For n square free, this is closely related to the class number of Q(sqrt{-n}); this is a result of Gauss. See Mathworld for a precise statement. This class number can then be rewritten in terms of the quadratic residue symbol. We can either use the class number formula to get an expresion as an infinite sum, or use Dirichlet's evaluation of L(1, chi) (same Wikipedia link) to give a finite expression.
When n is not square free, one can still give an answer in terms of the product of the class number and certain elementary correction factors, but the correction factors are so bad that no one wants to write them down. (By no one, I mean that the first half dozen papers I found on mathscinet wouldn't do it.)
Long answer: I did find a paper with all the details. See Theorem B of Bateman "On the representations of a number as the sum of three squares." Trans. Amer. Math. Soc. 71, (1951). 70--101. That's right, I won't write it down either :).
