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](https://en.wikipedia.org/wiki/Legendre%27s_three-square_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](https://mathoverflow.net/questions/77537/legendre-and-sums-of-three-squares)
* [Efficient computation of integer representation as sum of three squares](https://mathoverflow.net/questions/104322/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](https://oeis.org/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](https://oeis.org/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. 

[![enter image description here](https://i.stack.imgur.com/FqbQc.png)](https://i.stack.imgur.com/FqbQc.png)

Sometimes, there is no solution. How can I improve this to work quickly at $N \approx 10^6$ ?