Китайская теорема об остатках

Сколькими способами можно закрыть часть улиц вокруг Зелёного базара так, чтобы между любыми двумя перекрёстками по-прежнему можно было проехать, но ни одного кольца не осталось? Каждый такой набор улиц — остовное дерево района, и в круге радиусом полтора километра их число из 127 цифр. В число с плавающей точкой его не уместить, а считать сразу длинными числами долго. Поэтому его считают по частям: остаток от деления на одно простое число около 67 миллионов, на другое, на третье — и каждый раз только с маленькими числами. Китайская теорема об остатках гарантирует, что по двадцати восьми таким остаткам число восстанавливается однозначно. Ниже это происходит на настоящих улицах. Перенесите круг или измените радиус.

Каждый квадратик под картой — одна «цифра» числа в смешанной системе счисления; её дал очередной модуль. Пока цифры ненулевые, число ещё не собрано: модулей мало, и их произведение меньше ответа. Когда произведение модулей перерастает ответ, цифры становятся нулями — ответ найден, дальнейшие модули его только подтверждают. Строка «в double» — тот же определитель, посчитанный обычным методом Гаусса в числах с плавающей точкой: из ста с лишним цифр верны первые тринадцать-четырнадцать.

Идея

Если модули m₁, …, m_k попарно взаимно просты, то число от 0 до m₁·…·m_k − 1 однозначно определяется своими остатками от деления на каждый из них. Поэтому вычисление, где нужен большой точный ответ, можно провести отдельно по каждому маленькому модулю — с машинными числами, без длинной арифметики, хоть на разных компьютерах — и в конце собрать ответ из остатков. Собирает его алгоритм Гарнера: ответ записывается как c₀ + c₁·m₁ + c₂·m₁m₂ + …, и каждая «цифра» c_i находится по модулю m_i+1 одним проходом по уже найденным.

Сама теорема и её классическое доказательство таблицей остатков разобраны в главе «Задача Сунь-цзы» «Царицы наук». Там же формула, которая собирает ответ из «деталей», — числа, равного единице по одному модулю и нулю по остальным. Здесь нас интересует теорема как инструмент: как восстановить число быстро, сколько модулей нужно и где это окупается.

Алгоритм Гарнера

Формула с деталями требует перемножать большие числа: деталь для модуля m_i — это произведение всех остальных модулей. Харви Гарнер в 1959 году, работая над арифметикой в системе остаточных классов для вычислительных машин, предложил собирать число иначе — по цифрам смешанной системы счисления. Запишем

По модулю m₁ все слагаемые, кроме первого, нулевые: c₀ = r₁. По модулю m₂ нули — все слагаемые, начиная с третьего, поэтому r₂ ≡ c₀ + c₁m₁, и c₁ = (r₂ − c₀)·m₁⁻¹ mod m₂. И так далее: очередная цифра вычисляется по своему модулю из уже известных, а обратные элементы даёт расширенный алгоритм Евклида. Все числа в этих вычислениях меньше самого большого модуля, длинная арифметика нужна только в самом конце, чтобы сложить c_i, умноженные на произведения модулей.

Посмотрите, как это выглядит, на календаре. Год юлианского периода — число от 1 до 7980 — задаётся тремя остатками: индиктом по модулю 15, золотым числом по модулю 19 и солнечным кругом по модулю 28. Каждый ряд сцены — одна цифра Гарнера: прыжки по 1, потом по 15, потом по 15 · 19 = 285. Прыжок на произведение прежних модулей не меняет прежних остатков — в этом весь алгоритм.

Цифры c_i обладают удобным свойством, которое видно в первой сцене: если ответ меньше m₁·…·m_j, все цифры, начиная с c_j, — нули. Поэтому вычисление можно остановить, когда несколько новых модулей подряд дали нулевые цифры. Строго это ничего не доказывает — случайно нулевая цифра выпадает с вероятностью 1/m, — и если нужна гарантия, число модулей берут по заранее известной границе ответа.

Деревья, определитель и сколько нужно модулей

Откуда в первой сцене определитель. В 1847 году Густав Кирхгоф, рассчитывая токи в электрических цепях, нашёл формулу для числа остовных деревьев графа. Возьмём матрицу, где на диагонали стоят степени вершин, а для каждого ребра uv в клетках (u, v) и (v, u) — минус единица, и вычеркнем одну строку и такой же столбец. Определитель того, что осталось, равен числу остовных деревьев. Доказательство этой формулы длиннее, чем позволяет статья; его можно найти в учебниках по теории графов под названием «матричная теорема о деревьях». Та же матрица, только с проводимостями вместо единиц, — это система уравнений из статьи о методе Гаусса, где по улицам Алматы течёт ток.

Определитель целочисленной матрицы — целое число, и по модулю простого p его считает тот же метод Гаусса, только деление заменяется умножением на обратный элемент по модулю p. Чтобы восстановить ответ, нужно, чтобы произведение модулей превысило сам ответ, а значит, нужна заранее известная граница. Для деревьев она получается совсем просто.

Граница грубая. Для района у базара она обещает до 213 цифр, а в ответе 127; в сцене это видно по хвосту нулевых цифр Гарнера. Зато она честная: сколько модулей взять, известно до начала вычислений. Для произвольных матриц такую же роль играет неравенство Адамара — определитель не больше произведения длин строк.

Ещё одна тонкость — «невезучий» модуль. Метод Гаусса в сцене не переставляет строки: матрица Кирхгофа после вычёркивания положительно определена, и в обычных числах ведущий элемент нулём не бывает. А по модулю p бывает — с вероятностью около 1/p на шаг. Такой модуль просто пропускают и берут следующий: остатков от других простых хватит.

Код

Программа берёт с карты район радиусом километр вокруг Зелёного базара, нумерует перекрёстки так, чтобы соседи получали близкие номера, — тогда ненулевые элементы матрицы лежат в узкой ленте вдоль диагонали, — и считает определитель по модулю простых чисел меньше 2³¹. Произведение двух остатков меньше 2⁶², поэтому в Go хватает uint64, а math/big нужен только в самом конце, в алгоритме Гарнера. Последняя строка — проверка модулем, который в вычислении не участвовал.

district: 237 crossings, 328 streets, matrix 236 x 236, band 38
bound: at most 97 digits
11 primes below 2^31, the largest 2147483647, the smallest 2147483423
spanning trees: 34613761769290138902561752769643080652778611353227531802400
59 digits
check modulo 1000000007: 525194220 = 525194220

Вторая программа — календарь и общий случай. Юлианский период восстанавливается двумя способами: алгоритмом Гарнера и формулой, которую Жак де Бийи опубликовал в 1665 году, а Гаусс привёл в 1801-м: год периода равен (6916·i + 4200·g + 4845·s) mod 7980. Числа 6916, 4200 и 4845 — те самые «детали»: 6916 = 19 · 28 · 13 даёт остаток 1 при делении на 15 и делится на 19 и 28, и так далее.

year       indiction golden solar   Julian period  de Billy  mixed radix
2026 AD            4     13    19            6739      6739  4 + 12*15 + 23*285
1583 AD           11      7    24            6296      6296  11 + 1*15 + 22*285
1 AD               4      2    10            4714      4714  4 + 10*15 + 16*285
1 BC               3      1     9            4713      4713  3 + 10*15 + 16*285
4713 BC            1      1     1               1         1  1 + 0*15 + 0*285

Если модули не взаимно просты, решение есть не всегда. Сравнения x ≡ r₁ (mod m₁) и x ≡ r₂ (mod m₂) совместны ровно тогда, когда r₁ ≡ r₂ по модулю g = НОД(m₁, m₂), и тогда решение единственно по модулю НОК(m₁, m₂). Функция merge находит его расширенным алгоритмом Евклида: x = r₁ + m₁t, где t решает (m₁/g)·t ≡ (r₂ − r₁)/g по модулю m₂/g. Два сравнения x ≡ 2 (mod 4) и x ≡ 3 (mod 6) несовместны — первое требует чётного x, второе нечётного, — а x ≡ 2 (mod 4) и x ≡ 4 (mod 6) дают x ≡ 10 (mod 12).

Сколько это стоит

Восстановление по k модулям алгоритмом Гарнера — O(k²) операций с машинными числами плюс O(k²) на финальную сборку длинного числа. Это мало: основная работа — вычисления по каждому модулю, их k штук, и они друг от друга не зависят — их можно раздать разным ядрам или машинам. Число модулей растёт линейно с длиной ответа: на каждые 9 десятичных цифр — один 31-битный модуль.

Окупается ли это по сравнению с точным счётом дробями? Замер: тот же определитель для районов разного размера двумя способами — по модулям и обычным методом Гаусса на дробях, которые сами сокращаются.

В Go модули выигрывают в 9–15 раз на матрицах от семидесяти строк: на матрице из 457 строк — 62 миллисекунды против 586. А в Python — почти вничью, и на больших матрицах дроби даже чуть быстрее. Дело не в алгоритме, а в языке: арифметика длинных чисел внутри Fraction написана на C, а цикл по модулям — на Python, и его накладные расходы съедают выигрыш от маленьких чисел. Та же история, что со встроенным умножением в статье об умножении Карацубы: считать операции полезно, но решает замер.

Где это работает

• Точная линейная алгебра. Библиотека FLINT считает определитель целочисленной матрицы от 25 строк по многим модулям, а матрицы поменьше — методом Барейсса, без дробей. Тот же приём — у многочленов: наибольший общий делитель и разложение на множители удобнее искать по модулю простых чисел и потом поднимать результат.
• RSA. Закрытая операция c^d mod N считается отдельно по модулям p и q — с числами и показателями вдвое короче — и склеивается по китайской теореме. Возведение в степень стоит примерно куб длины, каждая половина — восьмую часть целого, и обе вместе — около четверти; поэтому в закрытом ключе RSA хранят не только d, но и p, q, d mod (p − 1), d mod (q − 1) и q⁻¹ mod p. У ускорения есть цена: если сбой оборудования испортит одну из половин, неверная подпись выдаёт множитель N через НОД, — поэтому реализации проверяют подпись перед тем, как её отдать.
• Умножение длинных чисел. Теоретико-числовое преобразование перемножает многочлены точно, но по модулю простого числа. Если коэффициенты произведения больше модуля, берут три простых вида c·2^k + 1, например 998 244 353 = 119·2²³ + 1, и склеивают три результата алгоритмом Гарнера.
• Системы остаточных классов. В представлении остатками сложение и умножение идут по каждому модулю независимо, без переносов между разрядами. Ради этого Гарнер и писал свою работу: такие системы изучали для вычислительных машин и применяют в специализированной аппаратуре цифровой обработки сигналов. Сравнивать числа и делить в остатках зато трудно — для этого и нужен перевод в обычную запись.
• Хронология. Юлианский период из второй сцены придумал Иосиф Скалигер в 1583 году, чтобы перевести даты разных календарей на один счёт; его первый год, 4713 год до н. э., — тот, где все три цикла равны единице. С 1849 года, после книги Джона Гершеля, астрономы считают от полудня 1 января этого года юлианские дни — и считают до сих пор.

Где спотыкаются

• Модули с общим делителем. Формула с деталями и алгоритм Гарнера требуют взаимно простых модулей: обратного к 4 по модулю 6 не существует. Для общего случая нужна проверка совместности и склейка через НОК, как в функции merge.
• Мало модулей. Если произведение модулей меньше ответа, Гарнер честно вернёт остаток ответа по этому произведению — правдоподобное и неверное число. Граница ответа нужна до вычислений.
• Отрицательный ответ. Определитель бывает отрицательным, а остатки — нет. Восстановленное x из промежутка [0, M) нужно понимать как число из [−M/2, M/2): если x > M/2, ответ равен x − M. Модулей тогда нужно на один больше.
• Переполнение. В Go и C произведение двух остатков по модулю около 2⁶³ не помещается в 64 бита. Модули до 2³¹ — самый простой способ обойтись без 128-битного умножения.

История

Задача о предметах, которые при счёте тройками, пятёрками и семёрками дают остатки 2, 3 и 2, есть в трактате «Сунь-цзы суань цзин», написанном в III–V веках; общий способ решать такие системы изложил Цинь Цзюшао в 1247 году. На языке сравнений теорему сформулировал и доказал Гаусс в «Арифметических исследованиях» 1801 года — там же, среди примеров, юлианский период. Как алгоритм для вычислительных машин она вернулась в 1959 году в работе Харви Гарнера о системе остаточных классов. Позже многомодульные вычисления стали стандартным приёмом компьютерной алгебры: так ищут наибольший общий делитель многочленов, где точные коэффициенты разрастаются, и так считают определители.

Задачи

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