Listing integers as the sum of three squares $m = x^2 + y^2 + z^2$

I would like to write an algorithm that lists all the ways a natural number $m \in \mathbb{N}$ can be written as the sum of three squares $m = x^2 + y^2 + z^2$ with integers $x,y,z \in \mathbb{Z}$.

Legendre's theorem says $m$ has such a representation if $m \neq 4^a(8b+7)$ with $a,b \in \mathbb{N}$. See on MathOverflow:

If I just wanted to count the number of representations as the sum of squares, we could find formulas in the Online Encyclopedia of Integer Sequences:

• A005875 - Theta series of simple cubic lattice; also number of ways of writing a nonnegative integer n as a sum of 3 squares (zero being allowed).
• A074590 - Number of primitive solutions to $n = x^2 + y^2 + z^2$ (i.e. with $\gcd(x,y,z) = 1$).

How do I list the sum of squares representations for each $m$? I would take a slower algorithm if it were easy to implement.

I considered just writing: $\boxed{m -x^2 -y^2 = z^2 }$ and looping over $0 < x < y < \sqrt{m}$ and checking it is perfect square.

Another possibility is keeping an array square writing square[x*x+y*y*z*z] += [[x,y,z]] storing triples in the appropriate place as I find them. That is something like $m^3$ time since I loop over $0 < x,y,z < m$.

Do any solutions jump out at you?

10/12/17 Here I take a naive approach:

N = int(M**0.5)

z = [ (a,b,c) for a in range(N) for b in range(N) for c in range(N)
if a**2 + b**2 + c**2 == M]

This rather slow, but I could get solutions in the range N=10000 in about 20 seconds. Here I run it 100 times within a half-hour. Sometimes, there is no solution. How can I improve this to work quickly at $N \approx 10^6$ ?

• 1. Are you sure you need to enumerate all ways explicitly? There can be exponentially many such representations, so listing all of them could be infeasible. Do you have a plan for how to deal with that? 2. About how big will your $m$ be? – D.W. Dec 5 '14 at 13:45
• @D.W. I am reading Local Statistics of Lattice Points on the Sphere on the arXiv. The list of integer triples on the sphere will help me implement they further calculations they describe in the paper. – john mangual Dec 5 '14 at 14:06
• I just noticed this cs.stackexchange.com/questions/2988/… – john mangual Dec 5 '14 at 14:39
• Notice that it is much quicker to calculate f (m) for many values m simultaneously. – gnasher729 Oct 14 '17 at 0:08

Your first approach is probably the simplest to implement. Its running time will be $O(m)$.

Another approach: you could loop over all $x$ such that $0 < x < \sqrt{m}$, and for each $x$ enumerate all ways to express $m-x^2$ as a sum of two squares. There's a lot of theory about how to express an integer as a sum of two squares (e.g., Lagrange's algorithm). However, this will be more complex to implement.

You are looking for solutions of $x^2 + y^2 + z^2 = m$ with no conditions for x, y, z, especially on the order. Instead solve for x < y < z and multiply the count by six (since x, y, z could be ordered in six ways), then solve $x^2 + 2y^2 = m$ for x ≠ y and multiply the count by three (you have solutions (x, y, y), (y, x, y) and (y, y, x), and then solve $3x^2 = m$ (solution x, x, x). I assume you allow x = 0, otherwise a slight change is needed. To find the count for all m, 0 ≤ m ≤ M = $10^6$:

Initialise an array with indices 0 to M to zeroes.
Loop for x = 0, 1, 2, ... as long as 3x^2 ≤ M
Loop for y = x + 1, x + 2, ... as long as x^2 + 2y^2 ≤ M
Loop for z = y + 1, y + 2, ... as long as x^2 + y^2 + z^2 ≤ M
Increase the count for x^2 + y^2 + z^2 by 6.
Loop for x = 0, 1, 2, ... as long as x^2 ≤ M
Loop for y = 0, 1, 2, ... as long as x^2 + 2y^2 ≤ M
Increase the count for x^2 + 2y^2 by 3 if x ≠ y, by 1 if x = y.

This is of course no help if you want the count for one single m (although it is faster than your original method by some constant factor, but it won't be anywhere as fast as something more clever), but I think it will be hard to beat if you want to find the count for many m. Quite similar to the sieve of Eratosthenes, which is very good at finding all primes in a large range.

Runtime will be $\Theta (M^{1.5})$ to find the count for all 0 ≤ m ≤ M. You can adapt this to find the counts for M - N ≤ m ≤ M, and that should work in $\Theta (M · max (1, N / M^{1/2}))$ which degenerates to $\Theta(M)$ for N = 0, which finds just one count.

It also finds the solutions for $M - M^{1/2} ≤ m ≤ M$ in about the same as for a single m.

PS. To find the count for a single solution, it is better to look for solutions with x > y > z: Choose $m^{1/2} ≥ x ≥ (m/3)^{1/2}$ (a smaller x cannot be the largest number in a solution). Then we need 0 ≤ y < x, and 0 ≤ z < y, so $x^2 + y^2 + 0^2 ≤ m ≤ x^2 + y^2 + (y-1)^2$. Let $m' = m - x^2$, then $y^2 ≤ m'$ and $y^2 + (y-1)^2 ≥ m'$, so $(m' / 2)^{1/2} ≤ y ≤ min (x, m'^{1/2})$. This minimises the number of pairs (x, y) to examine to about M / 18.

PS. The algorithm runs at full speed if for every pair (x, y) there are many values z. If you look for solutions in [a, b] then this is efficient as long as (b - a) is large compared to $b^{1/2}$, so you can divide the interval [a, b] in as many subintervals as you have processors, as long as the width of each subinterval is still much larger than $b^{1/2}$. If you modify the algorithm so that x > y > z, then you will have fewer pairs (x, y) and therefore more z values for each pair (x, y), so the overhead will be less.

• Could there be a way to run the same algorithm for all values of $a \leq m < b$ in parallel? If I compute all the squares once, I am really just populating all the sums of three possible values of that set that fall in a certain range. And I want to keep track of the possible combinations of (x,y,z) that do so. Is that less time / resources? – john mangual Oct 14 '17 at 19:26
• The total number of solutions for a ≤ m ≤ b is O ((b - a) b^(1/2)), so keeping track of solutions will take a lot of memory. – gnasher729 Oct 14 '17 at 23:11