В этой статье рассказывается о том, как отыскать точное решение целочисленной системы линейных алгебраических уравнений на порядки быстрее, чем это делает метод Гаусса и любые его аналоги. Алгоритм работает для любых линейных систем с целыми числами, но требует, чтобы система имела единственное решение, то есть матрица системы должна быть невырожденной. Алгоритм использует для своей работы P-адические аппроксимации. Но для понимания его работы не обязательно знать, что это такое.

Итак, задана система линейных уравнений Ax=b с целыми числами. Число переменных равно n. Определитель |А|≠0. Решать такую систему можно методом Диксона, который описан в статье

Dixon J. D. Exact solution of linear equations using P-adic expansions // Numer. Math., 40, 1982, P. 137-141.

Ссылка на статью даётся в комментариях к этому посту. Кому интересны все математические выкладки, читайте оригинал. Я опишу только сам метод, без строгих доказательств. Смысл метода в том, что мы выбираем какое-нибудь простое число p и решаем систему по модулю этого числа. Далее с помощью специальной процедуры мы получаем решение по модулю pm для некоторого m>0. Число m выбирается особым образом – так, чтобы информации хватило для восстановления всех цифр ответа. После этого из данного решения восстанавливается десятичная запись ответа в виде дроби. То есть алгоритм является приближённым, но приближение выполняется не в привычных нам дробях, а в расширении поля рациональных чисел.

Моя прямолинейная реализация данного метода, запущенная на задачах конкурса, работает всего 8 с, что на порядок быстрее, чем у победителя. Разумеется, можно написать и быстрее, чем я и предлагаю заняться читателю в свободное время. Алгоритм очень прост, но потребует чуть больше усилий, чем метод исключений Гаусса.

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

Описание метода

Задача. Требуется решить целочисленную систему линейных уравнений Ax=b (|А|≠0), найдя точное решение в виде вектора дробей x.

Когда я говорю, что что-то взято по модулю p, то имеется в виду, что берётся положительный модуль. То есть −20 по модулю 7 как бы может быть равно −6 (так привыкли чистые программисты), но на самом деле имеется в виду -20 ≡ 1 (mod 7).

Шаг 0. Выбираем какое угодно простое число p. Чем больше, тем лучше, но так, чтобы p2 умещалось в машинную арифметику вашего компьютера. Считываем исходные данные. Исходную матрицу и вектор будем называть A и b. Договоримся, что массивы и матрицы с длинными числами я буду называть обычным образом (например, A), а их копии, в которых числа взяты по модулю p – с нижним индексом (например, Ap).

Шаг 1.

Вычисляем

То есть находим матрицу, обратную к матрице AP по модулю p — такую матрицу Cp, что

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

Шаг 2. Вычисляем P-адическую аппроксимацию нашего решения с помощью процедуры Lifting ( m ). Смысл числа m будет раскрыт позже.

В строке 1 инициализируется будущая последовательность векторов bk (верхний индекс — это не степень, а номер вектора). В строке 3 берём очередной элемент последовательности векторов по модулю p (используем деление длинного на короткое). В строке 4 выполняется умножение / сложение только коротких чисел. А в строке 5 требуется умножать исходные числа матрицы A на короткие числа вектора xp и находить разность между векторами из длинных чисел.

Обратите внимание, что в строке 5 деление на короткое число p всегда целое, так как


то есть в 5-й строке выражение справа в скобках будет нацело делиться на p. Если не делится – у вас где-то косяк.

После выполнения этой процедуры считаем вектор z по формуле:

Нетрудно видеть, что

Поэтому фактически, z – решение нашей исходной системы по модулю pm.

Шаг 3. К сожалению, автор не пишет, что делать в случае нечётного значения m, поэтому в своей программе я делаю m чётным. Следующий процесс делает из полученного вектора z вектор обычных дробей, который и является ответом:

Обратите внимание на строки 2 и 3. Там массивы чисел u и v начинаются от «-1», а не от «0», как мы привыкли (здесь я следую обозначению автора). Алгоритм восстанавливает вектор x (ответ) покомпонентно (цикл в строке 1). В строке 10 i-я компонента вектора x есть дробь, которую нужно сократить. Эта строка и строка 6 — единственные, где задействовано деление длинного на длинное.

Как выбрать m?

В оригинальной статье автор пишет, что «гарантируется существование m=O(n⋅log(n)), которого будет достаточно для точного восстановления ответа». Там предлагается (и более подробно об этом написано здесь, спасибо за ссылку товарищу зщл) использовать евклидову норму:

В нашем случае (в конкурсе) я делаю по-другому. Нужно оценить размер числителя и знаменателя. Для этого нужно взять максимальное число из расширенной матрицы системы (точнее, его абсолютное значение), и возвести в степень n. Затем домножить на n!. Если мы возьмем от этого логарифм по основанию p, то это и даст ту степень числа p, которой будет достаточно чтобы полностью «покрыть» точность числителя или знаменателя. Удвоенный логарифм от этого будет «покрывать» и числитель и знаменатель вместе. В конкурсе, где матрица имеет порядок n=200 и максимальное число равно 231, мы получим 2 ⋅ logp (26200⋅200!)=482 (округлил до ближайшего чётного кверху). Здесь p=231−1.

Но это, товарищи, называется «Мегахак». Я не могу доказать, что данная оценка правильна, просто у меня она работает. Надёжнее будет все-таки вам использовать рекомендации автора статьи и пользоваться евклидовой нормой. Если следовать этой рекомендации, то для задачи конкурса нужно было брать m=500, что немногим больше.

Почему m должно быть чётным? Строго говоря, не должно. Но в строке 5 программы мы делим его на 2 (половина как бы отвечает за числитель, а половина — за знаменатель), если бы оно было нечётным, нужно было бы возводить неравенство в квадрат, либо возиться с корнями от длинных чисел… Каких-то конкретных советов я в исходной работе не нашёл.

Как выбрать p?

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

Но тут есть проблема. Может оказаться так, что ваша матрица невырождена, но определитель равен нулю по модулю p (в этом случае вы даже не сможете найти обратную матрицу). Что делать? Всё очень просто – выбрать другое простое число, например, предыдущее перед 231−1. А потом снова предыдущее и т. д. Гарантируется, что слишком долго выбирать не придётся. А на случайных тестах p=231−1 подходит всегда. Сами можете себя спросить: какова вероятность, что вам попадётся матрица, с определителем, кратным выбранному вами p? Если ваше p мало (например, p=2), то эта вероятность будет очень велика. Ну, вы поняли логику.

Сложность

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

Пример

Ребята, я не дам вам исходники своей программы и не нужно просить. Кого заинтересовало – без труда напишет программу самостоятельно. А кого нет – тому исходники и не нужны. Но я приготовил для вас пример, который поможет лучше понять логику метода и написать / отладить свою программу.

Пусть n=2. Нам дана матрица, вектор правой части и ответ (чтобы сверить его в конце).

Пусть p=5, чтобы вычисления не были слишком громоздкими. Считаем

Теперь нужно выбрать m. Используя наш «Мегахак», m=2⋅ log5 (272⋅2)=9.05. То есть m=10 (если округлить вверх). Это верхняя оценка на число итераций. Но для этого примера (я уже знаю заранее) достаточно брать m=6. Понятно, что если брать больше, то ответ все равно останется правильным. Итак, m=6. Выполняем второй шаг. Заливаю его целиком одной картинкой [ следуйте процедуре Lifting ( m ) ]:

Теперь вычисляем вектор z:

Итак, pm=56=15625 и pm/2=53=125. Восстанавливаем ответ [ прямо по процедуре Restore ( z, m ) ]. Сначала найдём первую компоненту ответа:

Теперь вторую:

Как нетрудно видеть, ответ сошёлся.

Спасибо за внимание

Если кто-то хочет выразить мне благодарность за проделанную работу, то лучшим «спасибо» будет написанная вами программа, реализующаю этот метод. Вы могли бы рассказать о своём опыте на своём блоге / сайте и сделать ссылку на эту страницу. Постарайтесь перейти через предел в 3 c на тестах конкурса. То есть сам код мне не нужен, а нужно только чтобы у вас получилось его написать.

Кстати, товарищи, вы должны понимать, что описанный в этой статье метод уступает более современным алгоритмам. Так что этот метод даже и не самый лучший. Но уж больно простой для своей скорости. В этом посте в комментариях мне дали ссылки на интересные статьи. Изучите их, если хотите продвинуться в этой области. Тем более, что конкурсы на решение убойных систем уравнений ещё будут. Когда? — пока не знаю…

Обсуждение статьи на форуме в этой теме.