Алгоритм Монтгомери: (A * B) % N

Правка ru3, от burunduk3, 2018-01-14 00:08:04

Вступление

Недавно мне рассказали про некое «преобразование Монтгомери», которое якобы сводит вычисления по произвольному модулю к вычислениям по модулю 2k. Звучит интересно: во всяких вычислениях по простому модулю самое долгое место — как раз взятие по модулю после каждой операции. Но как именно преобразование работает, не рассказали. Повезло, что я живу в эпоху интернета и любую интересующую меня тему могу с тем или иным успехом загуглить.

Гугл рассказал, что штука называется «Алгоритм Монтгомери» и описана на википедии и итмошных вики-конспектах. Оба описания мне не очень понравились, потому что (хоть и объясняют, что делать) показывают гипноформулу, которая магическим образом решает нашу проблему. С вики-конспектов есть ссылка на англоязычную статью, в которой хотя бы пояснено, почему формула работает. Впрочем, пока я с этим всем разбирался, выяснил ещё пару интересных моментов, так что всем привет и поехали.

О чём речь

Задача перед нами сегодня простая: умножить два числа по модулю третьего, то есть вычислить (A × B)%N. Вот только делить с остатком мы очень не любим: в целом инструкция деления в компьютерах выполняется не очень быстро, а уж если числа длинные, то не только работает долго, так ещё и писать это нудно и сложно.

И что теперь

Фокус, который предложил Монтгомери следующий: давайте считать не (A × B)%N, а (A × B × R - 1)%N. Казалось бы, добавилось ещё и деление по модулю, тоже операция не сахар. Однако, при некоторых R это почему-то значительно упростит нашу задачу.

Вот только чем это нам поможет умножить два числа? Предлагается сначала умножить числа на R, потом умножить друг с другом, потом умножить на R - 1. Действительно, перевод X' = (X × R)%N переводит сложение в сложение, а умножение — в то самое умножение с делением на R, которое было обещано раньше: A' × B' × R - 1 = A × R × B × R × R - 1 = A × B × R = (A × B)'. Википедия такое X' называет n-остатком числа X.

Оставим пока в стороне вопрос о том, как быстро переводить числа туда-сюда и разберёмся, как же делить на R по модулю N.

Гипноформула

Та самая, которая написана везде и считает что-то полезное. (Я немного наврал в ней, но потом исправлюсь.)

 

Обратите внимание, что в правой части обычное, целочисленное деление на R, никакого деления по модулю. Интересно, в какую сторону округлять? Ещё появился какой-то непонятный N', что это вообще такое.

На самом деле, формула не так уж и странна, как кажется на первый взгляд. Многие посетители сайта codeforces.com иногда решают задачи по спортивному программированию, в некоторых таких задачах надо выдать ответ по модулю 109 + 7, а порой в промежуточных вычислениях приходится делить на два. Как же разделить число на два по модулю 109 + 7? Если число чётное, то взять и разделить, к нечётному прибавить 109 + 7 (получится чётное) и тоже разделить.

В гипноформуле такой же фокус, просто обобщённый. Мы к числу T прибавим нечто, умноженное на N (остаток по модулю не изменился) и поделим на R. Если мы подберём нечто так, чтобы деление на R произошло без остатка (вот и ответ на вопрос про округление), то получим почти правильный ответ.

Где взять нечто

На самом деле, у этой формулы есть брат-близнец. Нам нужно вести одновременное вычисление по двум разным модулям — ровно это было в КТО. Кстати, пора выяснить, что такое N'. Википедия и конспекты предлагают брать его из формулы 1 = N × N' + R × R', то есть из расширенного алгоритма Евклида для N и R. Идея хорошая, но для понимания мне было проще считать, что N' = ( - N - 1)%R (эквивалентность этих двух N' остаётся читателю в качестве упражнения). То есть, числитель гипноформулы, вычисленный по модулю R выглядит как T + T × ( - N - 1 × N) = T - T × 1 = 0. И правда, делится на R. Кстати, это даёт нам первое ограничение на R: оно должно быть взаимно просто с N.

Где я соврал

По формуле получилось, что мы вычислили правильный ответ по модулю N. То есть, может быть нужное нам число. А может быть, с добавлением 15763 × N и теперь всё равно придётся брать по модулю. С тем же успехом можно было просто вычислить A × B и сказать «это ответ, держите». Но не совсем.

Давайте оценим что получилось, считая, скажем, что 0 ≤ T < N2 (как это обычно и бывает после умножения двух чисел). Как известно, (T × N')%R < R. Остаток по модулю вообще всегда меньше модуля. То есть, числитель меньше, чем N2 + N × R, а потом ещё и на R поделили. Уже видно, что если R > N, то результат будет меньше 2 × N. То есть после применения гипноформулы, конечно, надо ещё взять по модулю, но это «взятие» — одна проверка и, возможно, одно вычитание. Так нас ждёт успех и второе ограничение на R: оно должно быть больше, чем N.

Шило на мыло

Взятие остатка по модулю N заменилось на деление на R, и ещё взятие остатка по модулю R в числители формулы. Где же польза, если мы всего лишь одно деление заменили двумя другими. Фокус в том, что мы можем подобрать такое R, на которое делить — одно удовольствие. Степень двойки, например. Или, если у нас длинная арифметика с каким-то основанием, то степень этого основания. Тогда остатки по модулю и деления превратятся в сдвиги и выкидывание лишних цифр.

Что в итоге

Получилось, что аккуратно выбрав R (например, минимальное 2k > N) мы научились без деления вычислять (A × B × R - 1)%N. Для этого нам надо заранее вычислить N' — один раз для каждого модуля. Ещё нам может пригодиться RN = R%N, а в особо запущенных случаях — даже (R2)%N.

Дальше вики-конспекты утверждают, что для умножения это не очень полезно (наглая ложь — см. ниже), зато полезно для возведения в степень. И правда, при возведении в степень нам надо один раз перейти от числа к его n-остатку, возвести в степень его, по полученному n-остатку ответа получить итоговое значение. Настоящее деления для перевода туда-сюда использовалось всего дважды, а вот умножения внутри возведения в степень все быстренькие, монтгомериевые.

Ещё лучше дела обстоят, например, с тестом Миллера—Рабина, проверкой числа на простоту. Напомню, там надо было: * Разложить P - 1 = U × 2K * Сгенерировать случайное число X. (Кому оно надо? Сгенерируем сразу n-остаток.) * Возвести его в степень U. (Получим n-остаток возведённого, переводить в обычные числа лень.) * Сравнить с единицей. (Упс. Ну ладно, сравним с n-остатком единицы, то есть с RN.) * Возводить в квадрат, попутно сравнивая с 1 и  - 1. (Ну вы поняли, продолжаем держаться в n-остатках, сравниваем с RN и P - RN).

То есть, в некоторых случаях вообще не приходится переводить числа в n-остатки и обратно. Впрочем, так ли это сложно?

Скажем, по n-остатку X' надо вычислить обычный X = (X' × R - 1)%N. Надо ли брать по модулю? Конечно, нет. И умножать ничего не надо, и даже вычислять (R - 1)%N — лишнее. Потому что делить на R по модулю N — ровно то, что делает гипноформула, надо просто применить её ещё один раз к результату.

С переводом числа в его n-остаток сложнее, но не сильно. Давайте предподсчитаем (R2)%N, а потом скажем, что $(X \times R) = (X \times R^2 \times R^{-1})$, а это мы уже умеем считать по гипноформуле. То есть, в переводе чисел к n-остаткам и обратно тоже не нужно взятие остатка по модулю, достаточно применить умножение Монтгомери ещё раз.

Вот и сказочке конец

А кто слушал — может ещё почитать, как я простые числа ищу.

$ cat primes.in 
100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000010000
$ g++ -Ofast -Wall -Wextra -Werror -o primes -x c++ primes.c++
$ time ./primes < primes.in 2> /dev/null
36

real    0m1.233s
user    0m1.230s
sys     0m0.003s
Теги арифметика, модуль, умножение

История

 
 
 
 
Правки
 
 
  Rev. Язык Кто Когда Δ Комментарий
ru18 Русский burunduk3 2018-01-15 00:11:25 49 Мелкая правка: 'з формулы $1 = R \t' -> 'з формулы <strike>test</strike> $1 = R \t'
ru17 Русский burunduk3 2018-01-14 21:21:21 10
ru16 Русский burunduk3 2018-01-14 13:02:47 2
ru15 Русский burunduk3 2018-01-14 12:49:16 7 Мелкая правка: 'проблему. С вики-конспектов есть ссыл' -> 'проблему. На вики-конспектах есть ссыл'
ru14 Русский burunduk3 2018-01-14 12:48:33 1 Мелкая правка: 'лить.[cut]\n\n\nГугл' -> 'лить.[cut] \n\n\nГугл'
ru13 Русский burunduk3 2018-01-14 12:48:11 2 Мелкая правка: '.[cut]\n\nГугл р' -> '.[cut]\n\n\nГугл р'
ru12 Русский burunduk3 2018-01-14 12:42:10 2 Мелкая правка: ' загуглить[cut].\n\nГугл р' -> ' загуглить.[cut]\n\nГугл р'
ru11 Русский burunduk3 2018-01-14 12:41:48 5 Мелкая правка: ' загуглить.\n\nГугл ' -> ' загуглить[cut].\n\nГугл '
ru10 Русский burunduk3 2018-01-14 12:41:07 0 (опубликовано)
ru9 Русский burunduk3 2018-01-14 12:40:46 12 Мелкая правка: 'авнивая с $1$ и $-1$. (Ну вы п' -> 'авнивая с **1** и **-1**. (Ну вы п'
ru8 Русский burunduk3 2018-01-14 12:39:51 703 Мелкая правка: ' \times B) \% N$. Вот то' -> ' \times B)~\%~N$. Вот то'
ru7 Русский burunduk3 2018-01-14 00:18:31 2 Мелкая правка: 'зывается «Алгоритм Мо' -> 'зывается «алгоритм Мо'
ru6 Русский burunduk3 2018-01-14 00:17:59 2
ru5 Русский burunduk3 2018-01-14 00:16:32 15 Мелкая правка: 'до было:\n- Разлож' -> 'до было:\n\n- Разлож'
ru4 Русский burunduk3 2018-01-14 00:10:14 11 Мелкая правка: 'theorem). Кстати, пора выяс' -> 'theorem). Пришла пора выяс'
ru3 Русский burunduk3 2018-01-14 00:08:04 18
ru2 Русский burunduk3 2018-01-14 00:07:08 2 Мелкая правка: ' То есть, в числитель' -> ' То есть, числитель'
ru1 Русский burunduk3 2018-01-13 23:33:27 8776 Первая редакция (сохранено в черновиках)