A^2+B^2+C^2+D^2 = N给定一个整数N，打印出ABCD解方程的所有可能的整数值组合。
我猜我们可以比蛮力做得更好。
问问题
7628 次
4 回答
7
天真的蛮力是这样的：
n = 3200724;
lim = sqrt (n) + 1;
for (a = 0; a <= lim; a++)
for (b = 0; b <= lim; b++)
for (c = 0; c <= lim; c++)
for (d = 0; d <= lim; d++)
if (a * a + b * b + c * c + d * d == n)
printf ("%d %d %d %d\n", a, b, c, d);
不幸的是，这将导致超过一万亿次循环，效率并不高。
实际上，您可以通过在每个级别打折大量的不可能性来做得更好，例如：
#include <stdio.h>
int main(int argc, char *argv[]) {
int n = atoi (argv[1]);
int a, b, c, d, na, nb, nc, nd;
int count = 0;
for (a = 0, na = n; a * a <= na; a++) {
for (b = 0, nb = na - a * a; b * b <= nb; b++) {
for (c = 0, nc = nb - b * b; c * c <= nc; c++) {
for (d = 0, nd = nc - c * c; d * d <= nd; d++) {
if (d * d == nd) {
printf ("%d %d %d %d\n", a, b, c, d);
count++;
}
tot++;
}
}
}
}
printf ("Found %d solutions\n", count);
return 0;
}
它仍然是蛮力，但不像它那么粗暴，因为它知道何时尽早停止每个级别的循环。
在我的（相对）适度的盒子上，需要不到一秒（a）的时间来获得最多 50,000 个数字的所有解决方案。除此之外，它开始需要更多时间：
n time taken
---------- ----------
100,000 3.7s
1,000,000 6m, 18.7s
因为n = ten million，在我杀死它之前，它已经运行了大约一个半小时。
所以，我想说蛮力在某种程度上是完全可以接受的。除此之外，还需要更多的数学解决方案。
为了提高效率，您只能在d >= c >= b >= a. 那是因为您可以将这些组合中的所有解决方案构建成排列（在两个或多个a、b、c或的值d相同的情况下，可能会删除重复项）。
此外，d循环体不需要检查 的每个值d，只需检查最后一个可能的值。
在这种情况下获得结果1,000,000需要不到 10 秒而不是超过 6 分钟：
0 0 0 1000
0 0 280 960
0 0 352 936
0 0 600 800
0 24 640 768
: : : :
424 512 512 544
428 460 500 596
432 440 480 624
436 476 532 548
444 468 468 604
448 464 520 560
452 452 476 604
452 484 484 572
500 500 500 500
Found 1302 solutions
real 0m9.517s
user 0m9.505s
sys 0m0.012s
该代码如下：
#include <stdio.h>
int main(int argc, char *argv[]) {
int n = atoi (argv[1]);
int a, b, c, d, na, nb, nc, nd;
int count = 0;
for (a = 0, na = n; a * a <= na; a++) {
for (b = a, nb = na - a * a; b * b <= nb; b++) {
for (c = b, nc = nb - b * b; c * c <= nc; c++) {
for (d = c, nd = nc - c * c; d * d < nd; d++);
if (d * d == nd) {
printf ("%4d %4d %4d %4d\n", a, b, c, d);
count++;
}
}
}
}
printf ("Found %d solutions\n", count);
return 0;
}
而且，根据 DSM 的建议，d循环可以完全消失（因为d（负数折扣）只有一个可能的值并且可以计算出来），这将一百万个案例减少到两秒，而十百万案例到更易于管理的 68 秒。
该版本如下：
#include <stdio.h>
#include <math.h>
int main(int argc, char *argv[]) {
int n = atoi (argv[1]);
int a, b, c, d, na, nb, nc, nd;
int count = 0;
for (a = 0, na = n; a * a <= na; a++) {
for (b = a, nb = na - a * a; b * b <= nb; b++) {
for (c = b, nc = nb - b * b; c * c <= nc; c++) {
nd = nc - c * c;
d = sqrt (nd);
if (d * d == nd) {
printf ("%d %d %d %d\n", a, b, c, d);
count++;
}
}
}
}
printf ("Found %d solutions\n", count);
return 0;
}
(a) : 所有的时序都是在内部printf注释掉的情况下完成的，这样 I/O 就不会扭曲数字。
于 2012-07-31T04:45:50.940 回答
4
维基百科页面有一些有趣的背景信息，但拉格朗日的四平方定理（或者，更准确地说，巴切特定理 - 拉格朗日只证明了这一点）并没有真正详细说明如何找到所说的平方。
正如我在评论中所说，解决方案将很重要。本文讨论了四平方和的可解性。该文件称：
没有方便的算法（除了本文第二段中提到的简单算法）来查找表示计算所指示的其他解决方案，但也许这将通过给出什么样的解决方案的想法来简化搜索不存在。
还有一些与该主题相关的其他有趣事实。还有其他定理指出每个整数都可以写成四个特定平方倍数的和。例如，每个整数都可以写成 N = a^2 + 2b^2 + 4c^2 + 14d^2。有 54 个这样的情况适用于所有整数，Ramanujan 在 1917 年提供了完整的列表。
有关详细信息，请参阅模块化表单。除非你有一些数论背景，否则这并不容易理解。如果您可以概括 Ramanujan 的 54 种形式，您可能会更轻松。话虽如此，在我引用的第一篇论文中，有一个小片段讨论了一种可能找到所有解决方案的算法（即使我发现它有点难以理解）：
例如，据报道，1911 年，计算器 Gottfried Ruckle 被要求将 N = 15663 减少为四个平方和。他在 8 秒内产生了 125^2 + 6^2 + 1^2 + 1^2 的解，紧随其后的是 125^2 + 5^2 + 3^2 + 2^2。一个更困难的问题（反映在距离原始数字较远的第一项，后面的项相应地更大）花费了 56 秒：11399 = 105^2 + 15^2 + 8^2 + 5^2。一般来说，策略是首先将第一项设置为 N 以下的最大平方，并尝试将较小的余数表示为三个平方的和。然后将第一项设置为 N 以下的下一个最大正方形，依此类推。随着时间的推移，闪电计算器将熟悉将小数表示为平方和，这将加快处理速度。
（强调我的。）
该算法被描述为递归的，但它可以很容易地迭代实现。
于 2012-07-31T03:52:04.627 回答
2
似乎所有整数都可以通过这样的组合组成：
0 = 0^2 + 0^2 + 0^2 + 0^2
1 = 1^2 + 0^2 + 0^2 + 0^2
2 = 1^2 + 1^2 + 0^2 + 0^2
3 = 1^2 + 1^2 + 1^2 + 0^2
4 = 2^2 + 0^2 + 0^2 + 0^2, 1^2 + 1^2 + 1^2 + 1^2 + 1^2
5 = 2^2 + 1^2 + 0^2 + 0^2
6 = 2^2 + 1^2 + 1^2 + 0^2
7 = 2^2 + 1^2 + 1^2 + 1^2
8 = 2^2 + 2^2 + 0^2 + 0^2
9 = 3^2 + 0^2 + 0^2 + 0^2, 2^2 + 2^2 + 1^2 + 0^2
10 = 3^2 + 1^2 + 0^2 + 0^2, 2^2 + 2^2 + 1^2 + 1^2
11 = 3^2 + 1^2 + 1^2 + 0^2
12 = 3^2 + 1^2 + 1^2 + 1^2, 2^2 + 2^2 + 2^2 + 0^2
.
.
.
等等
当我在脑海中做了一些初步的工作时，我认为只有完美的正方形才会有超过 1 个可能的解决方案。但是，在将它们列出之后，在我看来，它们没有明显的顺序。但是，我想到了一种我认为最适合这种情况的算法：
重要的是使用 4 元组 (a, b, c, d)。在任何给定的 4 元组中，它是 a^2 + b^2 + c^2 + d^2 = n 的解，我们将为自己设置一个约束，即 a 始终是 4 中最大的，b 次之，并且依此类推，例如：
a >= b >= c >= d
另请注意，a^2 不能小于 n/4，否则平方和必须小于 n。
那么算法是：
1a. Obtain floor(square_root(n)) # this is the maximum value of a - call it max_a
1b. Obtain the first value of a such that a^2 >= n/4 - call it min_a
2. For a in a range (min_a, max_a)
在这一点上，我们已经选择了一个特定的 a，现在正在研究将 a^2 的差距缩小到 n - 即 (n - a^2)
3. Repeat steps 1a through 2 to select a value of b. This time instead of finding
floor(square_root(n)) we find floor(square_root(n - a^2))
等等等等。所以整个算法看起来像：
1a. Obtain floor(square_root(n)) # this is the maximum value of a - call it max_a
1b. Obtain the first value of a such that a^2 >= n/4 - call it min_a
2. For a in a range (min_a, max_a)
3a. Obtain floor(square_root(n - a^2))
3b. Obtain the first value of b such that b^2 >= (n - a^2)/3
4. For b in a range (min_b, max_b)
5a. Obtain floor(square_root(n - a^2 - b^2))
5b. Obtain the first value of b such that b^2 >= (n - a^2 - b^2)/2
6. For c in a range (min_c, max_c)
7. We now look at (n - a^2 - b^2 - c^2). If its square root is an integer, this is d.
Otherwise, this tuple will not form a solution
在步骤 3b 和 5b，我使用 (n - a^2)/3, (n - a^2 - b^2)/2。我们分别除以 3 或 2，因为元组中的值的数量尚未“固定”。
一个例子：
在 n = 12 上执行此操作：
1a. max_a = 3
1b. min_a = 2
2. for a in range(2, 3):
use a = 2
3a. we now look at (12 - 2^2) = 8
max_b = 2
3b. min_b = 2
4. b must be 2
5a. we now look at (12 - 2^2 - 2^2) = 4
max_c = 2
5b. min_c = 2
6. c must be 2
7. (n - a^2 - b^2 - c^2) = 0, hence d = 0
so a possible tuple is (2, 2, 2, 0)
2. use a = 3
3a. we now look at (12 - 3^2) = 3
max_b = 1
3b. min_b = 1
4. b must be 1
5a. we now look at (12 - 3^2 - 1^2) = 2
max_c = 1
5b. min_c = 1
6. c must be 1
7. (n - a^2 - b^2 - c^2) = 1, hence d = 1
so a possible tuple is (3, 1, 1, 1)
这是仅有的两个可能的元组 - 嘿！
于 2012-07-31T03:53:17.213 回答
0
nebffa 有一个很好的答案。一个建议：
step 3a: max_b = min(a, floor(square_root(n - a^2))) // since b <= a
max_c 和 max_d 也可以以相同的方式进行改进。
这是另一个尝试：
1. generate array S: {0, 1, 2^2, 3^2,.... nr^2} where nr = floor(square_root(N)).
现在的问题是从 sum(a, b,c,d) = N 的数组中找到 4 个数字；
2. according to neffa's post (step 1a & 1b), a (which is the largest among all 4 numbers) is between [nr/2 .. nr].
我们可以将 a 从 nr 循环到 nr/2 并计算 r = N - S[a]; 现在的问题是从 S 中找到 3 个数字 sum(b,c,d) = r = N -S[a];
这是代码：
nr = square_root(N);
S = {0, 1, 2^2, 3^2, 4^2,.... nr^2};
for (a = nr; a >= nr/2; a--)
{
r = N - S[a];
// it is now a 3SUM problem
for(b = a; b >= 0; b--)
{
r1 = r - S[b];
if (r1 < 0)
continue;
if (r1 > N/2) // because (a^2 + b^2) >= (c^2 + d^2)
break;
for (c = 0, d = b; c <= d;)
{
sum = S[c] + S[d];
if (sum == r1)
{
print a, b, c, d;
c++; d--;
}
else if (sum < r1)
c++;
else
d--;
}
}
}
运行时是O(sqare_root(N)^3).
这是在我的 VM 上运行 java 的测试结果（时间以毫秒为单位，结果# 是有效组合的总数，时间 1 有打印输出，时间 2 没有打印输出）：
N result# time1 time2
----------- -------- -------- -----------
1,000,000 1302 859 281
10,000,000 6262 16109 7938
100,000,000 30912 442469 344359
于 2012-07-31T07:49:17.973 回答