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.

enter image description here

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

  • $\begingroup$ 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? $\endgroup$
    – D.W.
    Dec 5, 2014 at 13:45
  • $\begingroup$ @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. $\endgroup$ Dec 5, 2014 at 14:06
  • $\begingroup$ I just noticed this cs.stackexchange.com/questions/2988/… $\endgroup$ Dec 5, 2014 at 14:39
  • $\begingroup$ Notice that it is much quicker to calculate f (m) for many values m simultaneously. $\endgroup$
    – gnasher729
    Oct 14, 2017 at 0:08

4 Answers 4


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.

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.

  • $\begingroup$ 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? $\endgroup$ Oct 14, 2017 at 19:26
  • $\begingroup$ 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. $\endgroup$
    – gnasher729
    Oct 14, 2017 at 23:11

Code :

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: 

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)


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:
        for c in range(b, rlim+1):
            abc = ab+sqval[c]
            if abc > lim:
            allz.append( (abc, a+b+c, a,b,c) )
# print(len(allz))

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.

  • $\begingroup$ Can you explain the algorithm in words? This is not a programming site. $\endgroup$ Mar 29, 2021 at 12:43
  • $\begingroup$ @YuvalFilmus as requested; it's very simple but such a massive speed improvement over that described in the question that it seemed worth answering. $\endgroup$
    – Joffan
    Mar 29, 2021 at 13:49
  • $\begingroup$ This is also essentially a brute-force approach. What is its running time as a function of $M$, the integer in whose number of representations as a sum of three squares we are interested? $\endgroup$ Mar 29, 2021 at 15:35
  • $\begingroup$ Yes, I agree that it is brute force. My recommendation here is clearly that given the apparent use case, it is better not to calculate for individual $M$ but to sieve for a set of values. It doesn't use any special properties of sums of squares. You might reasonably say this is engineering rather than science. :-) $\endgroup$
    – Joffan
    Mar 29, 2021 at 15:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.