Равные суммы одинаковых степеней
(эксперимент)

4 февраля 2022
(обновлено 9 февраля 2022)
Время на чтение: 36 минут

Эта статья дополняет и завершает серию постов в блоге о челлендже «Equal Sums of Like Powers». В ней я расскажу о том, как начиная с наивного решения, шаг за шагом была написана самая первая версия программы для нахождения решений диофантовых уравнений вида k.1.n. И начну я с исторической справки.

Но хочу заранее предупредить, что материал получился очень объёмным (TL;DR). Поэтому, если желаете, то можете сразу перейти к самой сути или финальным результатам.

Примерно в III веке нашей эры древнегреческий математик Диофант Александрийский написал книгу, которая называлась «Арифметика». Она была посвящена нахождению решений уравнений определённого вида, которые теперь называют диофантовыми. В общем виде эти уравнения можно записать так:

P(x1, x2, ..., xi) = 0
(1.1)

Здесь P — это некоторая целочисленная функция, например полином с целыми коэффициентами, а xi — неизвестные переменные. Особенность диофантовых уравнений в том, что их коэффициенты и переменные могут принимать только целые значения.

Математиков с очень давних времён интересовали эти изящные и простые, на первый взгляд, уравнения. Их изучением занимались такие известные учёные как Леонардо Фибоначчи, Франсуа Виет, Симон Стевин, Пьер де Ферма (считается одним из создателей теории чисел), Леонард Эйлер, а также очень многие современные математики.

Некоторые виды диофантовых уравнений очень хорошо изучены. Но тем не менее, общего алгоритма для их решения не существует. В 1900 году Давидом Гильбертом нахождение универсального метода определения разрешимости произвольного алгебраического диофантова уравнения было включено в список 23 важнейших проблем математики. И только в 1970 году Юрий Матиясевич завершил доказательство алгоритмической неразрешимости этой задачи, которое заняло у него около 20 лет!

• • •

Существует огромное количество вариантов диофантовых уравнений, отличающихся значением коэффициентов, их знаком, количеством переменных и разнообразием степеней. Однако в данной статье (и эксперименте) нас будет интересовать вполне конкретный вид. Потому что далее рассказ пойдёт об одной цепи конкретых событий. Для начала рассмотрим уравнение:

a1k + a2k + ... +amk = b1k + b2k + ... +bnk
(1.2)

где k, m, n, ai и bj — целые положительные числа. И для уравнения выполняются следующие условия:

a1 ≥ a2 ≥ ... ≥ am
(1.3)
b1 ≥ b2 ≥ ... ≥ bn
(1.4)
ai ≠ bj
(1.5)
a1 > 1
(1.6)
m ≤ n
(1.7)

Договоримся обозначать уравнения такого вида при помощи записи k.m.n, где k — степень уравнения, m — количество переменных в левой части, а n — количество переменных в его правой части.

В 1769 году Эйлером была выдвинута гипотеза, являвшаяся обобщением Великой теоремы Ферма. Эта гипотеза утверждала, что для любого натурального числа n > 2 никакую n-ю степень натурального числа нельзя представить в виде суммы (n − 1) n-х степеней других натуральных чисел. Другими словами, согласно этой гипотезе, следующие диофантовы уравнения не имеют решений в натуральных числах:

z3 = a3 + b3
(1.8)
z4 = a4 + b4 + c4
(1.9)
z5 = a5 + b5 + c5 + d5
(1.10)
. . .  

В 1966 году трое американских учёных Л. Ландер, Т. Паркин и Дж. Селфридж с помощью суперкомпьютера CDC 6600 нашли первый контрпример для n = 5, который опроверг гипотезу Эйлера.

В 1988 году Ноам Элкис с помощью эллиптических кривых нашёл контрпример для случая n = 4. А немного позднее в том же году Роджер Фрай его улучшил, найдя наименьшее возможное решение с помощью параллельных вычислений.

Любопытно, что для n = 6 и бо́льших значений такой пример до сих пор так и не найден.

Всё это безумно интересно! Но так как повествование уже затягивается, далее я лишь дам ссылку на страницу о Великой теореме Ферма (из которой вытекает справедливость для 1.8) и упомяну о гипотезе Ландера — Паркина — Селфриджа, а сам пойду дальше (желающие смогут всё узнать самостоятельно со страниц Википедии).

• • •

Итак. О каком же челлендже я говорил в самом начале? Я собираюсь повторить поиск минимального решения для уравнений 4.1.3 (1.9) и 5.1.4 (1.10), используя современное оборудование. Мне любопытно узнать, сколько времени займут эти вычисления на домашнем компьютере, производительность которого по современным мерками чуть выше средней.

Для этого я напишу универсальную программу (подходящую для любых уравнений вида k.1.n), перебирающую все решения, оптимизирую её алгоритм и выясню, за какое время будет найдено первое минимальное решение для двух интересующих меня уравнений на современном ПК с процессором Intel Core i5-10400, если производить вычисления на всех его 6 ядрах в 12 потоков (рабочая частота этого ЦПУ при такой нагрузке равна 4.0 GHz).

Начав с наивного решения, я буду шаг за шагом его улучшать, объясняя свои действия и приводя результаты этих улучшений. Измерение скорости работы программы будет выполняться на поиске решений для уравнения 4.1.4, так как его решения меньше и находятся значительно быстрее. Весь исходный код этого эксперимента расположен в основном репозитории проекта iLWN (проект Lab). Но здесь я буду использовать только Gists и только для основных частей кода программы.

Также для своего удобства я буду использовать библиотеку AML (класс файла, функции форматирования и работу с консолью). Но не буду давать дополнительных пояснений, так как думаю, что всё будет интуитивно понятно. Исходный код библиотеки снабжён комментариями и находится в публичном репозитории, поэтому при желании можно будет разобраться самостоятельно.

• • •

Итак, начнём с наивного решения. Для простоты я буду использовать 64-битные переменные, что позволит проверить все возможные комбинации до значения старшего коэффициента 49,796 включительно (для уравнения 4.1.3). Программа будет последовательно перебирать все возможные варианты коэффициентов правой части, вычислять сумму их степеней и затем проверять, существует ли такой коэффициент для левой части, степень которого равна этой сумме.

На мой взгляд, это очевидное и достаточно простое решение. Но чтобы делать меньше арифметических операций внутри цикла, сразу предлагаю простую оптимизацию: мы заранее вычислим нужные степени всех чисел от 1 до 49,796 и поместим их в массив. Сделаем это в отдельной функции CalcPowers():

В основной функции мы сначала инициализируем массив со степенями, затем массив коэффициентов правой части уравнения. Все возможные значения левой части помещаем в std::map для простоты проверки. Внутри основного цикла: считаем сумму степеней коэффициентов правой части уравнения, ищем полученное значение в std::map. Если оно найдено, значит мы нашли решение. Далее выбираем следующий набор коэффициентов и повторяем цикл.

Функция SearchForFactors() принимает 3 параметра: степень уравнения (power), количество коэффициентов правой части (count) и значение старшего коэффициента (hiFactor), с которого начинается проверка. Вот код этой функции:

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

Итак, у нас есть первая версия алгоритма! Запускаем и ждём результат: программе понадобилось примерно 150 секунд, чтобы найти первые 3 решения:

Searching for factors (4.1.4)
Factor limit is set to 46340
Solution found: 353=315+272+120+30
Solution found: 651=599+430+340+240
Solution found: 706=630+544+240+60

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

• • •

Чтобы убрать из вывода программы непримитивные решения, нам понадобится функция GetPrimes(), вычисляющая все простые числа до заданного. Используем простую реализацию решета Эратосфена (каждый бит массива bits соответствует одному нечётному числу, начиная с 3):

Создадим структуру Solution, в которой будем хранить решение:

Также создадим класс контейнера Solutions, в который будем добавлять все найденные решения:

Проверку решения выполняет функция Solutions::IsPrimitive(). В ней мы проверяем делимость каждого коэффициента решения на простые числа, не превышающие значения самого младшего коэффициента. Если мы сможем найти такое простое N, которое будет делить без остатка каждый коэффициент решения, то это будет означать, что решение не примитивное. Такое решение мы не будем добавлять в контейнер.

Теперь проверка решения в главной функции выглядит так:

• • •

Вывод непримитивных решений мы отфильтровали. Однако, есть ещё одна, не совсем очевидная проблема. Лучше будет показать её на примере вывода для уравнения 3.1.3. Собираем новую версию программы и запускаем:

Searching for factors (3.1.3)
Factor limit is set to 49999
Solution found: 6=5+4+3
Solution found: 9=8+6+1
Solution found: 20=17+14+7
Solution found: 19=18+10+3
Solution found: 28=21+19+18
Solution found: 25=22+17+4
Solution found: 29=27+15+11
Solution found: 41=33+32+6
Solution found: 46=37+30+27
Solution found: 46=37+36+3
Solution found: 41=40+17+2

Непримитивные решения пропали из вывода. Но, как мы видим, сейчас программа считает «старшим» наибольший коэффициент в правой части. Тогда как на самом деле старший коэффициент уравнения — тот, что находится в его левой части. Из-за этого не только нарушается порядок вывода решений, но и могут быть пропущены некоторые из них, когда программа достигнет ограничения по максимальному значению старшего коэффициента.

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

Пусть теперь k[0] будет хранить коэффициент левой части уравнения. Тогда условие главного цикла изменится на while (k[0] <= maxFactor), а проверка решения заменится условием if (sum == powers[k[0]]). Вот полный код изменённой функции:

Из-за внесённых изменений время работы программы для уравнения 4.1.4 увеличилось как минимум в несколько раз (см. подробности далее). Но теперь её вывод (для 3.1.3) выглядит так, как должен:

Searching for factors (3.1.3)
Factor limit is set to 49999
Solution found: 6=5+4+3
Solution found: 9=8+6+1
Solution found: 19=18+10+3
Solution found: 20=17+14+7
Solution found: 25=22+17+4
Solution found: 28=21+19+18
Solution found: 29=27+15+11
Solution found: 41=33+32+6
Solution found: 41=40+17+2
Solution found: 44=41+23+16
Solution found: 46=37+30+27
Solution found: 46=37+36+3
• • •

Программа работает корректно, хоть и очень медленно. Теперь можно приступать к оптимизации алгоритма. Но прежде, чем я начну, я хочу добавить вывод результатов в файл, заменить вызовы printf() на aux::Printf(), добавить периодический вывод прогресса и возможность прерывания работы программы. И конечно, нужно добавить вывод времени работы для того, чтобы можно было оценить влияние оптимизаций.

Я не буду вставлять здесь получившийся код, так как он теперь стал значительно больше. Посмотреть содержимое файла можно тут. Воспользовавшись кнопкой History (История) можно сравнить новую версию файла с предыдущей.

• • •

Теперь моя задача — ускорить работу программы. И первое, о чём стоит подумать, это о падении скорости после отказа от использования std::map. Первая версия программы примерно за 150 секунд работы успела выдать 3 решения и дойти до старшего коэффициента 630. Новая же версия, проработав 10 минут, выдала лишь одно решение, добравшись только до значения 448. Так почему же это произошло?

В первой версии мы перебирали только коэффициенты правой части, а в новой — все, то есть на 1 больше. Для уравнения 4.1.4, в котором всего 5 коэффициентов, при значении старшего коэффициента, равном 20, количество итераций основного цикла выросло в 5.75 раза, при коэффициенте 30 — в 8.25 раза, при 50 — в 13.25 раза, при 100 — в 25.75 раза! Легко заметить, что влияние дополнительного коэффициента нелинейно и тем сильнее, чем дальше поиск. Поэтому несмотря на то, что цикл стал проще — в нём теперь нет вызова std::map::find(), — такое значительное увеличение количества итераций имеет крайне негативный эффект.

Хорошая новость в том, что мы можем перебирать на один коэффициент меньше и при этом «искать» последний тем же методом, что и раньше. И для этого нам даже не будет нужен std::map. Если мы вычтем из значения левой части уравнения степени каждого коэффициента правой части, кроме последнего, то полученное значение будет степенью последнего коэффициента. Всё, что нам останется сделать — найти такое целое число, степень которого была бы равна полученной разности. И если такое число существует, то вместе с остальными коэффициентами оно даст нам решение.

Так как теперь в сумме степеней будет на одно слагаемое меньше, то в качестве бонуса мы также получим более высокое значение переменной hiFactor. Максимальное значение старшего коэффициента (для уравнения 4.1.3) увеличится с 49,796 до 55,108!

Так как массив степеней упорядочен по возрастанию, мы можем использовать бинарный поиск, чтобы найти число по его степени. По сложности это и будет примерно соответствовать std::map::find(). И так как коэффициенты решения должны располагаться в невозрастающем порядке, то искомое значение не будет превышать значения предпоследнего коэффициента, что сделает поиск ещё быстрее. И теперь начало основного цикла нашего алгоритма будет выглядеть следующим образом:

Новый оптимизированный алгоритм работает даже лучше, чем я ожидал. До значения старшего коэффициента, равного 700, на этот раз программа добралась всего за 101 секунду! Это примерно в 1.49 раза быстрее, чем исходный вариант с std::map (который закончил работу даже раньше, ещё на значении 630)! И более, чем в 6 раз быстрее варианта, предшествующего оптимизации!

• • •

Перед тем, как двигаться дальше, я сделал небольшой рефакторинг: вынес вывод прогресса в отдельную функцию UpdateProgress(), чтобы сократить тело главного цикла и сделать код понятнее. Изменение незначительное, поэтому останавливаться на нём буду.

Далее нам следует разделить основной цикл, в котором сейчас перебираются все коэффициенты, на два цикла: внешний, в котором мы будем инкрементировать старший коэффициент, и внутренний, в котором мы будем перебирать коэффициенты правой части уравнения.

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

После перечисленных изменений время работы программы (для уравнения 4.1.4, до k[0] == 700) сократилось до 91 секунды, то есть ещё в 1.11 раза. Совсем немного, но тоже неплохо. Вот так теперь выглядит начало цикла:

А так выглядит его конец:

• • •

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

Например, когда в левой части значение коэффициента равно 100, нет смысла начинать перебор с единицы: для первого коэффициента правой части все значения, меньшие 71, не имеют смысла, так как если даже все остальные коэффициенты будут равны 70, то сумма их степеней всё равно окажется меньше значения левой части!

Также возникает ситуация, когда сумма степеней коэффициентов правой части уже превышает значение левой. В таких случаях перебор нужно прерывать и переходить к следующему значению коэффициента. Именно это мы сейчас и сделаем, изменив всего 2 строчки в конце внутреннего цикла:

Эта очень простая оптимизация сократила время работы программы до 62.2 секунды, или ещё в 1.46 раза!

• • •

Далее я снова сделал небольшой рефакторинг. Вынес открытие файла в отдельную функцию OpenLogFile(). А проверку и вывод решения — в функцию OnSolutionFound(). Вот листинг этой функции:

Теперь главный цикл функции SearchForFactors() стал ещё компактнее:

Вот здесь находится исходный файл со всеми сделанными изменениями.

И перед тем, как продолжить оптимизацию алгоритма, мне бы хотелось кое-что ещё изменить в выводе программы. Речь идёт о том, как выглядят одинаковые коэффициенты правой части уравнения. Вот, например, что сейчас выдаёт программа для уравнения вида 4.1.6:

Searching for factors (4.1.6)
Factor limit is set to 43826
Solution found: 3=2+2+2+2+2+1
Solution found: 7=6+4+4+4+4+3
Solution found: 9=7+6+6+6+4+2
Solution found: 9=8+7+2+2+2+2
Solution found: 13=12+8+7+6+2+2
Solution found: 17=14+14+8+6+6+1

Идея заключается в объединении одинаковых коэффициентов и записи их в виде ...+n*k+..., где n — количество повторений, а k — значение коэффициента. Так, первое решение из вывода, показанного выше, можно будет записать как 3=5*2+1, а второе — как 7=6+4*4+3. Это позволит сделать вывод решений более компактным для уравнений с большим количеством коэффициентов в правой части. Внесём необходимые изменения в функцию OnSolutionFound():

И теперь тот же самый вывод выглядит вот так:

Searching for factors (4.1.6)
Factor limit is set to 43826
Solution found: 3=5*2+1
Solution found: 7=6+4*4+3
Solution found: 9=7+3*6+4+2
Solution found: 9=8+7+4*2
Solution found: 13=12+8+7+6+2*2
Solution found: 17=2*14+8+2*6+1
• • •

Снова возвращаемся к оптимизации. Мы уже пропускаем значения коэффициентов, которые больше, чем необходимо. Пора сделать то же самое для значений, которые меньше. Как я говорил выше, нет смысла всегда начинать перебор с единицы. Когда мы обнаруживаем, что коэффициент достиг верхнего предела своего значения, мы должны увеличить на 1 коэффициент слева от него, и начать перебор «сначала».

Каждый раз это начальное значение будет разным: оно зависит от суммы степеней предыдущих коэффициентов. И так как эту сумму мы знаем, то несложно вычислить, с какого значения нужно начать перебор, чтобы не проверять те решения, при которых сумма в правой части всегда будет меньше, чем необходимо. Это должно быть минимально возможное значение, при котором сумма степеней коэффициентов, стоящих слева, и степеней максимальных значений коэффициентов правее, включая искомый, будет превышать значение левой части уравнения.

Чтобы это реализовать, нам понадобится считать сумму степеней всех коэффициентов правой части. Но одно сделанное недавно изменение будет мешать. Помните, что мы сделали, когда перестали учитывать в сумме правой части значение степени последнего коэффициента? Мы стали при инициализации массива степеней передеавать в функцию CalcPowers() значение count — 1 вместо count. Вот это изменение придётся откатить, снова понизив максимальное значение старшего коэффициента (то есть переменной maxFactor).

Начало внешнего цикла немного изменилось. Теперь мы быстро пропускаем малые значения первого коэффициента правой части, начиная во внутреннем цикле сразу с нужного значения:

В конце внутреннего цикла, где мы выбирали следующий набор, измения будут более существенные. Добавим новую переменную idx и проинициализируем её нолём. В случае, когда коэффициент будет сброшен в единицу, мы запомним индекс самого левого сброшенного коэффициента. Выйдя из цикла инкремента коэффициентов, мы проверим значение idx и если оно будет отлично от 0, то мы вычислим «нужные» значения всех сброшенных коэффициентов.

Чтобы найти минимальное значение коэффициента, большее 1, такое, что будет возможно решение, мы можем просто увеличивать значение в цикле на 1 и каждый раз проверять сумму. Но можно поступить оптимальнее: прибавлять в каждой итерации не 1, а половину значения предыдущего коэффициента, уменьшая это приращение вдвое каждую итерацию.

Приведу пример. Пусть у нас уравнение 4-й степени и 4 коэффициента в его правой части. И только что был обнулён второй. Нам нужно подобрать начальное значение для этого коэффициента. Положим, что значение коэффициента в левой части равно 300, а значение первого коэффициента в правой — 240. Тогда минимальное значение, удовлетворяющее условию должно быть равно 200.

И чтобы его найти, мы последовательно проверим значения 1+120 (мало, принимаем), 121+60 (мало), 181+30 (много, пропускаем), 181+15 (мало), 196+7 (много), 196+3 (мало), 199+1 (мало). Таким образом, всего за 7 итераций мы увеличили начальное значение коэффициента с единицы до 200, с которого и начнём перебор второго коэффициента. Вот код с изменениями:

С этой оптимизацией скорость работы программы значительно выросла, так как мы отбросили огромное количество наборов коэффициентов, которые заведомо не могли дать решение. Время работы для уравнения вида 4.1.4 до достижения k[0] == 700 сократилось до 4.65 секунды, то есть примерно в 13.4 раза относительно предыдущего раунда оптимизации.

• • •

Есть ещё один интересный приём. Но сначала вопрос. Что именно, какую операцию больше всего раз (и дольше) мы делаем в программе? На что уходит больше всего времени? Правильный ответ — поиск (проверка) значения последнего коэффициента в масиве степеней! Именно эта операция выполняется в каждом цикле перебора коэффициентов. Её оптимизация должна дать хорошее ускорение. Поэтому при переходе от последовательного перебора к бинарному поиску в начале статьи мы получили очень сильный эффект.

Но как можно уменьшить количество итераций цикла для бинарного поиска? ... Никак! Но можно сделать так, чтобы частенько обходиться вообще без него! Секрет здесь в том, что намного чаще искомого значения в массиве нет и поиск напрасен. И если бы можно было сразу понять, что такого значения в массиве нет, мы бы сэкономили кучу времени. И такой способ есть!

Мы можем заранее вычислить хеши для всех значений в массиве степеней и поместить их в очень компактный массив (размер которого будет сопоставим с размером кеша L1/L2 процессора). Проверка значения в такой хеш-таблице — это очень быстро. Но так как значений в хеш-таблице меньше, чем возможных значений коэффициентов, то на каждый элемент хеш-таблицы будет приходиться несколько из значений массива степеней. Но это совсем не страшно.

Таким образом, если мы обнаружим в хеш-таблице хеш нужного значения степени, то это будет означать, что искомое значение коэффициента может существовать и нам нужно выполнить поиск в массиве как обычно. А вот если в хеш-таблице значения не окажется, это будет гарантировать нам то, что искомого значения не существует! Ниже приведён код класса HashTable, в котором за проверку значения отвечает функция Exists():

В начале функции SearchForFactors(), сразу после инициализации массива степеней, инициализируем нашу хеш-таблицу:

Теперь внесём изменения в начало внутреннего цикла. Непосредственно перед поиском исходного числа для lastFP (степень последнего коэффициента правой части) добавим проверку хеша:

Благодаря сделанным изменениям удалось сократить время работы программы ещё сильнее. Для того же уравнения вида 4.1.4 до достижения старшим коэффициентом значения k[0] == 700 теперь требуется всего 0.55 секунды! Таким образом, работа программы стала быстрее в 8.5 раза относительно предыдущего раунда оптимизации. А относительно самого первого (наивного) решения, которому требовалось более 150 секунд, разница составила более 270 раз!

На этом, пожалуй, я закончу. Вот ссылка на финальный исходный файл программы. Я продолжу дорабатывать её код: добавлю поддержку 128-битной арифметики и многопоточность. Но это произоизойдёт уже за рамками этой статьи. Весь исходный код можно будет найти в репозитории iLWN (проект ESLP).

• • •

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

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

Тем не менее, полученная программа проста и универсальна: её можно использовать для поиска решений любых диофантовых уравнений вида k.1.n. А при необходимости её можно сравнительно легко доработать для решения любых уравнений вида k.m.n.

• • •

В качестве эпилога я покажу, как можно немного сократить количество проверяемых комбинаций для уравнений вида 4.1.3 и 4.1.4. Речь пойдёт о модульной арифметике. И сначала рассмотрим, может ли быть коэффициент в левой части уравнения чётным. Вполне очевидно, что:

Z ≡ 0 (mod 2)  ⇒  Z4 ≡ 0 (mod 16)
(2.1)

Также несложно убедиться, что:

Z ≢ 0 (mod 2)  ⇒  Z4 ≡ 1 (mod 16)
(2.2)

Рассмотрим диофантово уравнение вида 4.1.3 (для 4.1.4 все рассуждения будут аналогичны):

z4 = a4 + b4 + c4
(2.3)

Если z чётно, а также чётны a, b, и c, то очевидно, что такое решение не примитивно, так как каждый коэффициент делится на 2 без остатка. Так как нас не итересуют подобные решения, то предположим, что z — чётное, и хотя бы один из коэффициентов правой части нечётный. В таком случае, по свойству сравнений, учитывая (2.2), мы получим, что:

a4 + b4 + c4 ≡ 1, 2, или 3 (mod 16)
(2.4)

Так как z чётно, то выполняется (2.1), следовательно, при таких условиях левая часть уравнения (2.3) не будет равна правой, а значит такой набор коэффициентов решением не будет. Таким образом, можно сделать вывод, что при чётных значениях z непримитивных решений не существует. Это означает, что перебирая все возможные наборы коэффициентов, мы можем пропустить те из них, в которых старший коэффициент чётный. И проверять только те наборы, в которых он нечётный. Это должно сократить количество всех проверяемых комбинаций вдвое.

Рассмотрим теперь делимость старшего коэффициента на 5. Аналогично (2.1) и (2.2), мы можем легко убедиться, что:

Z ≡ 0 (mod 5)  ⇒  Z4 ≡ 0 (mod 5)
(2.5)
Z ≢ 0 (mod 5)  ⇒  Z4 ≡ 1 (mod 5)
(2.6)

Выполняя рассуждения аналогичные тем, что мы делали для делимости на 2, мы можем прийти к следующему выводу. Если (z, a, b, c) является непримитивным решением уравнения (2.3), то:

a4 + b4 + c4 ≢ 0 (mod 5)
(2.7)

Утверждение (2.7) справедливо для уравнений вида 4.1.3 и 4.1.4 (но не для 4.1.5 и других). И означает, что если коэффициент в левой части (z) делится на 5 без остатка, то не существует таких целых чисел a, b и c, которые могли бы дать непримитивное решение. Поэтому при переборе мы можем пропускать не только чётные значения старшего коэффициента, но также и те, что кратны 5.

• • •

Update Я решил дополнить статью результатами эксперимента, чтобы Вам не пришлось искать в блоге соответствующие записи (в них не будет ничего, о чём бы я не рассказал в этой статье). Но напомню, что решение, описанное в данной статье, для эксперимента было существенно доработано.

Была реализована 128-битная арифметика (на интринсиках c++), а также многопоточность. В эксперименте участвовали две версии программы: первая была без дополнительных оптимизаций (кроме исключения старших коэффициентов, которые чётны или кратны пяти), и вторая — специализированная, в которой были реализованы все имеющиеся возможности для уменьшения количества перебираемых вариантов для уравнения вида 4.1.3.

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

Решение для уравнения 4.1.3 было найдено универсальной программой за 204 часа и 48 минут (то есть практически за 8.5 суток). Это было очень долго, но я решил не останавливать программу и всё-таки дождаться, когда она завершит нужный объём вычислений.

Вторая оптимизированная версия со специализированным для 4.1.3 алгоритмом нашла первое минимальное решение всего за 2 минуты и 28 секунд! А в течение суток нашла и второе! Таким образом, разница в скорости их работы составила около 5 тысяч раз! Но интересно даже не это.

В 1988 году Роджер Фрай использовал для вычислений суперкомпьютер Connection Machine. Его архитектура (CM-1/CM-2) совершенно не похожа на то, как устроены современные RISC-процессоры. Однако в то время эта машина обладала довольно большой вычислительной мощностью. Благодаря техническому прогрессу современные процессоры общего назначения в десятки и сотни раз превосходят по производительности суперкомпьютеры конца 80-х годов 20-го века.

Ниже перечислены найденные программой минимальные решения для уравнения 4.1.3 (3.1 и 3.2), а также уравнения 5.1.4 (3.3):

4224814 = 4145604 + 2175194 + 958004
(3.1)
28130014 = 27676244 + 13904004 + 6738654
(3.2)
1445 = 1335 + 1105 + 845 + 275
(3.3)