将给定数字表示为四个方格的总和

Question

我正在寻找一种算法,将一个给定的数字表示为(最多)四个方格的总和.

例子

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

我的方法

取平方根并重复重复此余数:

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

但是当N为23时,即使有解决方案,这也会失败:

       3 ² + 3 ² + 2 ² + 1 ²

题

1. 有没有其他算法可以做到这一点？

2. 它总是可能吗？

Answer 1

总是可能吗？

是的,拉格朗日的四方定理指出:

每个自然数可以表示为四个整数平方的总和.

它已经在几个方面得到了证明.

算法

有一些更智能的算法,但我会建议以下算法:

将数字分解为素数因子.它们不一定是素数,但它们越小越好:所以素数最好.然后解决下面这些因素中的每一个的任务,并将任何得到的4个方格与先前找到的4个方格与欧拉的四方方案同一性组合.

         (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. 给定一个数n(上面提到的因素之一),得到不大于n的最大正方形,并看看n减去这个正方形是否可以使用勒让德的三平方定理写成三个正方形的总和:它是可能,当且仅当此数字不是以下形式时:

           4 ^a(8b + 7)

   如果找不到合适的方块,请尝试下一个较小的方块,直到找到一个.它保证会有一个,大多数都会在几次重试中找到.

2. 尝试以与步骤1中相同的方式找到实际的第二个平方项,但现在使用费马定理测试其在两个平方的总和上的可行性,其在扩展中意味着:

   如果所有的素因子Ñ全等3模4发生于偶数指数,然后Ñ可表示为两个平方之和.反过来也有.

   如果找不到合适的方块,请尝试下一个较小的方块,直到找到一个.它保证会有一个.

3. 现在我们在减去两个方格后有一个余数.尝试减去第三个方格,直到产生另一个方格,这意味着我们有一个解决方案.通过首先分解最大的平方除数可以改善该步骤.然后,当识别出两个平方项时,每个平方项可以再次乘以该平方除数的平方根.

这大致就是这个想法.为了找到素因子,有几种解决方案.下面我将使用Eratosthenes的Sieve.

这是JavaScript代码,因此您可以立即运行它 - 它将生成一个随机数作为输入并将其显示为四个方块的总和:

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; 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+'²' ).join(' + '));

Michael Barr关于这个主题的文章可能代表了一种更具时间效率的方法,但该文本更多的是作为证据而不是算法.但是,如果您需要更高的时间效率,可以考虑使用更高效的分解算法.