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

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

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

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

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

Ввести натуральное 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=0; 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 count = 0;
   6 
   7     scanf("%d", &N);
   8     for(x=0; 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                     count++;
  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",count);
  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 count = 0;
   6 
   7     scanf("%d", &N);
   8     for(x=0; 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                     count++;
  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",count);
  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=0; 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 count = 0;
   6 
   7     scanf("%d", &N);
   8     for(x=0; 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                     count++;
  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",count);
  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 count = 0;
   6 
   7     scanf("%d", &N);
   8     for(x=0; 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                     count++;
  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",count);
  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 count = 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                     count++;
  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",count);
  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 count = 0;
   6 
   7     scanf("%d", &N);
   8     for(x=0; (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                     count++;
  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",count);
  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 Y, Z;
   6     long int count = 0;
   7 
   8     scanf("%d", &N);
   9     for(x=0; x*x*4<N; x++);
  10     for(Y=x; x*x<=N; x++) {
  11         while(Y && (Y-1)*(Y-1)*3>=N-x*x) Y--;
  12         for(y=Z=Y; y<=x && x*x+y*y<=N; y++) {
  13             while(Z && (Z-1)*(Z-1)*2>=N-x*x-y*y) Z--;
  14             for(z=t=Z; z<=y && x*x+y*y+z*z<=N; z++) {
  15                 while(t*t>N-x*x-y*y-z*z) t--;
  16                 count++;
  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",count);
  23 }

Пришлось дополнительно отследить нулевые значения переменных на случай, когда левая часть неравенства вида $$2(Z-1)^2>=N-x^2-y^2$$ увеличится после уменьшения $$Z$$.


CategoryArticle

FrBrGeorge/News/2017-02-01 (last edited 2017-02-02 13:05:54 by FrBrGeorge)