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 &lt; x &lt; y &lt; \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 &lt; x,y,z &lt; 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 &lt; x &lt; \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.