algorithm - Express a given number as a sum of four squares

Question

I am looking for an algorithm that expresses a given number as a sum of (up to) four squares.

Examples

???????120 = 8² + 6² + 4² + 2²
???????6 = 0² + 1² + 1² + 2²
???????20 = 4² + 2² + 0²+ 0²

My approach

Take the square root and repeat this repeatedly for the remainder:

while (count != 4) {
    root = (int) Math.sqrt(N)
    N -= root * root
    count++
} 

But this fails when N is 23, even though there is a solution:

???????3² + 3²+ 2² + 1²

Question

1. Is there any other algorithm to do that?

2. Is it always possible?

Answer 1

###Always possible?

Yes, the Lagrange's four square theorem states that:

every natural number can be represented as the sum of four integer squares.

It has been proved in several ways.

###Algorithm

There are some smarter algorithms, but I would suggest the following algorithm:

Factorise the number into prime factors. They don't have to be prime, but the smaller they are, the better: so primes are best. Then solve the task for each of these factors as below, and combine any resulting 4 squares with the previously found 4 squares with the Euler's four-square identity.

?????????(a² + b² + c² + d²) (A² + B² + C² + D²) =
???????????????(aA + bB + cC + dD)² +
???????????????(aB ? bA + cD ? dC)² +
???????????????(aC ? bD ? cA + dB)² +
???????????????(aD + bC ? cB ? dA)²

1. Given a number n (one of the factors mentioned above), get the greatest square that is not greater than n, and see if n minus this square can be written as the sum of three squares using the Legendre's three-square theorem: it is possible, if and only when this number is NOT of the following form:

????????4^a(8b+7)

If this square is not found suitable, try the next smaller one, ... until you find one. It guaranteed there will be one, and most are found within a few retries.

2. Try to find an actual second square term in the same way as in step 1, but now test its viability using Fermat's theorem on sums of two squares which in extension means that:

   if all the prime factors of n congruent to 3 modulo 4 occur to an even exponent, then n is expressible as a sum of two squares. The converse also holds.

If this square is not found suitable, try the next smaller one, ... until you find one. It's guaranteed there will be one.

3. Now we have a remainder after subtracting two squares. Try subtracting a third square until that yields another square, which means we have a solution. This step can be improved by first factoring out the largest square divisor. Then when the two square terms are identified, each can then be multiplied again by the square root of that square divisor.

This is roughly the idea. For finding prime factors there are several solutions. Below I will just use the Sieve of Eratosthenes.

This is JavaScript code, so you can run it immediately -- it will produce a random number as input and display it as the sum of four squares:

function divisor(n, factor) {
    var divisor = 1;
    while (n % factor == 0) {
        n = n / factor;
        divisor = divisor * factor;
    }
    return divisor;
}

function getPrimesUntil(n) {
    // Prime sieve algorithm
    var range = Math.floor(Math.sqrt(n)) + 1;
    var isPrime = Array(n).fill(1);
    var primes = [2];
    for (var m = 3; m < range; m += 2) {
        if (isPrime[m]) {
            primes.push(m);
            for (var k = m * m; k <= n; k += m) {
                isPrime[k] = 0;
            }
        }
    }
    for (var m = range + 1 - (range % 2); m <= n; m += 2) {
        if (isPrime[m]) primes.push(m);
    }
    return {
        primes: primes,
        factorize: function (n) {
            var p, count, primeFactors;
            // Trial division algorithm
            if (n < 2) return [];
            primeFactors = [];
            for (p of this.primes) {
                count = 0;
                while (n % p == 0) {
                    count++;
                    n /= p;
                }
                if (count) primeFactors.push({value: p, count: count});
            }
            if (n > 1) {
                primeFactors.push({value: n, count: 1});
            }
            return primeFactors;
        }
    }
}

function squareTerms4(n) {
    var n1, n2, n3, n4, sq, sq1, sq2, sq3, sq4, primes, factors, f, f3, factors3, ok,
        res1, res2, res3, res4;
    primes = getPrimesUntil(n);
    factors = primes.factorize(n);
    res1 = n > 0 ? 1 : 0;
    res2 = res3 = res4 = 0;
    for (f of factors) { // For each of the factors:
        n1 = f.value;
        // 1. Find a suitable first square
        for (sq1 = Math.floor(Math.sqrt(n1)); sq1>0; sq1--) {
            n2 = n1 - sq1*sq1;
            // A number can be written as a sum of three squares
            // <==> it is NOT of the form 4^a(8b+7)
            if ( (n2 / divisor(n2, 4)) % 8 !== 7 ) break; // found a possibility
        }
        // 2. Find a suitable second square
        for (sq2 = Math.floor(Math.sqrt(n2)); sq2>0; sq2--) {
            n3 = n2 - sq2*sq2;
            // A number can be written as a sum of two squares
            // <==> all its prime factors of the form 4a+3 have an even exponent
            factors3 = primes.factorize(n3);
            ok = true;
            for (f3 of factors3) {
                ok = (f3.value % 4 != 3) || (f3.count % 2 == 0);
                if (!ok) break;
            }
            if (ok) break;
        }
        // To save time: extract the largest square divisor from the previous factorisation:
        sq = 1;
        for (f3 of factors3) {
            sq *= Math.pow(f3.value, (f3.count - f3.count % 2) / 2); 
            f3.count = f3.count % 2;
        }
        n3 /= sq*sq;
        // 3. Find a suitable third square
        sq4 = 0;
        // b. Find square for the remaining value:
        for (sq3 = Math.floor(Math.sqrt(n3)); sq3>0; sq3--) {
            n4 = n3 - sq3*sq3;
            // See if this yields a sum of two squares:
            sq4 = Math.floor(Math.sqrt(n4));
            if (n4 == sq4*sq4) break; // YES!
        }
        // Incorporate the square divisor back into the step-3 result:
        sq3 *= sq;
        sq4 *= sq;
        // 4. Merge this quadruple of squares with any previous 
        // quadruple we had, using the Euler square identity:
        while (f.count--) {
            [res1, res2, res3, res4] = [
                Math.abs(res1*sq1 + res2*sq2 + res3*sq3 + res4*sq4),
                Math.abs(res1*sq2 - res2*sq1 + res3*sq4 - res4*sq3),
                Math.abs(res1*sq3 - res2*sq4 - res3*sq1 + res4*sq2),
                Math.abs(res1*sq4 + res2*sq3 - res3*sq2 - res4*sq1)
            ];
        }
    }
    // Return the 4 squares in descending order (for convenience):
    return [res1, res2, res3, res4].sort( (a,b) => b-a );
}

// Produce the result for some random input number
var n = Math.floor(Math.random() * 1000000);
var solution = squareTerms4(n);
// Perform the sum of squares to see it is correct:
var check = solution.reduce( (a,b) => a+b*b, 0 );
if (check !== n) throw "FAILURE: difference " + n + " - " + check;
// Print the result
console.log(n + ' = ' + solution.map( x => x+'2' ).join(' + '));