The first proof was given by Lagrange

In 1770 Lagrange has shown that any natural number can be written as a sum of at most four squares, a conjecture dating back to the days of Diophantus. In this section we show how this decomposition can effectively be computed even for large numbers.

Backtracking: simple but not efficient

The simplest solution one can think of works by subtracting a trial square from the number to be decomposed and then trying to decompose the resulting number into a sum of three squares. Using the same idea, this problem is then reduced to the representation of two squares. The final step of this algorithm asks to decompose an integer into just one square which boils down to checking whether a number D is a square. This can easily be accomplished by adopting the Newton iteration
xn+1 = (xn + D / xn) / 2
for finding the square root of D to the integers. Unfortunately, the backtracking algorithm does not perform very well in practice. Nevertheless we will use one piece of this algorithm, namely the procedure to find the integer square root x of a number D which is defined as the largest x such that x2 ≤ D.

Reducing the problem to the sum of three squares

In a first step we will reduce the problem to the representation as a sum of three squares and later to the representation of a prime as a sum of two squares. An important result is that any integer n which is not of the form 4a(8m + 7) can be written as a sum of three squares. We will therefore first check whether n, the number to be decomposed, is representable as a sum of three squares. If this is not the case then n - (2a)2 is not of the forbidden form and we will obtain a representation as a sum of four squares. Note that one could also try to subtract the largest possible square from n that makes sure that n is not of the forbidden in order to have a smaller number to be represented as a sum of three squares.

Representing a prime as the sum of two squares

Before we can tackle the problem of representation as a sum of three squares we need to solve the simpler problem of representing a prime p, p congruent 1 modulo 4, as a sum of two squares (note that numbers congruent 3 modulo 4 cannot be represented as a sum of two squares since any square is either congruent 0 modulo 4 or 1 modulo 4). One can show that any prime p congruent 1 modulo 4 is the sum of two squares. The calculation of this representation works in two steps: computing an imaginary unit and reduction to the final result. The following pseudo-code gives an efficient algorithm for decomposing a prime congruent 1 modulo 4 into a sum of two squares. It is a simplified variant of Cornacchia's algorithm but avoids computing the (integer) square root. 
integer[] decomposePrime(integer p) { // (p is prime) and (p = 1 mod 4)
   if (p = 5 mod 8) b := 2;
   else { b := 3; while (b(p - 1) / 2 = 1 mod p) b := nextPrime(b); }
   // nextPrime(b) returns the smallest prime larger than b
   b := b(p - 1) / 4 mod p;
   // b is now an imaginary unit, i.e. b2 = -1 mod p
   a := p;
   while (b2 > p) {
	a, b := b, a mod b;
   }
   return [b, a mod b]
}

Implementation notes

The Jacobi symbol is calculated by the following algorithm: 
integer Jacobi(integer a, integer b) {
  sign := 1;
  while (a > 1) {
    if (a = 0 mod 4) a := a / 4;
    else if (a = 0 mod 2) {
      if ((b = 3 mod 8) or (b = 5 mod 8)) sign := -sign;
      a := a / 2;
    }
    else {
      if ((a = 3 mod 4) and (b = 3 mod 4)) sign := -sign;
      t := a; a := b mod a; b := t;
    }
  }
  if (a = 0) return -1; // ensure termination in decomposePrime
  else return sign;
}
The integer square root which was already mentioned in the backtracking algorithm can be calculated as follows: 
integer isqrt(integer n) {
  if (n ≤ 1) return n;
  y = 2(floor(log2(n)) + 1) / 2; // easily computed by counting the bits of n
  do {
    x = y;
    y = (x + n / x) / 2; // truncating integer division
  } while (x > y)
  return x;
}
The real implementation (both in Python and Java) uses a much more efficient approach which is also a variant of the Newton iteration from above but uses in its recursive branch an idea from Paul Zimmermann (Karatsuba Square Root).

Representing a number as a sum of three squares

For this problem we use the approach from Rabin and Shallit which tries to employ the algorithm from the previous section. Let n be the number to be written as a sum of three squares. Without loss of generality we can assume that n is not divisible by four and that n is not congruent 7 mod 8. Obviously the success of this method depends on the existence of an appropriate prime p. Rabin and Shallit conjecture that for sufficiently large n an appropriate prime can always be found. Indeed, numerical experiments suggest that 9634 is the largest integer which cannot be represented as a square plus a prime congruent 1 mod 4 or as a square plus twice a prime congruent 1 mod 4.

Putting it all together

We can now put all ideas from above together to obtain the following algorithm which represents a positive integer as a sum of at most four squares. 
integer[] decompose(integer n) {
  v = 0;
  while (n = 0 mod 4) { v := v + 1; n := n / 4; }
  if (n = 7 mod 8) { d := 1; n := n - 1; } else d := 0;
  if (n is special case) { // table lookup - see cases below
    n := x2 + y2 + z2; // table lookup
    return 2v * [d, x, y, z];
  }
  if (n = 3 mod 8) {
    find odd x such that p := (n - x2) / 2 is prime;
    // sequential search, start with largest possible x
    [y, z] := decomposePrime(p);
    return 2v * [d, x, y + z, abs(y - z)];
  }
  find even x such that p := n - x2 is prime;
  // sequential search, start with largest possible x
  [y, z] := decomposePrime(p);
  return 2v * [d, x, y, z];
}

Implementation notes

home