Sum of Three Squares and Factors of Two

Background

I've been working a lot recently with the sums of squares (I have a previous post about finding the Exhaustive List of the Sums of Four-Squares). This work has been for some research and also for Project Euler. However, I came across an interesting result: If $n$ is even and $a^2 + b^2 + c^2 = n^2$, then $a,b,$ and $c$ are also even.

To kind of tip-toe around the issue, we know that there is at least one solution $a,b,c$ such that this is true: $a=n, b=c=0$. Others may and do exist, but we won't fret too much right now about what they are. All we need to know is that there is at least one solution. Since all natural numbers are either even or odd, we have four possibilities: all $a,b,c$ are odd, two are odd and one is even, one is odd and two are even, or all three are even.

We can easily rule out the case for which one is odd and three are odd. Because we know $n$ is even, $n^2$ must also be even. However, there are an odd number of odds on the other side (remember an odd squared is still odd), which means the entire other side is odd. And, since an odd cannot equal an even, we have a contradiction. However, in order to eliminate the case where two are odd, we must do a little more work.

Assume, without loss of generality, that $a$ and $b$ are odd and $c$ is even and we will prove a contradiction. Let $d,e,f,g\in\mathbb{N}$ be such that $$\begin{array}{c} a = 2d + 1 \\ b = 2e + 1 \\ c = 2f \\ n = 2g \end{array} $$ We know that. $$n^2 = (2g)^2 = 4g^2$$ And we know that $$\begin{array}{r l} n^2 &= a^2 + b^2 + c^2 \\ &= (2d+1)^2 + (2e+1)^2 + (2f)^2 \\ & = 4d^2 + 4d + 1 + 4e^2 + 4e + 1 + 4f^2 \\ &= 4(d^2 + d + e^2 + e + f^2) + 2 \end{array}$$ Putting the two together: $$4g^2 = 4(d^2 + d + e^2 + e + f^2) + 2$$ Clearly, the left hand side is a multiple of $4$, however, the righthand side is not. Therefore, these two cannot be equal and we cannot have two odd squares to sum. We won't try to eliminate the last possibility (that all three be even) because we know that at least one solution exists and it has to fall into one of those categories.

Now, this probably doesn't seem like that interesting of a result. And honestly speaking, it's really not particularly interesting. However, we can generalizes this:

Theorem

If $2^k\vert n$ for some $k\in\mathbb{N}$ and $a^2 + b^2 + c^2 = n^2$ for $a,b,c\in \mathbb{N}$, then $2^k\vert a,b,c$.

Proof

For $k=1$, this simplifies to the above proof and $k=0$ is trivial. As such, we have base cases and we will proceed by induction. So, we assume that if $2^{k-1}\vert n$ for some $k\in \mathbb{N}$ and $a^2 + b^2 + c^2 = n^2$, then $2^{k-1}\vert a,b,c$. Now assume that $2^k\vert n$ and that $a^2 + b^2 + c^2 = n^2$. We know that there exists $m$ such that $n=2^km$. We also can show that $$n = 2^km = 2^{k-1}(2m)$$ which means that $2^{k-1}\vert n$. And, since $a^2 + b^2 + c^2 = n^2$ by our induction hypothesis, we can conclude that $2^{k-1}\vert a,b,c$. This implies there exists $i,j,l\in \mathbb{N}$ such that $$\begin{array}{c} a = 2^{k-1}i \\ b = 2^{k-1}j \\ c = 2^{k-1}l \\ \end{array}$$ We can now classify $a,b,c$ based on the respective parity of $i,j,l$ (because again, $i,j,l$ must either be even or odd). We, again, have four possibilities: all three are odd, two are odd and one is even, one is odd and two are even, and all three are even.

Case 1

Assume that all three are odd. Therefore, there exists $i_0, j_0, l_0$ such that $$\begin{array}{c} i = 2i_0 + 1 \\ j = 2j_0 + 1 \\ l = 2l_0 + 1 \end{array}$$ and furthermore $$\begin{array}{c} a = 2^{k-1}(2i_0 + 1) = 2^ki_0 + 2^{k-1} \\ b = 2^{k-1}(2j_0 + 1) = 2^kj_0 + 2^{k-1} \\ c = 2^{k-1}(2l_0 + 1) = 2^kl_0 + 2^{k-1} \\ \end{array}$$ Consider $$\begin{array}{r l} n^2 &= a^2 + b^2 + c^2 \\ & = (2^ki_0 + 2^{k-1})^2 + (2^kj_0 + 2^{k-1})^2 + (2^kl_0 + 2^{k-1})^2 \\ &= 2^{2k}(i_0^2+i_0+j_0^2+j_0+l_0^2+l_0) + 3(2^{2k-2}) \\ \end{array}$$ Note that this is not divisible by $2^{2k}$ because of the last term. On the contrary, note that $$n^2 = (2^km)^2 = 2^{2k}m^2$$ and therefore $2^{2k}\vert n^2$; a contradiction. Therefore, all three cannot be odd.

Case 2

Assume that two are odd and one is even. Therefore there exists $i_0,j_0, l_0$ such that $$\begin{array}{c} i = 2i_0 + 1 \\ j = 2j_0 + 1 \\ l = 2l_0 \end{array}$$ As before, we can substitute and write $a,b,c$ in terms of $i_0, j_0,l_0$ and as a result $n^2$ in terms of $i_0, j_0,l_0$. $$n^2 = 2^{2k}(i_0^2+i_0+j_0^2+j_0+l_0^2) + 2(2^{2k-2})$$ Note that this is not divisible by $2^{2k}$ because of the last term only having a factor of $2^{2k-1}$. But again, because know that $2^{2k}\vert n^2$ we have a contradiction. Therefore, two cannot be odd.

Case 3

Assume that one is odd and two are even. Therefore there exists $i_0,j_0, l_0$ such that $$\begin{array}{c} i = 2i_0 + 1 \\ j = 2j_0 \\ l = 2l_0 \end{array}$$ Once again, we can substitute and write $n^2$ in terms of $i_0, j_0,l_0$. $$n^2 = 2^{2k}(i_0^2+i_0+j_0^2+l_0^2) + 2^{2k-2}$$ Note that this is not divisible by $2^{2k}$ because of the last term. And again, because know that $2^{2k}\vert n^2$ we have a contradiction. Therefore, one cannot be odd.

Case 4

That only leaves the option that all three are even. Since $i,j,k$ are even, and we have $i_0,j_0, l_0$ such that $$\begin{array}{c} i = 2i_0 \\ j = 2j_0 \\ l = 2l_0 \end{array}$$ we can now write $$\begin{array}{c} a = 2^{k-1}(2i_0) = 2^ki_0 \\ b = 2^{k-1}(2j_0) = 2^kj_0 \\ c = 2^{k-1}(2l_0) = 2^kl_0 \\ \end{array}$$ and we have shown that $2^k\vert a,b,c$!

Implications

What does this mean, though? It means that if we want to find the three integers whose squares sum to $n^2$(such as lattice points on a sphere), then we can omit multiples of two from our search and work with smaller numbers. This saves a great deal of computing time.

For instance, if we know that $2^k\vert n$ and we want to find $a,b,c$ such that $n^2 = a^2+b^2+c^2$, we can really just search for $d,e,f$ such that $m^2=d^2+e^2+f^2$ ($n = 2^km$). We jump back to the original problem by multiplying by $2^{2k}$ on both sides: $$\begin{array}{r l} m^2 & = d^2 + e^2 + f^2\\ 2^{2k}m^2 & = 2^{2k}(d^2 + e^2 + f^2) \\ (2^km)^2&= 2^{2k}d^2 + 2^{2k}e^2 + 2^{2k}f^2 \\ n^2 & = (2^kd)^2 + (2^ke)^2 + (2^kf)^2\\ n^2 & = a^2 + b^2 + c^2 \end{array} $$ Where clearly $d,e,f$ is the respective result of the division of $a,b,c$ by $2^k$.

Exhaustive List of the Sums of Four-Squares

I'd recently been asked to assist on some research on equilateral triangles in $\mathbb{Z}^4$. The project required coming up with a list of all of the lattice points $(i,j,k,l)$ on a hypersphere of radius $\sqrt{d}$ for $d\in\mathbb{N}$. Basically, this reduces to the problem of finding all sets of four integers $(i,j,k,l)$ such that $i^2+j^2+k^2+l^2=d$.

This is an old problem and in 1770 Lagrange proved that every natural number $d$ can be represented as the sum of four squares^[1]. So we know that we should always find at least one. Additionally, in 1834 Jacobi determined how many solutions exist^[2]: $$N(d)=8\sum_{4\nmid m|d}m$$ We can use that to check to see if our program can find all the points (at least, the right number of points). Of course for the smaller values of $d$, a complete check on the lattice can be performed from $-\sqrt{d}$ to $\sqrt{d}$ in all four cardinal directions. However, this has a complexity of $O(d^2)$.

from math import sqrt

def sum_squares(d):

   sqrt_d = int(sqrt(d))

   ## This list will hold all of the points
   points = []

   ## Nested for loops to check
   ## Add 1 because of the exclusivity of the bound of range() function
   for i in range(-sqrt_d,sqrt_d + 1):
      for j in range(-sqrt_d,sqrt_d + 1):
         for k in range(-sqrt_d,sqrt_d + 1):
            for l in range(-sqrt_d,sqrt_d + 1):
               if i**2 + j**2 + k**2 + l**2 == d:
                  points.append([i,j,k,l])

   return points

However, this gets really slow really quickly! Especially if I wanted to run this for $d$ in the tens of thousands. Better results can be obtained when we realize that our $l$ loop is completely unnecessary. Since we would already have a guesses for $i, j,$ and $k$ and we know $d$, we can solve for $l$ and check to see if it is an integer. This reduces complexity to $O(d^{1.5})$. We have to do a few things to make sure that the sqrt() function only takes a positive number and that we include both the positive and negative root, but the updated code is as follows.

from math import sqrt

def sum_squares(d):

   ## This list will hold all of the points
   points = []

   ## Nested for loops to check
   ## Add 1 because of the exclusivity of the bound of range() function
   for i in range(-sqrt_d,sqrt_d + 1):
      for j in range(-sqrt_d,sqrt_d + 1):
         for k in range(-sqrt_d,sqrt_d + 1):
            l_squared = d - i**2 - j**2 - k**2
            if l_squared >= 0:
               l = sqrt(l_squared)
               if l.is_integer():
                  points.append([i,j,k,int(l)])
                  if l != 0:
                     points.append([i,j,k,-int(l)])

   return points

Even better results can be obtained if you only check the first orthant and permute the sign of the results into the other 15 orthants. This will run in slightly more than $1/16$ the amount of time, but will still have the same complexity. However, depending on how many first orthant points we find, the permutation can take a significant amount of time and I will address the permutations more in a different post. Looking in the first orthant only, we eliminate the negative root for $l$ in the last line of the nested loops.

from math import sqrt

## This list contains the 16 sign permutations
signs = [[ 1, 1, 1, 1], [ 1, 1, 1,-1], [ 1, 1,-1, 1], [ 1, 1,-1,-1],
         [ 1,-1, 1, 1], [ 1,-1, 1,-1], [ 1,-1,-1, 1], [ 1,-1,-1,-1],
         [-1, 1, 1, 1], [-1, 1, 1,-1], [-1, 1,-1, 1], [-1, 1,-1,-1],
         [-1,-1, 1, 1], [-1,-1, 1,-1], [-1,-1,-1, 1], [-1,-1,-1,-1]]

def sum_squares(d):

   ## Add 1 because of the exclusivity of the bound of range() function
   sqrt_d = int(sqrt(d)) + 1

   ## This list will hold the points in the first orthant
   points = []

   ## This list will hold the points from all orthants
   final = []

   ## Nested for loops to check
   for i in range(sqrt_d):
      for j in range(sqrt_d):
         for k in range(sqrt_d):
            l_squared = d - i**2 - j**2 - k**2
            if l_squared >= 0:
               l = sqrt(l_squared)
               if l.is_integer():
                  points.append([i,j,k,int(l)])

   ## For each point found in the first orthant
   for m in range(len(points)):

      ## This list will hold all sign permutations for the mth point
      ## in the first orthant
      point_signs = []

      ## For each of the possible sign permutations
      for i in range(len(signs)):
         tmp = [points[m][0] * signs[i][0], points[m][1] * signs[i][1],
                points[m][2] * signs[i][2], points[m][3] * signs[i][3]]

         ## To avoid errors with points on the axes
         ## i.e. 0 * 1 = 0 * -1
         if tmp not in point_signs:
            point_signs.append(tmp)

      final = final + point_signs

   return final

However, this is still has the same complexity. For larger values of $d$ something different must be done to reduce complexity! One such way to reduce complexity is to search for only the basic representation of the lattice points in the first orthant. The basic representation is the representation $(i,j,k,l)$ such that $i\geq j\geq k\geq l \geq 0$. For example, rather than finding $\{(2,1,1,0), (2,1,0,1), (2,0,1,1), \dots\}$, we only need to find $(2,1,1,0)$ and permute the order to produce the rest.

This further restricts the bounds on our for loops. For instance, we know that $\sqrt{d/4} \leq i \leq \sqrt{d}$. The upper bound is the same as before, however we have now bounded $i$ from below! To see why, we set $i=j=k=l=i_{min}$. This is $i$ at it's lowest possible value and $j,k,$ and $l$ at their largest. This yields $$d=i^2+j^2+k^2+l^2 = 4{i_{min}}^2\implies i_{min} = \sqrt{d/4} $$ Each of the subsequent indices are then bounded above by both the previous index and the difference between d and the sum of the squares of all the previous indices. For instance, $k \leq j$ and $k \leq d - (i^2 + j^2)$. A stronger statement is therefore $k \leq \mathrm{min}(j, d - (i^2 + j^2))$.

from math import sqrt

signs = [[ 1, 1, 1, 1], [ 1, 1, 1,-1], [ 1, 1,-1, 1], [ 1, 1,-1,-1],
         [ 1,-1, 1, 1], [ 1,-1, 1,-1], [ 1,-1,-1, 1], [ 1,-1,-1,-1],
         [-1, 1, 1, 1], [-1, 1, 1,-1], [-1, 1,-1, 1], [-1, 1,-1,-1],
         [-1,-1, 1, 1], [-1,-1, 1,-1], [-1,-1,-1, 1], [-1,-1,-1,-1]]

orders = [[0,1,2,3], [0,1,3,2], [0,2,1,3], [0,2,3,1], [0,3,1,2], [0,3,2,1],
          [1,0,2,3], [1,0,3,2], [1,2,0,3], [1,2,3,0], [1,3,0,2], [1,3,2,0],
          [2,0,1,3], [2,0,3,1], [2,1,0,3], [2,1,3,0], [2,3,0,1], [2,3,1,0],
          [3,0,1,2], [3,0,2,1], [3,1,0,2], [3,1,2,0], [3,2,0,1], [3,2,1,0]]

def sum_squares(d):
   
   ## Add 1 because of the exclusivity of the bound of range() function
   sqrt_d = int(sqrt(d)) + 1
   sqrt_d4 = int(sqrt(d/4))
   
   ## This list will hold the basic points in the first orthant
   points = []

   ## This list will hold the points from all orthants
   final = []
   
   ## Nested for loops to check
   for i in range(sqrt_d4, sqrt_d):
      sum_i = d - i**2
      for j in range(min(i,int(sqrt(sum_i)))+1):
         sum_j = sum_i - j**2
         for k in range(min(j,int(sqrt(sum_j)))+1):
            sum_k = sum_j - k**2
            l = int(sqrt(sum_k))
            if l**2 == sum_k and l <= k:
               points.append([i,j,k,l])

   ## For each basic point found in the first orthant
   for m in range(len(points)):
      ## This list will hold all order permutations for the
      ## ith sign permutation
      point_orders = []

      ## For each of the possible order permutations
      for n in range(len(orders)):
         tmp = [points[m][orders[n][0]],
                points[m][orders[n][1]],
                points[m][orders[n][2]],
                points[m][orders[n][3]]]

         ## To avoid duplicates who have the same value in two directions
         ## i.e. [2,1,1,0] = [2,1,1,0]
         if tmp not in point_orders:
            point_orders.append(tmp)

      for n in range(len(point_orders)):
         
         ## This list will hold all sign permutations for the mth point
         ## in the first orthant
         point_signs = []
 
         ## For each of the possible sign permutations
         for o in range(len(signs)):
            tmp = [point_orders[n][0] * signs[o][0],
                   point_orders[n][1] * signs[o][1],
                   point_orders[n][2] * signs[o][2],
                   point_orders[n][3] * signs[o][3]]
 
            ## To avoid duplicates with points on the axes
            ## i.e. 0 * 1 = 0 * -1
            if tmp not in point_signs:
               point_signs.append(tmp)

         final = final + point_signs
            
   return final

But this is still pretty much the same complexity. I still wanted to get something faster. I did some research online and was unable to find anything. The closest I came was an article in Communications on Pure and Applied Mathematics by Rabin and Shallit^[3] that used randomized algorithms to produce a single lattice point such that $x^2+y^2+z^2+w^2=d^2$. While fast (complexity of $O(\log^2 d)$), this method produces only one lattice points and there is no way to guarantee which point you get.

Finally, I stumbled upon a post on the Computer Science Stack Exchange^[4] where a user suggested this nearly $O(d)$ algorithm. The way this algorithm works is that is populates a list that is $d+1$ long with pairs of numbers whose squares sum to the index of that particular element. So for the $0$^th index we will have the pair $[0,0]$ because $0^2+0^2=0$, for the $5$^th index we will have the pair $[2,1]$ because $2^2+1^2=5$, and for the $25$^th index we will have the pairs $[5,0]$ and $[4,3]$ because $5^2+0^2=4^2+3^2=25$. Now, we know that the pairs in the $0$^th location and the pairs in the $d$^th location will have a four-square sum of $d$ and the pairs in the $5$^th location and the pairs in the $d-5$^th location will have a four-square sum of $d$. We do a meet in the middle loop through this list and we have our basic representations!

For example, consider $d=6$. We have the following list of pairs which we can generation in $O(d)$ time. $$\begin{array}{|c|c|c|c|c|c|c|c|} \hline Index & 0 & 1 & 2 & 3 & 4 & 5 & 6 \\ \hline Pairs & [0,0] & [1,0] & [1,1] & \phantom{[0,0]} & [2,0] & [2,1] & \phantom{[0,0]}\\ \hline \end{array}$$ Now look at the pair in index $1$, $[1,0]$ and the pair in index $5$, $[2,1]$. Because $1+5=6$, we know the sum of the squares of all four of the numbers in these two pairs will sum to $6$: $2^2+1^2+1^2+0^2=6$. A wonderfully constructed list! The code for this method is below.

from math import sqrt

signs = [[ 1, 1, 1, 1], [ 1, 1, 1,-1], [ 1, 1,-1, 1], [ 1, 1,-1,-1],
         [ 1,-1, 1, 1], [ 1,-1, 1,-1], [ 1,-1,-1, 1], [ 1,-1,-1,-1],
         [-1, 1, 1, 1], [-1, 1, 1,-1], [-1, 1,-1, 1], [-1, 1,-1,-1],
         [-1,-1, 1, 1], [-1,-1, 1,-1], [-1,-1,-1, 1], [-1,-1,-1,-1]]

orders = [[0,1,2,3], [0,1,3,2], [0,2,1,3], [0,2,3,1], [0,3,1,2], [0,3,2,1],
          [1,0,2,3], [1,0,3,2], [1,2,0,3], [1,2,3,0], [1,3,0,2], [1,3,2,0],
          [2,0,1,3], [2,0,3,1], [2,1,0,3], [2,1,3,0], [2,3,0,1], [2,3,1,0],
          [3,0,1,2], [3,0,2,1], [3,1,0,2], [3,1,2,0], [3,2,0,1], [3,2,1,0]]

def sum_squares(d):
   ## Element i of sum_pairs will be populated with all pairs [a,b]
   ## such that a, b < sqrt(d) and a*a + b*b = i
   sum_pairs = [[] for i in range(d+1)]

   ## Fill sum_pairs
   for a in range(0,int(sqrt(d)) + 1):
      for b in range(0,min(a + 1, int(sqrt(d - a**2)) + 1)):
         i = a*a + b*b
         if i <= d:
            sum_pairs[i].append([a,b])
         else:
            break

   ## This list will be populated with all points [i,j,k,l]
   ## such that i >= j >= k >= l >= 0 whose squares sum to d
   points = []

   ## This list will hold the points from all orthants
   final = []
   
   for i in range(0,d//2 + 1):
      for j in range(0,len(sum_pairs[i])):
         for k in range(0,len(sum_pairs[d-i])):
            tmp = [sum_pairs[i][j][0], sum_pairs[i][j][1],
                   sum_pairs[d-i][k][0], sum_pairs[d-i][k][1]]
            tmp.sort()
            if tmp not in points:
               points.append(tmp)

   ## For each basic point found in the first orthant
   for m in range(len(points)):
      ## This list will hold all order permutations for the
      ## ith sign permutation
      point_orders = []

      ## For each of the possible order permutations
      for n in range(len(orders)):
         tmp = [points[m][orders[n][0]],
                points[m][orders[n][1]],
                points[m][orders[n][2]],
                points[m][orders[n][3]]]

         ## To avoid duplicates who have the same value in two directions
         ## i.e. [2,1,1,0] = [2,1,1,0]
         if tmp not in point_orders:
            point_orders.append(tmp)

      for n in range(len(point_orders)):
         
         ## This list will hold all sign permutations for the mth point
         ## in the first orthant
         point_signs = []
 
         ## For each of the possible sign permutations
         for o in range(len(signs)):
            tmp = [point_orders[n][0] * signs[o][0],
                   point_orders[n][1] * signs[o][1],
                   point_orders[n][2] * signs[o][2],
                   point_orders[n][3] * signs[o][3]]
 
            ## To avoid duplicates with points on the axes
            ## i.e. 0 * 1 = 0 * -1
            if tmp not in point_signs:
               point_signs.append(tmp)

         final = final + point_signs
            
   return final

Ultimately, we've found a very quick way to find all lattice points on the hypersphere of radius $\sqrt{d}$. For the methods that we had to permute (either sign or position), most of the computation time is spent permuting. Therefore, I will spend some time working on a faster way to permute these basic representations while being mindful of repeats. Those methods will be discussed in another post.
[1] https://en.wikipedia.org/wiki/Lagrange%27s_four-square_theorem
[2] https://en.wikipedia.org/wiki/Jacobi%27s_four-square_theorem
[3] M.O. Rabin, J.O. Shallit, Randomized Algorithms in Number Theory, Communications on Pure and Applied Mathematics (1986), no. S1, pp S239-S256.
[4]http://cs.stackexchange.com/questions/2988/how-fast-can-we-find-all-four-square-combinations-that-sum-to-n