Lagrange
Decomposing a non-negative integer into a sum of at most four squares
For the integer expression you can use parentheses and the usual operators + - * / % (modulo), ^ (exponentiation) and ! (factorial). In addition the following functions are predefined:
- prime(n): Computes the smallest (probable) prime p ≥ n
- sqrt(n): Computes y such that y2 ≤ n < (y + 1)2
- ip(n): Computes the n-th prime (ip(1) = 2)
- pp(n): Computes the product of the first n (probable) primes
- d(min, mod, rem): Computes the smallest odd (probable) prime p ≥ min which is congruent rem modulo mod. Make sure that mod ≥ 1 and that gcd(mod, rem) = 1.
Explanation
The following explains in more details the mathematical foundation and some implementation aspects of the applet.
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) {
assert (p is prime) and (p = 1 mod 4)
if (p = 5 mod 8) b := 2;
else {
b := 3;
// nextPrime(b) returns the smallest prime larger than b
while (b^((p - 1) / 2) = 1 mod p) b := nextPrime(b);
}
b := b^((p - 1) / 4) mod p;
// b is now an imaginary unit, i.e. b^2 = -1 mod p
a := p;
while (b^2 > p) {
a, b := b, a mod b;
}
return [b, a mod b]
}
Implementation notes
- The computation of b(p - 1) / 2 mod p can be replaced by the Jacobi symbol (b / p) which is more efficient to compute.
- If the Jacobi symbol is used, special care needs to be taken to ensure termination of the first while loop in the case that p is not a prime. This can happen as in practice p is only probable prime but not necessarily guaranteed to be prime. A simple solution is to let the Jacobi symbol (b / p) return -1 in case gcd(b, p) > 1. Even more trouble can arise when p is not only composite but even a square number. In this case, the Jacobi symbol is always 1 and finding the smallest prime factor might take a long time. Therefore it makes sense to check whether p is actually a square after too many primes have been checked.
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;
// computed by counting the bits of n
y = 2^((floor(log2(n)) + 1) / 2);
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.
- If n = 3 mod 8 we try to represent n as x2 + 2p with prime p = 1 mod 4 and write p = y2 + z2. This yields the representation n = x2 + (y + z)2 + (y - z)2.
- If n = 1 mod 4 or n = 2 mod 4 we try to represent n as x2 + p with prime p = 1 mod 4 and write p = y2 + z2. This yields the representation n = x2 + y2 + z2.
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 := x^2 + y^2 + z^2; // table lookup
return 2^v * [d, x, y, z];
}
if (n = 3 mod 8) {
find odd x such that p := (n - x^2) / 2 is prime;
// sequential search, start with largest possible x
[y, z] := decomposePrime(p);
return 2^v * [d, x, y + z, abs(y - z)];
}
find x such that p := n - x^2 is prime;
// sequential search, start with largest possible x
[y, z] := decomposePrime(p);
return 2^v * [d, x, y, z];
}
Implementation notes
- In a practical implementation the primality test will be a probabilistic one, i.e. it might fail to reject a number which is in fact composite. This is not a problem since in this case the function decomposePrime will most certainly fail to produce the correct result but this is easy to verify.
- The special cases are the following ones:
- 1 = 12
- 2 = 12 + 12
- 3 = 12 + 12 + 12
- 10 = 12 + 32
- 34 = 32 + 32 + 42
- 58 = 32 + 72
- 85 = 62 + 72
- 130 = 32 + 112
- 214 = 32 + 62 + 132
- 226 = 82 + 92 + 92
- 370 = 82 + 92 + 152
- 526 = 62 + 72 + 212
- 706 = 152 + 152 + 162
- 730 = 12 + 272
- 1414 = 62 + 172 + 332
- 1906 = 132 + 212 + 362
- 2986 = 212 + 322 + 392
- 9634 = 562 + 572 + 572