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

Question

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[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$ ?

Answer 1

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.

Answer 2

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.

PS. You loop over x, and for each x you want to count the pairs (y, z) such that $y^2+z^2=m-x^2$. Consider that a^2 modulo 16 is always one of 0, 1, 4 or 9. So $y^2+z^2$ modulo 16 is one of 0, 1, 2, 4, 5, 8, 9, 10, 13. For seven possible values of $m-x%2$ modulo 16 we don’t need to do anything. For the others we can significantly reduce the number of values y to examine. And we can consider other modulo values, significantly reducing the amount of work needed.

Answer 3

Code :

N=int(M**0.5)
z=[]
for a in range(int(N/3**0.5),N+1):
    for b in range(int(((M-a**2)/2)**0.5),min(a+1,int((M-a**2)**0.5)+1)): 
        if int((M-a**2-b**2)**0.5) == (M-a**2-b**2)**0.5: 
            z.append((a,b,(M-a**2-b**2)**0.5))

We assume that a>=b>=c => a>=(M/3)**0.5,

Since b>=c : b²+c²=M-a² => 2*b²>= (M-a²) => sqrt(b>= ((M-a²)/2)).

You fix a and b and then you test if c = (M-a²-b²)**0.5) is integer. You loop twice : First for a 2/3*sqrt(M) elements and secondly for b sqrt(M)*(sqrt(2)-1)=> your complexity is O(M)

Answer 4

I came across this question while considering this Mathematics SE question, for which one solution is numbers that are the sum of $k$ squares in multiple ways, with the added criteria of a common sum of square roots and at least one distinct $k$-tuplet. Since three squares was the simplest option, I wrote a Python program to investigate these; for large limits it is far quicker to use a sieve approach, if you can stand the large storage requirement.

This is the portion of code that generate the various sums of squares; note that they are in non-decreasing order, so a general case of $a<b<c$ would imply different coordinates by re-ordering.

lim = int(input('Search limit: '))
rlim = int(lim**0.5)
sqval = {i:i*i for i in range(rlim+1)}

# sieve
allz = []
for a in range(1,rlim+1):
    for b in range(a,rlim+1):
        ab = sqval[a]+sqval[b]
        if ab > lim:
            break
        for c in range(b, rlim+1):
            abc = ab+sqval[c]
            if abc > lim:
                break
            allz.append( (abc, a+b+c, a,b,c) )
# print(len(allz))
allz.sort()

The process simply iterates over the possible square root values until the sum of squares goes out of range, and stores all results (i.e. all ways to generate all sums in the range).

List creation time for a sum-of-squares limit of 100,000 was about 5 seconds, most of which was the final sort. It would not be hard to modify this to generate values in a range, rather than starting at zero; you would just need to find/update the right start value for $c$ in the generating loop.