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

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

示例

       120 = 82 + 62 + 42 + 22
       6 = 02 + 12 + 12 + 22
       20 = 42 + 22 + 02+ 02

我的方法

对该数字求平方根,然后重复此步骤处理余数:

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

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

       32 + 32+ 22 + 12

问题

  1. 是否有其他算法可以解决这个问题?

  2. 这个问题总是有解吗?


请参见 http://cs.stackexchange.com/q/68501/755,https://dev59.com/PZ3ha4cB1Zd3GeqPcv6l,http://mathoverflow.net/q/259152/37212,http://math.stackexchange.com/q/366673/14578,http://math.stackexchange.com/q/483101/14578,http://mathoverflow.net/q/110239/37212。 - D.W.
5个回答

12

###始终可能吗?

是的,拉格朗日四平方和定理表明:

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

已经有几种证明方法。

###算法

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

将数字分解为质因数。它们不必是质数,但越小越好:所以质数是最好的。然后按下面的方式为每个因子解决任务,并使用欧拉四平方恒等式将任何结果为4个平方的值与之前找到的4个平方的值相结合。

(a2 + b2 + c2 + d2) (A2 + B2 + C2 + D2) =
               (aA + bB + cC + dD)2 +
               (aB − bA + cD − dC)2 +
               (aC − bD − cA + dB)2 +
               (aD + bC − cB − dA)2

  1. 给定一个数 n(上面提到的因素之一),找出不大于 n 的最大平方数,检查 n 减去这个平方数是否可被表示为三个平方数的和,使用勒让德三平方定理:当且仅当这个数不是以下形式时,它是可能的:

        4a(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. 尝试像第一步一样寻找实际的二次方项,但现在使用Fermat's theorem on sums of two squares来测试其可行性,其扩展意味着:
如果所有与3模4同余的质因数都出现了偶次幂,则可以表示为两个平方数的和。逆命题也成立。
3. 如果找不到适合的二次方项,请尝试下一个更小的...直到找到一个。保证会有一个。
4. 现在我们用完两个平方数后还剩下一个余数。尝试减去第三个平方数,直到得到另一个平方数,这意味着我们有一个解。这一步可以通过首先分解最大平方数因子来改进。然后,在确定了两个平方项后,每个平方项都可以再乘以该平方因子的平方根。
这大致是个想法。关于找质因数,有几种解决方案。下面我将只使用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 + 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+'²' ).join(' + '));

Michael Barr关于这个主题的文章可能代表了一种更加节省时间的方法,但这篇文章更多地是作为证明而不是算法。然而,如果您需要更高的时间效率,可以考虑使用更有效的分解算法。


我认为你的筛法代码有缺陷。它会留下很多偶数。 - Minh Nghĩa
@MinhNghĩa,我无法重现出会生成“许多偶数”的情况。例如,getPrimesUntil(100).primes 返回 [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]。你调用 getPrimesUntil 函数时使用了哪个参数,以及在结果的 primes 属性中找到了哪个非质数? - trincot
我使用 getPrimesUntil(1562),它返回 ... 31 37 44 46 ...。很奇怪,当数字改变时,质数也会改变。我注意到的第一个问题是循环没有排除所有偶数。它从 var m = 3 开始。 - Minh Nghĩa
感谢报告这个错误,@MinhNghĩa。现在应该已经被纠正了。 - trincot

3

有没有免费的链接可以找到上述算法的实现细节? - Narendra Modi

1
这里是一个解决方案,简单的 4 循环。
max = square_root(N)
for(int i=0;i<=max;i++)
  for(int j=0;j<=max;j++)
     for(int k=0;k<=max;k++)
          for(int l=0;l<=max;l++)
              if(i*i+j*j+k*k+l*l==N){
                        found
                       break;
                }

所以您可以测试任何数字。如果总和超过了,您可以在两个循环后使用中断条件来中断它。


如果你要使用线性搜索,至少更好地限制内部循环,例如 for (int j=i; j <= sqrt(N-i*i); ++j) 等等。而且最内层循环根本不需要是一个循环 - 只需测试余数是否为完全平方数即可。 - Toby Speight
这是一个有效的解决方案。 - Tatyana Molchanova

0
const fourSquares = (n) => {
    const result = [];
        for (let i = 0; i <= n; i++) {
          for (let j = 0; j <= n; j++) {
              for (let k = 0; k <= n; k++) {
                  for (let l = 0; l <= n; l++) {
                      if (i * i + j * j + k * k + l * l === n) {
                          result.push(i, j, k, l);
                          return result;
                      }
                  }
              }
          }
      }
    return result;
  }

运行时间太长了

const fourSquares = (n) => {
    const result = [];
        for (let i = 0; i <= n; i++) {
          for (let j = 0; j <= (n - i * i); j++) {
              for (let k = 0; k <= (n - i * i - j * j); k++) {
                  for (let l = 0; l <= (n - i * i - j * j - k * k); l++) {
                      if (i * i + j * j + k * k + l * l === n) {
                          result.push(i, j, k, l);
                          return result;
                      }
                  }
              }
          }
      }
    return result;
  }

0
const fourSquares = (n) => {
const result = [];
    for (let i = 0; i * i <= n; i++) {
        for (let j = 0; j * j <= n; j++) {
            for (let k = 0; k * k <= n; k++) {
                for (let l = 0; l * l <= n; l++) {
                    if (i * i + j * j + k * k + l * l === n) {
                        result.push(i, j, k, l);
                        return result;
                    }
                }
            }
        }
    }
return result;

}

const fourSquares = (n) => {
let a = Math.sqrt(n);
let b = Math.sqrt(n - a * a);
let c = Math.sqrt(n - a * a - b * b);
let d = Math.sqrt(n - a * a - b * b - c * c);
if (n === a * a + b * b + c * c + d * d) {
    return [a, b, c, d];
}

}


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接