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:

• Legendre and sums of three squares

• Efficient computation of integer representation as sum of three squares

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[xx+yyzz] += [[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$ ?

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.

StackExchange Computer Science Q#33887, answer score: 4