Как найти все возможные значения четырех переменных при квадрате суммы в N?
A^2+B^2+C^2+D^2 = N Дано целое число Nраспечатать все возможные комбинации целочисленных значений ABCD которые решают уравнение.
Я предполагаю, что мы можем сделать лучше, чем грубая сила.
4 ответа
На странице Википедии есть некоторая интересная справочная информация, но теорема Лагранжа о четырех квадратах (или, точнее, теорема Бахе - Лагранж только доказал это) не дает подробных сведений о том, как найти указанные квадраты.
Как я уже сказал в своем комментарии, решение будет нетривиальным. В этой статье обсуждается разрешимость четырехугольных сумм. В документе утверждается, что:
Не существует удобного алгоритма (кроме простого, упомянутого во втором абзаце этой статьи) для поиска дополнительных решений, на которые указывает вычисление представлений, но, возможно, это упростит поиск, давая представление о том, какие виды решений и не существует.
Есть несколько других интересных фактов, связанных с этой темой. Существуют другие теоремы, которые утверждают, что каждое целое число может быть записано как сумма четырех конкретных кратных квадратов. Например, каждое целое число может быть записано как N = a^2 + 2b^2 + 4c^2 + 14d^2. Существует 54 подобных случая, которые верны для всех целых чисел, и Рамануджан предоставил полный список в 1917 году.
Для получения дополнительной информации см. Модульные формы. Это нелегко понять, если у вас нет опыта в теории чисел. Если бы вы могли обобщить 54 формы Рамануджана, вам может быть легче с этим. С учетом вышесказанного в первой статье, которую я цитирую, есть небольшой фрагмент, в котором обсуждается алгоритм, который может найти любое решение (даже если я нахожу его немного сложным для подражания):
Например, в 1911 году сообщалось, что калькулятора Готфрида Ракла попросили уменьшить N = 15663 как сумму четырех квадратов. Он получил решение 125^2 + 6^2 + 1^2 + 1^2 за 8 секунд, а затем сразу 125^2 + 5^2 + 3^2 + 2^2. Более сложная проблема (отраженная первым термином, который находится дальше от исходного числа, с соответственно более поздними терминами) заняла 56 секунд: 11399 = 105^2 + 15^2 + 8^2 + 5^2. В общем, стратегия состоит в том, чтобы установить первый член, который будет наибольшим квадратом ниже N, и попытаться представить меньший остаток как сумму трех квадратов. Затем первый член устанавливается на следующий по величине квадрат ниже N и так далее. Со временем калькулятор молний привыкнет к выражению небольших чисел в виде сумм квадратов, что ускорит процесс.
(Акцент мой.)
Алгоритм описывается как рекурсивный, но его можно легко реализовать итеративно.
Наивная грубая сила будет примерно такой:
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;
}
Это все еще грубая сила, но не такая грубая, поскольку она понимает, когда нужно останавливать каждый уровень зацикливания как можно раньше.
На моем (относительно) скромном ящике, который занимает менее секунды (а), чтобы получить все решения для чисел до 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 в этом случае требуется менее десяти секунд, а не более шести минут:
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;
}
(а): Все время сделано с внутренним printf закомментировано так, что ввод / вывод не искажает цифры.
Кажется, что все целые числа могут быть получены такой комбинацией:
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)
Это единственные два возможных кортежа - эй престо!
Небфа имеет отличный ответ. одно предложение:
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)).
теперь задача состоит в том, чтобы найти 4 числа из массива, сумма (a, b,c,d) = N;
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]; теперь вопрос состоит в том, чтобы найти 3 числа из S суммы (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),
Вот результат теста, выполняющего 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