一尘不染

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

algorithm

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

例子

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

我的方法

取平方根,其余部分重复此步骤:

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

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

3 2 + 3 2 + 2 2 + 1 2

  1. 还有其他算法可以做到吗?

  2. 总是有可能吗?


阅读 392

收藏
2020-07-28

共1个答案

一尘不染

永远有可能吗?

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

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

已经通过多种方式证明了这一点。

算法

有一些更聪明的算法,但是我建议使用以下算法:

将数字分解为主要因子。它们不必是素数,但是它们越小越好:素数是最好的。然后针对以下每个因素解决任务,并将得到的4个正方形与先前找到的4个正方形(具有Euler的四正方形标识)组合

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

  1. 给定一个数 n (上述因素之一),得到不大于 n 的最大平方,然后看一下用Legendre的三平方定理n 减去该平方是否可以写成三个平方的和:当且仅当此数字不是以下格式时,才有可能:

4 一个(8B + 7)

如果找不到合适的正方形,请尝试下一个较小的正方形,直到找到一个。它保证会有一个,并且大多数都可以在几次重试中找到。

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

如果所有与3模4一致的 n 的质数因子均出现偶数指数,则 n 可表示为两个平方之和。反之亦成立。

如果找不到合适的正方形,请尝试下一个较小的正方形,直到找到一个。保证会有一个。

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

这大概是个主意。为了找到主要因素,有几种解决方案。下面我将仅使用Eratosthenes筛

这是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(' + '));

迈克尔·巴尔(MichaelBarr)撰写的有关该主题的文章可能代表了一种更省时的方法,但是本文的目的更多是作为一种证明,而不是一种算法。但是,如果您需要更高的时间效率,则可以考虑将其与更高效的分解算法一起使用。

2020-07-28