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);
残念ながら、これは1兆を超えるループになり、過度に効率的ではありません。
次のような方法で、各レベルで不可能な数を大幅に割り引くことで、実際にはそれよりも大幅に優れた成果を上げることができます。
#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;
}
それはまだブルートフォースですが、ループの各レベルをできるだけ早く停止するタイミングを理解している限り、それほどブルートではありません。
私の（比較的）控えめなボックスでは、50,000までの数のすべてのソリューションを取得するのに1秒もかかりません（a） 。それを超えて、それはより多くの時間を要し始めます：
n time taken
---------- ----------
100,000 3.7s
1,000,000 6m, 18.7s
なぜならn = ten million、私がそれを殺す前に、それは約1時間半続いていました。
ですから、ブルートフォースはある程度までは完全に許容できると思います。それを超えて、より多くの数学的解決策が必要になるでしょう。
さらに効率を上げるために、これらのソリューションをチェックできるのはd >= c >= b >= a。これは、これらの組み合わせからすべてのソリューションを順列に構築できるためです（ 、、、、またはの2つ以上の値が同一である場合、重複が削除される可能性がありますa）。bcd
さらに、dループの本体は、のすべての値をチェックする必要はなくd、最後の可能な値だけをチェックします。
その場合の結果の取得には、1,000,0006分以上ではなく10秒未満かかります。
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（（負の数を割り引く）の可能な値は1つだけであり、計算できるため）。これにより、100万ケースが2秒に短縮され、10百万ケースからはるかに扱いやすい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）：I / Oが数字を歪めないように、すべてのタイミングは内側をprintfコメントアウトして行われます。
于 2012-07-31T04:45:50.940 に答える
4
ウィキペディアのページにはいくつかの興味深い背景情報がありますが、ラグランジュの4平方定理（より正確には、バチェットの定理-ラグランジュはそれを証明しただけです）は、その正方形を見つける方法について実際には詳しく説明していません。
コメントで述べたように、解決策は簡単ではありません。この論文では、4乗和の可解性について説明します。論文は次のように主張している：
表現の計算によって示される追加のソリューションを見つけるための便利なアルゴリズム（このペーパーの2番目の段落で述べた単純なアルゴリズム以外）はありませんが、おそらくこれにより、ソリューションの種類と存在しない。
このトピックに関連する他のいくつかの興味深い事実があります。すべての整数は、4つの特定の平方の倍数の合計として記述できるという他の定理があります。たとえば、すべての整数はN = a ^ 2 + 2b ^ 2 + 4c ^ 2 + 14d^2と書くことができます。すべての整数に当てはまるこのようなケースは54あり、ラマヌジャンは1917年に完全なリストを提供しました。
詳細については、モジュラー形式を参照してください。数論のバックグラウンドがなければ、これを理解するのは簡単ではありません。ラマヌジャンの54の形式を一般化できれば、これでもっと楽になるかもしれません。そうは言っても、私が引用する最初の論文には、すべての解決策を見つけることができるアルゴリズムについて説明している小さなスニペットがあります（私はそれを理解するのは少し難しいと思いますが）：
たとえば、1911年に、電卓GottfriedRuckleがN=15663を4つの正方形の合計として減らすように求められたことが報告されました。彼は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の下の最大の正方形に設定することから始め、小さい余りを3つの正方形の合計として表すことを試みます。次に、最初の項は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）を使用することです。a ^ 2 + b ^ 2 + c ^ 2 + d ^ 2 = nの解である任意の4タプルでは、aが常に4の中で最大であり、bが次にあり、などなど：
a >= b >= c >= d
また、a^2はn/4未満にすることはできません。そうでない場合、2乗の合計は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)
これらは2つの可能なタプルだけです-ちょっとプレスト！
于 2012-07-31T03:53:17.213 に答える
0
nebffaには素晴らしい答えがあります。1つの提案：
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を実行したテスト結果です（ミリ秒単位の時間、result＃は有効な組み合わせの合計数、印刷出力ありの時間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 に答える