О теореме Лагранжа

Статья, посвящённая поиску решений теоремы Лагранжа о сумме четырёх квадратов с примерами на Си.

Небольшая статья, в которой последовательно строится поиск решения теоремы Лагранжа о том, что любое натуральное число можно представить в виде суммы квадратов четырёх целых чисел. Статья сопровождается примерами программ на языке Си. Последний пример не только весьма эффективен, но и полностью базируется на целочисленной арифметике.

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

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

Ввести натуральное N и найти для него такие целые неотрицательные $$x,y,z$$ и $$t$$, чтобы $$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 внутри циклов выполнится $$(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$$, то решения, опять-таки, нет. Таким образом, нижние оценки для всех переменных:

  1. $$(sqrt(N))/4<=x<=sqrt(N)$$

  2. $$(sqrt(N-x^2))/3<=y<=min(x,sqrt(N-x^2))$$

  3. $$(sqrt(N-x^2-y^2))/2<=z<=min(y,sqrt(N-x^2-y^2))$$

  4. $$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$$.


CategoryArticle

FrBrGeorge/News/2017-02-01 (последним исправлял пользователь FrBrGeorge 2017-02-02 16:05:54)