О теореме Лагранжа
Статья, посвящённая поиску решений теоремы Лагранжа о сумме четырёх квадратов с примерами на Си.
Небольшая статья, в которой последовательно строится поиск решения теоремы Лагранжа о том, что любое натуральное число можно представить в виде суммы квадратов четырёх целых чисел. Статья сопровождается примерами программ на языке Си. Последний пример не только весьма эффективен, но и полностью базируется на целочисленной арифметике.
Измышления о теореме Лагранжа с примерами программ
Известно, что любое натуральное число можно представить в виде суммы не более чем четырех квадратов неотрицательных целых чисел (теорема Лагранжа).
Ввести натуральное N и найти для него такие целые неотрицательные $$x,y,z$$ и $$t$$, чтобы $$x^2+y^2+z^2+t^2=N$$.
Очевидный способ:
Кстати, оператор if внутри циклов выполнится $$(sqrt(N)+1)^4$$ раз (больше, чем $$N^2$$). Проверим это:
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 long int loops = 0;
6
7 scanf("%d", &N);
8 for(x=1; x*x <= N; x++)
9 for(y=0; y*y <= N; y++)
10 for(z=0; z*z <= N; z++)
11 for(t=0; t*t <= N; t++) {
12 loops++;
13 if(x*x+y*y+z*z+t*t == N)
14 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
15 }
16 printf("%ld\n", loops);
17 }
Для N=20 получим число 625=54, для 10000 — 104060401
Что-то не так: слишком много одинаковых вариантов. Заметим, что (например, для N==11) это всё перестановки одного и того же сочетания. Как добиться, чтобы выводилась только одна перестановка? Из всех перестановок упорядочена (например, по возрастанию) только одна, так что можно смело рассматривать только варианты $$x<=y<=z<=t$$. Любой другой вариант будет перестановкой одного из этих.
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 long int loops = 0;
6
7 scanf("%d", &N);
8 for(x=1; x*x <= N; x++)
9 for(y=0; y <= x; y++)
10 for(z=0; z <= y; z++)
11 for(t=0; t <= z; t++) {
12 loops++;
13 if(x*x+y*y+z*z+t*t == N)
14 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
15 }
16 printf("%ld\n",loops);
17 }
В этом алгоритме сложность пониже: для 20 — 70, для 10000 — 4598126, для 100000 — 428761520 (уже можно дождаться конца 
). Как и следовало ожидать, это число сочетаний с повторениями $$C_k^m$$, где $$m=sqrt(N)+1$$, а $$k=4$$.
На самом деле мы всё ещё «решаем не ту задачу», потому что находим все решения, а не одно — первое попавшееся. Поскольку в Си нет удобного способа выйти сразу из всех циклов после нахождения первого же сочетания, проще всего оформить весь алгоритм в виде функции. Другой вариант — ввести дополнительное условие в цикл.
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t, todo;
5
6 scanf("%d", &N);
7 todo = 1;
8 for(x=1; x*x <= N && todo; x++)
9 for(y=0; y <= x && todo; y++)
10 for(z=0; z <= y && todo; z++)
11 for(t=0; t <= z && todo; t++)
12 todo = x*x+y*y+z*z+t*t != N;
13 x--; y--; z--; t--;
14 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
15 }
Обратите внимание на то, что
условие выхода формируется и записывается в переменную todo, а сам условный оператор оказался не нужен;
- все слагаемые увеличились на 1 перед проверкой условия, так что эту единицу пришлось из них вычесть.
В дальнейшем мы для простоты переключимся на задачу поиска всех решений: метод превращения её решения в решение задачи «поиск первого» нам уже известен.
Между тем алгоритм всё ещё остаётся сильно избыточным. Чтобы в этом убедиться, добавим вывод x,y,z и t в предпоследнюю программу.
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 long int loops = 0;
6
7 scanf("%d", &N);
8 for(x=1; x*x <= N; x++)
9 for(y=0; y <= x; y++)
10 for(z=0; z <= y; z++)
11 for(t=0; t <= z; t++) {
12 loops++;
13 printf("%2d %2d %2d %2d\n", x, y, z, t);
14 if(x*x+y*y+z*z+t*t == N)
15 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
16 }
17 printf("%ld\n", loops);
18 }
Вот что получится для числа 10:
10 0 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1 2 0 0 0 2 1 0 0 2 1 1 0 2 1 1 1 2 2 0 0 2 2 1 0 2 2 1 1 2²+2²+1²+1²=10 2 2 2 0 2 2 2 1 2 2 2 2 3 0 0 0 3 1 0 0 3²+1²+0²+0²=10 3 1 1 0 3 1 1 1 3 2 0 0 3 2 1 0 3 2 1 1 3 2 2 0 3 2 2 1 3 2 2 2 3 3 0 0 3 3 1 0 3 3 1 1 3 3 2 0 3 3 2 1 3 3 2 2 3 3 3 0 3 3 3 1 3 3 3 2 3 3 3 3 35
Изрядная часть строк в начале, и ещё большая — в конце — представляют собой набор чисел, которые вообще не могут быть решением. В самом деле, если "3 1 0 0" — решение, то все последующие сочетания ("3 1 1 0", "3 1 1 1", "3 2 0 0" и т. д.) проверять вообще не надо: сумма их квадратов заведомо больше.
Разберёмся с этим путём отсечения заведомо неподходящих частей цикла (пока только некоторых).
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 long int loops = 0;
6
7 scanf("%d", &N);
8 for(x=1; x*x <= N; x++)
9 for(y=0; y <= x && x*x+y*y <= N; y++)
10 for(z=0; z <= y && x*x+y*y+z*z <= N; z++)
11 for(t=0; t <= z && x*x+y*y+z*z+t*t <= N; t++) {
12 loops++;
13 if(x*x+y*y+z*z+t*t == N)
14 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
15 }
16 printf("%ld\n", loops);
17 }
Оценка сложности улучшилась, хотя и не слишком радикально (20: 70 → 32, 10000: 4598126 → 1426253, 100000: 428761520 → 132866304, т. е. при больших N раза в три).
Дело в том, что мы всё ещё рассматриваем «ненужные» варианты слагаемых — на этот раз только те, чьи суммы квадратов заведомо меньше N. Рассмотрим наш пример для числа 10. Если "2 1 1 1" не является решением в силу того, что сумма квадратов этих чисел меньше 10, то не являются решениями и "2 1 1 0", "2 1 0 0" и другие предыдущие сочетания алгоритме.
И вот тут мы подходим к самому интересному месту рассуждений. А как вообще сделать так, чтобы не проверять в цикле ненужные предыдущие сочетания?
Очевидно, недостаточно одной проверки, а не получился ли у нас такой кандидат в решения, после которого поиск не имеет смысла. Верхнюю и нижнюю границы такого набора надо научиться вычислять (или подбирать по предыдущему набору).
Поможет в этом следующее соображение. Будем обозначать целую часть некоторого $$k$$ в виде $$|__k__|$$. Если $$N$$ — полный квадрат, то "$$sqrt(N)$$ 0 0 0" — правильный ответ. Если нет, то всё равно $$x$$ не может быть больше $$sqrt(N)$$, потому что $$(|__sqrt(N)__|+1)^2>N$$. А каково наименьшее возможное значение $$x$$? Вспомним, что в наиболее эффективном алгоритме (выдающем ровно одно размещение из нескольких возможных) x, y, z, t нестрого упорядочены по убыванию. Следовательно, если $$x$$ настолько мало, что $$4x^2 < N$$, ни при каких $$y,z,t <= x$$ правильного ответа не получится.
Это даёт нам нижнюю оценку для $$x$$: $$x>=sqrt(N)/4$$.
Точно так же, если $$3y^2 < N-x^2$$, то решения, опять-таки, нет. Таким образом, нижние оценки для всех переменных:
$$(sqrt(N))/4<=x<=sqrt(N)$$
$$(sqrt(N-x^2))/3<=y<=min(x,sqrt(N-x^2))$$
$$(sqrt(N-x^2-y^2))/2<=z<=min(y,sqrt(N-x^2-y^2))$$
$$sqrt(N-x^2-y^2-z^2)<=t<=min(z,sqrt(N-x^2-y^2-z^2))$$
Последнее неравенство примечательно. Если немного подумать, станет очевидно, что внутренний цикл — лишний: если x, y и z уже посчитаны, единственное t, которое надо проверить — это $$|__sqrt(N-x^2-y^2-z^2)__|$$. Иначе говоря, ответ есть только когда $$N-x^2-y^2-z^2$$ — это полный квадрат, не превосходящий $$z^2$$ (чтобы не нарушать порядок).
Получившееся решение (с верхними и нижними границами слагаемых) выглядит так:
1 #include <stdio.h>
2 #include <math.h>
3
4 void main() {
5 int N, x, y, z, t;
6 long int loops = 0;
7
8 scanf("%d", &N);
9 for(x=(int)(sqrt(N)/4); x*x <= N; x++)
10 for(y=(int)(sqrt((N-x*x))/3); y <= x && x*x+y*y <= N; y++)
11 for(z=(int)(sqrt((N-x*x-y*y))/2); z <= y && x*x+y*y+z*z <= N; z++)
12 if(N-x*x-y*y-z*z <= z*z) {
13 loops++;
14 t=(int)(sqrt(N-x*x-y*y-z*z));
15 if(x*x+y*y+z*z+t*t == N)
16 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
17 }
18 printf("%ld\n", loops);
19 }
Хотелось бы поточнее оценить сложность этого алгоритма, потому что она невероятно понизилась! В самом деле:
- 20: 70 → 32 → 3
- 10000: 4598126 → 1426253 → 7442
- 100000: 428761520 → 132866304 → 229493
1000000: 7194448 - 10000000: 226897862
За такое можно простить использование плавающей арифметики и функции квадратного корня. Тем более что корень нужен только для упрощения записи алгоритма, к нему можно «подобраться» перед началом цикла целочисленно. Правда, это займёт какое-то время, по-хорошему, стоит подсчитать, сколько мы «подбираемся» от 0 до минимально допустимого значения, но т. к. квадрат слагаемого сравнивается с N, в деле явно замешаны логарифмы, и грубая оценка не изменится:
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 long int loops = 0;
6
7 scanf("%d", &N);
8 for(x=1; (x+1)*(x+1) < N/4; x++);
9 for(; x*x <= N; x++) {
10 for(y=0; (y+1)*(y+1) < (N-x*x)/3; y++);
11 for(; y <= x && x*x+y*y <= N; y++) {
12 for(z=0; (z+1)*(z+1) < (N-x*x-y*y)/2; z++);
13 for(; z<= y && x*x+y*y+z*z <= N; z++)
14 if(N-x*x-y*y-z*z <= z*z) {
15 loops++;
16 for(t=0; (t+1)*(t+1) < (N-x*x-y*y-z*z); t++);
17 if(x*x+y*y+z*z+t*t == N)
18 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
19 }
20 }
21 }
22 printf("%ld\n", loops);
23 }
Избавляться от вещественной арифметики надо не только потому, что она медленная, но и потому, что она не по-настоящему вещественная. В самом деле, не всякий алгоритм вычисления квадратного корня гарантирует, что $$sqrt(100)$$ — это $$10.0$$, а не $$9.99999(9)…$$, ведь это одинаковые вещественные числа. Однако преобразование (int) вполне способно превратить 10.0 в 10, а 9.9999999 — в 9.
И всё-таки было бы неплохо подбираться к нужному значению переменной откуда-то ближе, чем от 0. Допустим, для $$x=k$$ нижняя граница $$y$$ составляет $$m$$. Тогда для следующего $$x=k+1$$ нижняя граница $$y$$ будет чуть меньше $$m$$:
Для $$x=k$$: $$y^2=(N-k^2)/3=N/3-k^2/3$$
Для $$x=k+1$$: $$y^2=(N-(k+1)^2)/3=N/3-k^2/3-2k/3-1/3=m-(2k+1)/3$$
Осталось придумать, откуда начинать подбор внутренних переменных y, z и t перед самым первым запуском цикла, а также когда значение внешних переменных (y и z соответственно) скачкообразно меняется на следующем витке более внешнего цикла (по x и по y).
Поскольку подбор теперь идёт в убывающую сторону, можно начать со значения соответствующей внешней переменной:
1 #include <stdio.h>
2
3 void main() {
4 int N, x, y, z, t;
5 int miny, minz;
6 long int loops = 0;
7
8 scanf("%d", &N);
9 for(x=1; x*x*4 < N; x++);
10 for(miny=x; x*x <= N; x++) {
11 while(miny && (miny-1)*(miny-1)*3 >= N-x*x) miny--;
12 for(y=minz=miny; y <= x && x*x+y*y <= N; y++) {
13 while(minz && (minz-1)*(minz-1)*2 >= N-x*x-y*y) minz--;
14 for(z=t=minz; z <= y && x*x+y*y+z*z <= N; z++) {
15 while(t*t > N-x*x-y*y-z*z) t--;
16 loops++;
17 if(x*x+y*y+z*z+t*t == N)
18 printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
19 }
20 }
21 }
22 printf("%ld\n", loops);
23 }
пришлось дополнительно отследить нулевые значения переменных на случай, когда левая часть неравенства вида $$2(minz-1)^2>=N-x^2-y^2$$ увеличится после уменьшения $$minz$$.