FrBrGeorge/FourSquaresSum

Измышления о теореме Лагранжа с примерами программ

Известно, что любое натуральное число можно представить в виде суммы не более чем четырех квадратов неотрицательных целых чисел (теорема Лагранжа).

Ввести натуральное N и найти для него такие целые неотрицательные ${\displaystyle x,y,z}$ и ${\displaystyle t}$, чтобы ${\displaystyle x^2+y^2+z^2+t^2=N}$.

Очевидный способ:

   1 #include <stdio.h >
   2 
   3 void main() {
   4     int N, x, y, z, t;
   5 
   6     scanf("%d", &N);
   7     for(x=1; x*x <= N; x++)
   8         for(y=0; y*y <= N; y++)
   9             for(z=0; z*z <= N; z++)
  10                 for(t=0; t*t <= N; t++)
  11                     if(x*x+y*y+z*z+t*t == N)
  12                         printf("%d²+%d²+%d²+%d² = %d\n", x, y, z, t, N);
  13 }

Кстати, оператор if внутри циклов выполнится ${\displaystyle {(\sqrt{N}+1)}^4}$ раз (больше, чем ${\displaystyle 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=5⁴, для 10000 — 104060401

Что-то не так: слишком много одинаковых вариантов. Заметим, что (например, для N==11) это всё перестановки одного и того же сочетания. Как добиться, чтобы выводилась только одна перестановка? Из всех перестановок упорядочена (например, по возрастанию) только одна, так что можно смело рассматривать только варианты . Любой другой вариант будет перестановкой одного из этих.

   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 (уже можно дождаться конца ). Как и следовало ожидать, это число сочетаний с повторениями ${\displaystyle C_k^m}$, где ${\displaystyle m=\sqrt{N}+1}$, а ${\displaystyle 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" и другие предыдущие сочетания алгоритме.

И вот тут мы подходим к самому интересному месту рассуждений. А как вообще сделать так, чтобы не проверять в цикле ненужные предыдущие сочетания?

Очевидно, недостаточно одной проверки, а не получился ли у нас такой кандидат в решения, после которого поиск не имеет смысла. Верхнюю и нижнюю границы такого набора надо научиться вычислять (или подбирать по предыдущему набору).

Поможет в этом следующее соображение. Будем обозначать целую часть некоторого ${\displaystyle k}$ в виде ${\displaystyle \lfloor k\rfloor }$. Если ${\displaystyle N}$ — полный квадрат, то "${\displaystyle \sqrt{N}}$ 0 0 0" — правильный ответ. Если нет, то всё равно ${\displaystyle x}$ не может быть больше ${\displaystyle \sqrt{N}}$, потому что ${\displaystyle {(\lfloor \sqrt{N}\rfloor +1)}^2>N}$. А каково наименьшее возможное значение ${\displaystyle x}$? Вспомним, что в наиболее эффективном алгоритме (выдающем ровно одно размещение из нескольких возможных) x, y, z, t нестрого упорядочены по убыванию. Следовательно, если ${\displaystyle x}$ настолько мало, что , ни при каких правильного ответа не получится.

Это даёт нам нижнюю оценку для ${\displaystyle x}$: ${\displaystyle x\ge \frac{\sqrt{N}}{4}}$.

Точно так же, если , то решения, опять-таки, нет. Таким образом, нижние оценки для всех переменных:

Последнее неравенство примечательно. Если немного подумать, станет очевидно, что внутренний цикл — лишний: если x, y и z уже посчитаны, единственное t, которое надо проверить — это ${\displaystyle \lfloor \sqrt{N-x^2-y^2-z^2}\rfloor }$. Иначе говоря, ответ есть только когда ${\displaystyle N-x^2-y^2-z^2}$ — это полный квадрат, не превосходящий ${\displaystyle 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 }

Избавляться от вещественной арифметики надо не только потому, что она медленная, но и потому, что она не по-настоящему вещественная. В самом деле, не всякий алгоритм вычисления квадратного корня гарантирует, что ${\displaystyle \sqrt{100}}$ — это ${\displaystyle 10.0}$, а не ${\displaystyle 9.99999\left(9\right)\dots }$, ведь это одинаковые вещественные числа. Однако преобразование (int) вполне способно превратить 10.0 в 10, а 9.9999999 — в 9.

И всё-таки было бы неплохо подбираться к нужному значению переменной откуда-то ближе, чем от 0. Допустим, для ${\displaystyle x=k}$ нижняя граница ${\displaystyle y}$ составляет ${\displaystyle m}$. Тогда для следующего ${\displaystyle x=k+1}$ нижняя граница ${\displaystyle y}$ будет чуть меньше ${\displaystyle m}$:

• Для ${\displaystyle x=k}$: ${\displaystyle y^2=\frac{N-k^2}{3}=\frac{N}{3}-\frac{k^2}{3}}$

• Для ${\displaystyle x=k+1}$: ${\displaystyle y^2=\frac{N-{(k+1)}^2}{3}=\frac{N}{3}-\frac{k^2}{3}-2\frac{k}{3}-\frac{1}{3}=m-\frac{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 }

пришлось дополнительно отследить нулевые значения переменных на случай, когда левая часть неравенства вида ${\displaystyle 2{(minz-1)}^2\ge N-x^2-y^2}$ увеличится после уменьшения ${\displaystyle minz}$.