Применение элементарной теории чисел в программировании

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

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

От читателя этого текста предполагается знание алгебры в объеме первого семестра мехмата. Ключевые понятия и теоремы: кольцо вычетов по модулю m, Малая теорема Ферма, группа, теорема Лагранжа, прямая сумма (прямое произведение) колец и групп. Очень важным утверждением является Китайская теорема об остатках и соответствующий алгоритм, которые, к сожалению, не рассматриваются в обязательном курсе алгебры. К важнейшим алгоритмам относятся расширенный алгоритм Евклида (для заданных целых чисел m и n найти их наибольший общий делитель d и его выражение в виде линейной комбинации исходных чисел d = u m + v n) и алгоритм быстрого возведения в степень. Они рассматривались в предыдущем семестре и здесь повторяться не будут. Из расширенного алгоритма Евклида непосредственно получается алгоритм нахождения обратного элемента в кольце вычетов по модулю m. Отметим, что при программировании на языке Java мы будем использовать class BigInteger (принадлежащий пакету java.math), который входит в состав JDK 1.1. Все упомянутые алгоритмы реализованы в языке Java и являются методами класса BigInteger.

Основным объектом, с которым мы будем работать, является кольцо вычетов по модулю m (фактор кольцо кольца целых чисел по главному идеалу, порожденному элементом m). Будем обозначать его Z_m. Заметим, что в программировании практически всегда работают не с целыми числами, а с элементами колец вычетов. То, что называется целым числом в программировании (тип "int"в СИ), представляет собой на самом деле элемент кольца вычетов по модулю 2³² (или 2⁶⁴ на современных процессорах). "Настоящие" целые числа при перемножении и возведении в степень очень быстро растут, и работать с ними невозможно.

Элементы кольца вычетов Z_m, т.е. классы эквивалентности, представляются обычно числами 0, 1, ..., m-1. Напомним, что класс эквивалентности состоит из чисел, попарная разность которых делится на m. Система представителей классов (по одному из каждого класса) называется "системой остатков". Наряду с вышеприведенной системой остатков, состоящей из неотрицательных целых чисел, применяется также симметричная система, состоящая из отрицательных и неотрицательных чисел, не превосходящих по модулю m/2. Правда, симметричной она будет только в случае нечетного m. При четном m обычно выбирают -m/2 в качестве представителя класса {..., -m/2, m/2, 3m/2, ...}, так что количество отрицательных представителей на единицу больше, чем положительных. К примеру, байты, т.е. элементы кольца вычетов по модулю 256, можно интерпретировать либо как неотрицательные числа в диапазоне 0...255, либо как числа со знаком в диапазоне -128...127. Этому соответсвуют целочисленные типы unsigned char и char в СИ.

Напомним необходимые нам результаты из элементарной теории чисел и алгебры.

Малая теорема Ферма. Пусть p -- простое число. Тогда для всякого целого числа b, отличного от нуля, справедливо сравнение

b^(p-1) == 1 (mod p).

Малая теорема Ферма является непосредственным следствием теоремы Лагранжа (порядок любого элемента группы делит порядок группы) и того факта, что кольцо Z_p в случае простого p является полем, т.е. все его ненулевые элементы принадлежат группе обратимых элементов. Порядок группы обратимых элементов кольца Z_p равен p-1.

Китайская теорема об остатках. Пусть m₁, m₂, ..., m_k -- попарно взаимно простые целые числа, r₁, r₂, ..., r_k -- произвольные целые числа. Тогда найдется целое цисло x такое, что

x == r₁ (mod m₁),
x == r₂ (mod m₂),
   . . .
x == r_k (mod m_k).

В более математической формулировке Китайская теорема об остатках выглядит следующим образом: пусть

m = m₁ m₂ ... m_k.

Тогда кольцо вычетов по модулю m изоморфно прямой сумме (или прямому произведению) колец вычетов Z_mi:

Z_m == Z_m1 + Z_m2 + ... + Z_mk

Доказательсто является элементарным упражнением на тему понятия фактор-кольца. Более содержательным является построение алгоритма, который по заданным модулям m₁, m₂, ..., m_k и остаткам r₁, r₂, ..., r_k вычисляет x. Этот алгоритм сводится несложным образом к расширенному алгоритму Евклида (вернее, к поиску обратного элемента в кольце вычетов по модулю m). Нам он, однако, не понадобится, и здесь мы его рассматривать не будем.

При всей своей простоте Китайская теорема об остатках очень важна, так как позволяет свести изучение кольца вычетов по модулю m, где m -- произвольное целое число, к изучению колец вычетов по модулю p^s, где p -- простое число. Действительно, любое целое число представимо в виде

m = m₁ m₂ ... m_k,

где каждое из чисел m_i есть степень простого числа. Тогда Z_m изоморфно прямой сумме колец Z_mi (примарных колец).

Пример на применение Китайской теоремы об остатках. Пусть m = 3 * 5 * 7 = 105. Сколько квадратных корней из единицы в Z_m?

Решение: 2³ = 8 корней. При представлении Z₁₀₅ в виде прямой суммы Z₃ + Z₅ + Z₇ они соответствуют тройкам

    (1, 1, 1)
    (1, 1, -1)
    (1, -1, 1)
    (1, -1, -1)
    (-1, 1, 1)
    (-1, 1, -1)
    (-1, -1, 1)
    (-1, -1, -1)

Определение. Функцией Эйлера φ(m) называется порядок группы обратимых элементов кольца Z_m.

Мы будем обозначать группу обратимых элементов кольца Z_m через U_m. Таким образом, φ(m) = |U_m|.

Другое, менее "математическое", определение функции Эйлера -- это количество целых чисел в интервале 0...m, взаимно простых с m.

Наша ближайшая цель -- вычисление функции Эйлера. Пусть

m = p₁^e₁ p₂^e₂ ... p_k^e_k

разложение m на простые множители.

Лемма 1. φ(m) = φ(p₁^e₁) φ(p₂^e₂) ... φ(p_k^e_k)

Доказательство. По Китайской теореме об остатках, кольцо Z_m изоморфно прямой сумме

Z_m = Z_p1^e1 + Z_p2^e2 + ... + Z_pk^ek

Отсюда вытекает, что группа обратимых элементов U_m изоморфна прямому произведению групп U_pi^ei:

U_m = U_p1^e1 × U_p2^e2 × ... × U_pk^ek

Утверждение леммы следует из того, что порядок прямого произведения групп равен произведению порядков сомножителей.

Лемма 2. Пусть m = p^k, где p -- простое число. Тогда

φ(m) = (p - 1) p^{k - 1}.

Доказательство. В кольце Z_p^k необратимы элементы, которые делятся на p, их число равно p^k / p = p^k-1. Следовательно, число обратимых элементов равно p^k - p^k-1 = (p - 1 ) p^{k - 1}.

Итак, мы получили формулу для функции Эйлера. Пусть

m = p₁^e₁ p₂^e₂ ... p_k^e_k

представление числа m в виде произведения степеней простых чисел. Тогда функция Эйлера от m выражается в виде

φ(m) = (p₁ - 1) p₁^{e₁ - 1} (p₂ - 1) p₂^{e₂ - 1} ... (p_k - 1) p_k^{e_k - 1}

Следствие 3. Пусть m = p q, где p и q -- простые числа. Тогда

φ(m) = (p - 1) (q - 1).

Обобщением малой теоремы Ферма является следующая теорема.

Теорема Эйлера. Пусть m -- произвольное целое число, и пусть x взаимно просто с m. Тогда

x^φ(m) == 1 (mod m)

В случае простого m получается малая теорема Ферма. Доказательство аналогично (следствие теоремы Лагранжа о порядке элементов конечной группы).

Следствие 4. Пусть m = p q. Тогда для всякого x (не обязательно взаимно простого с m) и для всякого целого числа h справедливо сравнение

x^{h φ(m) + 1} == x (mod m)

Следствие 4 вытекает из Китайской теоремы об остатках и аналогичного утверждения для простого m (форма малой теоремы Ферма: x^p == x (mod p)).

Мы сформулировали все утверждения, которые нам понадобятся при рассмотрении схемы кодирования с открытым ключом RSA.

Кодирование с открытым ключом, схема RSA

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

Итак, имеем два отображения:

E: S --> T
D: T --> S

где S -- множество всевозможных незашифрованных сообщений, T -- множество зашифрованных сообщений. Буква "E" -- первая буква слова "Encoding", буква "D" -- первая буква слова "Decoding". Отображение

E: s --> t

переводит исходное сообщение s в зашифрованное сообщение t, отображение

D: t --> s

переводит зашифрованное сообщение t обратно в s. Тот факт, что D является декодирующей процедурой, на математическом языке означает, что композиция отображений D E является тождественным отображением: для всякого s справедливо

D(E(s)) = s

или

D E = 1 (тождественное отображение в S).

Все это справедливо для любой схемы асимметричного кодирования. Перейдем непосредственно к схеме RSA, названной так по первым буквам фамилий ее авторов -- Rumley, Shamir, Adleman. Отметим сразу, что схема RSA обладает двумя дополнительными очень полезными свойствами.

1. Множество исходных сообщений S совпадает с множеством закодированных сообщений T; в качестве этого множества используется кольцо вычетов по модулю m, где m -- произведение двух больших простых чисел (десятичная запись m имеет длину не меньше 200).
2. Не только DE = 1, но и ED = 1! Таким образом, D и E -- два взаимно обратных отображения. Это позволяет владельцу секретной процедуры декодирования D применять ее для кодирования. При этом все могут раскодировать это сообщение, используя открытую процедуру E, но только владелец секретной процедуры D может послать его. Такая "обратная" схема применения открытого ключа позволяет удостоверить отправителя сообщения. В практических применениях (для аутентификации отправителя) обратная схема даже более важна, чем прямая.

Итак, в схеме RSA в качестве множества исходных и зашифрованных сообщений используется кольцо вычетов Z_m, где

m = p q

произведение двух больших простых чисел (длина десятичной записи каждого из чисел p и q не меньше 100). Всякое сообщение представляется в виде элемента из Z_m. (Любое ссобщение -- это последовательность битов, которую можно рассмотреть как большое целое число. Если длина сообщения больше, чем длина двоичной записи m, то оно разбивается на блоки, и каждый блок шифруется отдельно.)

Число m открытое, однако разложение m на множители -- секретное. Разложение позволяет вычислить функцию Эйлера (следствие 3):

φ(m) = (p - 1) (q - 1)

Нетрудно показать, что знание функции Эйлера дает возможность разложить число на множители, так что сложность задачи взламывания открытого ключа равна сложности задачи разложения на множители. Математики верят, что это действительно сложная задача, хотя никаких удовлетворительных оценок снизу в настоящее время не получено. (И вряд ли это NP-полная задача.)

Построение кодирующей процедуры E

E: s --> s^e (mod m)

Для практического вычисления применяется алгоритм быстрого возведения в степень.

Построение декодирующей процедуры D

e d == 1 (mod φ(m))

Это легко делается с помощью расширенного алгоритма Евклида. Пара (m, d) является секретным ключом. Отображение D состоит в возведении в степень d в кольце вычетов по модулю m.

D: t --> t^d (mod m)

Покажем, что отображение D является левым обратным к E, т.е. для всякого ссобщения s выполняется равенство D(E(s)) = s. Имеем

D(E(s)) == D(s^e) == (s^e) ^d == s^ed (mod m)

Так как e d == 1 (mod φ(m)), имеем

e d = 1 + h φ(m)

По следствию 4,

s^ed = s^{1 + h φ(m))} == s (mod m)

Итак, DE = 1. Аналогично доказывается, что ED = 1.

Суммируем все вышесказанное.

Рассматривается множество сообщений Z_m, где m -- произведение двух больших простых чисел: m = p q. Число m является открытым, но его разложение на множители -- секретным. Знание разложения позволяет вычислить функцию Эйлера φ(m) = (p - 1)(q - 1). Случайным образом выбирается обратимый элемент e в кольце вычетов по модулю φ(m). Для него вычисляется (с помощью расширенного алгоритма Евклида) обратный элемент d в кольце вычетов по модулю φ(m). Отображение E задается парой (m, e) и состоит в возведении в степень e по модулю m:

E(s) = s^e (mod m).

Отображение D задается парой (m, d) и состоит в возведении в степень d по модулю m:

D(t) = t^d (mod m).

Эти два отображения взаимно обратны. Пара (m, e) является открытым ключом (public key), пара (m, d) является секретным ключом (private key).

Пример. Рассмотрим пример с небольшими числами, чтобы только проиллюстрировать схему RSA. В реальных приложениях используют большие целые числа, порядка 200-400 десятичных цифр.

Пусть m = 11*13 = 143. Вычислим функцию Эйлера φ(m) = 10*12 = 120. Выберем e = 113, тогда d = 17 -- обратный к e элемент в кольце Z₁₂₀. Действительно,

113 * 17 = 1921 = 120 * 16 + 1.

Пара (143, 113) составляет открытый ключ, пара (143, 17) -- секретный ключ. Отображение E состоит в возведении в степень 113 по модулю 143, отображение D -- в степень 17 по модулю 143. Рассмотрим произвольное сообщение s = 123. Тогда

E(123) == 123¹¹³ (mod 143) == 41.

Таким образом, 41 -- это закодированное сообщение. Применим к нему декодирующую процедуру:

D(41) == 41¹⁷ (mod 143) == 123.

Мы получили исходное сообщение.

Алгоритмические задачи, связанные со схемой RSA

В связи со схемой RSA возникает ряд алгоритмических задач.

1. Для генерации ключей нам надо уметь генерировать большие простые числа. Близкой задачей является проверка простоты целого числа.

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

Опыт общения со студентами первого-второго курсов мехмата свидетельствует, что обычно никто не может предложить ничего более интересного, чем пробные деления на все нечетные числа (или все простые числа) до корня квадратного. Мы рассмотрим здесь несколько простых, но очень изящных алгоритмов проверки простоты и факторизации (разложения на множители): вероятностный тест простоты Рабина, алгоритмы факторизации Полларда (их называют также методами Монте-Карло).

Вероятностный тест простоты Рабина

Простейший тест проверки простоты числа m состоит в проверке малой теоремы Ферма. Выберем произвольное целое число b (например, b = 2), и возведем его в степень m - 1 по модулю m. Если мы получим не единицу, то по малой теореме Ферма число m составное. Беда состоит в том, что если

b^{m - 1} == 1 (mod m)

то ничего нельзя сказать об m. Древние греки ошибочно полагали, что все числа m, удовлетворяющие обращению малой теоремы Ферма для основания 2, простые: если

2^{m - 1} == 1 (mod m),

то m -- простое число. Минимальный контрпример к этому утверждению был найден только в XVII веке:

2³⁴⁰ == 1 (mod 341),

но число 341 -- не простое, 341 = 11 * 31. (Действительно, 2³⁴⁰ = (2¹⁰)³⁴ = 1024³⁴, но 1024 = 3 * 341 + 1 == 1 (mod 341), поэтому 1024³⁴ == 1 (mod 341).)

То, что 341 не удовлетворяет малой теореме Ферма, может быть показано с помощью других оснований:

3³⁴⁰ == 56 (mod 341)

Тем не менее существуют числа, которые не являются простыми, но которые ведут себя как простые в малой теореме Ферма. Такие числа называются кармайкловыми.

Определение. Число m называется кармайкловым, если оно не простое и для всякого b, взаимно простого с m, выполняется утверждение малой теоремы Ферма:

b^m-1 == 1 (mod m)

Минимальные кармайкловы числа -- это 561, 1105, 1729, ...

Несложно доказать следующее утверждение.

Предложение 5. Пусть

m = p₁^e₁ p₂^e₂ ... p_k^e₂ ---

представление целого числа m в виде произведения степеней простых. Число m является кармайкловым тогда и только тогда, когда

1) для всякого i показатель степени e_i = 1;
2) k >= 3;
3) для всякого i число p_i - 1 делит m - 1.

Доказательство. Докажем только обратную, наиболее интересную импликацию. Пусть число m удовлетворяет условиям 1-3. Рассмотрим произвольное b, взаимно простое с m. По Китайской теореме об остатках, кольцо Z_m представляется в виде прямой суммы

Z_m == Z_p1 + Z_p2 + ... + Z_pk.

При этом изоморфизме элемен b представляется в виде строки

b == (b₁, b₂, ..., b_k)

Тогда

b^(m-1) == (b₁^(m-1), b₂^(m-1), ..., b_k^(m-1)).

По малой теореме Ферма, для всякого i

b_i^(m-1) == 1 (mod p_i),

поскольку (m-1) делится на (p_i-1). Поэтому

b^(m-1) == (1, 1, ..., 1),

т.е. b^(m-1) == 1 (mod m).

Пример. Покажем, что число 561 является кармайкловым. Действительно, 561 = 3 * 11 * 17. Имеем

(3 - 1) | 560, (11 - 1) | 560, (17 - 1) | 560.

Следовательно, число 561 удовлетворяет условиям предложения 5.

Итак, для кармайкловых чисел тест простоты, основанный на теореме Ферма, не работает. Тем не менее его модификация, предложенная Рабином, применима к любым целым числам.

Тест Рабина является вероятностным. Это означает, что он использует датчик случайных чисел и, таким образом, работает не детерминированно. Для входного целого числа m тест Рабина может выдать один из следующих двух ответов.

1. Число m является составным.
2. Не знаю.

Таким образом, если мы применим 100 раз тест Рабина к числу m и получим 100 ответов "не знаю", то можно с большой вероятностью утверждать, что число m простое. Более точно, вероятность получения ста ответов "не знаю" для составного числа m не превышает (1/4)¹⁰⁰, т.е. практически равна нулю. Тем не менее тест Рабина не предъявляет доказательства того, что число m простое.

Перейдем непосредственно к изложению теста Рабина. Мы проверяем простоту входного числа m. Допустим сразу, что число m нечетное. (Существует только одно четное простое число -- 2.) Тогда число m-1 четное. Представим его в виде

m - 1 = 2^t * s

где s -- нечетное число. Выберем случайное число b такое, что b ≠ 0, b ≠ 1 (mod m) При выборе b используется датчик случайных чисел.

Используя алгоритм быстрого возведения в степень по модулю m, вычислим следующую последовательность элементов кольца Z_m:

x₀ == b^s (mod m), (1)
x₁ == x₀ x₀ (mod m),
x₂ == x₁ x₁ (mod m),
...
x_t == x_{t - 1} x_{t - 1} == b^{m - 1} (mod m)

(На каждом шаге мы возводим в квадрат число, полученное на предыдущем шаге.) Тест Рабина выдает ответ "m -- составное число" в случае, если

• x_t ≠ 1 (mod m), или
  
• в последовательности x₀, x₁, x₂, ..., x_t имеется фрагмент вида
  ..., *, 1, ...
  где звездочкой обозначено число, отличное от единицы или минус единицы по модулю m.

Теорема 6 (законность теста Рабина).

1. Если тест Рабина выдает ответ "m -- составное число", то m действительно является составным.
2. Вероятность ответа "не знаю" для составного числа m не превосходит 1/4.

Доказательство. Докажем только первое утверждение. Если x_t ≠ 1 (mod m), то m не удовлетворяет малой теореме Ферма и, следовательно, не является простым. Если же последовательность (1) содержит фрагмент ..., a, 1, ..., где a ≠ +- 1 (mod m), то имеем

a² == 1 (mod m), a ≠ 1, a ≠ -1 (mod m)

Если бы m было простым, то кольцо Z_m являлось бы полем. Но в любом поле есть только два квадратных корня из единицы: это единица и минус единица. (По теореме Безу, число корней многочлена не превосходит его степени, квадратные корни из единицы -- это корни многочлена x² - 1.) Следовательно, число m не является простым.

Алгоритмы факторизации целых чисел

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

Пример. В средние века существовала гипотеза, что все числа Ферма

F_k = 2^2k + 1

простые. Действительно,

F₁ = 5, F₂ = 17, F₃ = 257, F₄ = 65537 --

простые числа. Опроверг эту гипотезу Эйлер, разложив на множители число

F₅ = 2³² + 1 = 4294967297.

Попробуйте сделать то же самое, не используя компьютер (как в свое время Эйлер).

Указание. Идея Эйлера проста -- достаточно найти два разных числа, квадраты которых совпадают по модулю m:

a² == b² (mod m), a ≠ +-b (mod m)

Тогда наибольший общий делитель чисел (a - b) и m нетривиален. Действительно,

0 == a² - b² == (a - b)(a + b) (mod m)

Произведение (a - b)(a + b) делится на m, но сомножители не делятся. Следовательно, НОД(a - b, m) нетривиален.

Алгоритмы факторизации Полларда (методы Монте-Карло)

Мы рассмотрим 2 весьма изящных метода факторизации, предложенных Поллардом. Они очень просты и позволяют быстро извлечь из составного числа все его небольшие простые делители.

Метод Монте-Карло 1

Время работы этого метода порядка корень квадратный из минимального простого числа, делящего m. То есть в худшем случае, когда m есть произведение двух простых чисел примерно одного порядка, число m раскладывается этим методом на множители за время корень четвертой степени из m. Алгоритм очень прост, его запись намного короче, чем объяснение, почему он работает. Слова "Монте-Карло" присутствуют в его названии потому, что работа алгоритма зависит от случайного выбора начального числа.

Рассмотрим отображение f кольца Z_m в Z_m:

f: Z_m --> Z_m
f(x) = x² + 1 (mod m)

Выберем случайным образом элемент b₀ кольца Z_m. Рассмотрим бесконечную последовательность элементов кольца Z_m:

b₀, b₁ = f(b₀), b₂ = f(b₁), b₃ = f(b₂), ... (2)

Последовательность представляет собой орбиту элемента b₀ при отображении f. Поскольку все элементы последовательности принадлежат конечному множеству, последовательность циклическая -- точнее, она содержит начальный апериодический отрезок и далее бесконечно повторяющийся период.

Пусть p -- делитель числа m. Рассмотрим элементы последовательности (2) по модулю p (т.е. образы элементов b_i при каноническом эпиморфизме Z_m --> Z_p):

c₀ == b₀(mod p), c₁ == b₁(mod p), c₂ == b₂(mod p), ... (3)

Так как в Z_p меньше элементов, чем в Z_m, то с большой вероятностью период последовательности (3) меньше, чем период последовательности (2). Следовательно, найдется пара индексов i, j такая, что

c_i = c_j, b_i ≠ b_j.

Это означает, что

b_i == b_j(mod p), b_i ≠ b_j(mod m).

Отсюда вытекает, что (b_i - b_j) делится на p, но не делится на m. Следовательно, НОД(b_i - b_j, m) нетривиален, и нам удалось разложить m на множители.

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

Можно предложить 2 способа решения задачи поиска цикла в последовательности. Первый способ наиболее простой. Второй чуть-чуть сложнее, но зато более быстрый.

Способ 1

Выполняется следующая последовательность сравнений:

b₀ <---> b₁
b₁ <---> b₃
b₂ <---> b₅
b₃ <---> b₇
b₄ <---> b₉
. . .
b_i <---> b_2i+1
. . .

Выпишем алгоритм нахождения делителя.

алгоритм факторизация1(вход: целое число m, 
|                      выход: целое число d): успех
| дано: целое число m
| надо: получить нетривиальный делитель d числа m
| возвращаемое значение: true, если удалось разложить,
|                        false в противном случае
начало
| maxSteps := 1000000                 // Максимальное число шагов
| step := 0
|
| b0 := случайное число в интервале 0..m
| b1 := mod(b0 * b0 + 1, m)
| d := gcd(b1 - b0, m)
|
| цикл пока step < maxSteps && d == 1 // Пока НОД тривиален
| |выполнять
| | b0 = mod(b0 * b0 + 1, m)          // Применяем отображение f
| | b1 = mod(b1 * b1 + 1, m);         // один раз к b0 и дважды
| | b1 = mod(b1 * b1 + 1, m)          // к b1
| | d := gcd(b1 - b0, m)
| | step := step + 1
| конец_цикла
|
| вернуть (d != 1)                    // Успех := d нетривиален
конец_алгоритма

Способ 2

Выполняется следующая бесконечная последовательность сравнений

b₀ <---> b₁

b₁ <---> b₂
b₁ <---> b₃

b₂ <---> b₄
b₂ <---> b₅
b₂ <---> b₆
b₂ <---> b₇

b₄ <---> b₈
b₄ <---> b₉
b₄ <---> b₁₀
. . .
b₄ <---> b₁₅

b₈ <---> b₁₆
b₈ <---> b₁₇
. . .
b₈ <---> b₃₁

. . .

Вся последовательность сравнений разбивается на серии. В очередной серии мы сравниваем элемент b_s, где s -- степень двойки, последовательно с элементами b_2s, b_2s+1, b_2s+2, ..., b_4s-1. Серия содержит 2s сравнений.

Выпишем алгоритм.

алгоритм факторизация2(вход: целое число m, 
|                      выход: целое число d): успех
| дано: целое число m
| надо: получить нетривиальный делитель d числа m
| возвращаемое значение: true, если удалось разложить,
|                        false в противном случае
начало
| maxSteps := 19            // Максимальная длина серии 2^19
| step := 0
|
| b0 := случайное число в интервале 0..m
| b1 := mod(b0 * b0 + 1, m)
| a := b1                   // Первый элемент серии
| seriesLength := 1         // Длина серии
| d := gcd(b1 - b0, m)
|
| цикл пока step < maxSteps && d == 1   // пока НОД тривиален
| | выполнять
| | Инвариант: 
| |     b0 -- элемент последовательности с индексом,
| |           равным нулю или степени двойки
| |     a  -- элемент, индекс которого равен удвоенному индексу
| |           элемента b0 (или 1, если индекс b0 равен 0)
| |     seriesLength == удвоенному индексу элемента a
| | d := gcd(b1 - b0, m)
| | len := 0
| |
| | цикл пока d == 1 и len < seriesLength 
| | | выполнять
| | | b1 = mod(b1 * b1 + 1, m);
| | | d := gcd(b1 - b0, m)
| | | len := len + 1
| | конец_цикла
| |
| | b0 := a
| | a := b1
| | seriesLength := seriesLength * 2
| конец_цикла
|
| вернуть (d != 1)          // Успех := d нетривиален
конец_алгоритма

Метод Монте-Карло 2

Пусть m -- целое число, которое мы раскладываем на множители. Оно представимо в виде произведения степеней простых чисел

m = p₁^e₁ p₂^e₂ ... p_k^e_k

Предположим, что p₁-1 представимо в виде произведения степеней простых чисел, причем каждая из этих степеней не очень велика. Более точно, существует N такое, что

p₁-1 = q₁^a₁ q₂^a₂ ... q_r^a_r,
q₁^a₁ < N, q₂^a₂ < N, ..., q_r^a_r < N.

Рассмотрим всевозможные максимальные степени простых чисел, не превосходящие N. Например, пусть N = 20, тогда рассматриваются степени простых 16, 9, 5, 7, 11, 13, 17, 19. Обозначим эти степени простых через t₁, t₂, ..., t_s. Выберем произвольное целое число b. Рассмотрим последовательность

b₀ = b, b₁ = b₀^t₁ (mod m),
b₂ = b₁^t₂ (mod m), ..., b_s = b_s-1^t_s (mod m)

Каждый раз, вычислив b_i, вычисляем одновременно

НОД(b_i - 1, m).

Утверждается, что с большой вероятностью на каком-то шаге этот НОД будет нетривиальным делителем N. Действительно, покажем, что

p | НОД(b_s - 1, m).

Действительно,

b_s = b^{t₁ t₂ ... t_s}

и, поскольку по предположению, p₁ - 1 | t₁ t₂ ... t_s, то есть t₁ t₂ ... t_s = (p₁ - 1)g, то

b_s = b^{t₁ t₂ ... t_s} = b^{(p₁ - 1)g} = (b^{(p1 - 1)})^g = 1 (mod p₁)

по малой теореме Ферма. Значит, b_s - 1 делится на p₁, число m также делится на p₁, следовательно, НОД(b_s - 1, m) делится на p₁.

Проиллюстрируем алгоритм на простом примере. Возьмем N = 20. Выпишем все степени простых, не превосходящие 20:

t₁ = 16, t₂ = 9, t₃ = 5, t₄ = 7,
t₅ = 11, t₆ = 13, t₇ = 17, t₈ = 19.

Попытаемся разложить на множители число m = 41779 = 41 * 1019. Выберем b = 2. Последовательно вычисляем

2¹⁶ (mod 41779) = 23757, gcd(23757 - 1, 41779) = 1,
23757⁹ (mod 41779) = 7970, gcd(7970 - 1, 41779) = 1,
7970⁵ (mod 41779) = 33580, gcd(33580 - 1, 41779) = 41.

Мы получили нетривиальный делитель на третьем шаге, поскольку 41-1 = 8*5 делит (t₁ t₂ t₃) = 16 * 9 * 5.

Мощность алгоритма зависит от числа N -- чем больше оно, тем большие числа можно разложить с помощью этого алгоритма. Работа алгоритма разбивается на 2 шага. Сначала мы генерируем все максимальные степени простых чисел, не превосходящие N. Этот шаг выполняется только один раз и не зависит от входного числа m, поэтому сгенерированные степени можно, к примеру, записать в файл и в дальнейшем использовать многократно. Затем мы выбираем случайным образом число b и вычисляем указанную выше последовательность степеней b. Для каждой степени b_i вычисляется НОД(b_i - 1, m). Алгоритм завершается успешно, если вычисленный НОД нетривиален. Алгоритм можно убыстрить, если вычислять НОД не на каждом шаге, а, скажем, на каждом сотом шаге. При этом на промежуточных шагах последовательно вычисляется произведение

(b_i - 1)(b_i+1 - 1)(b_i+2 - 1)... (b_i+99 - 1) = n (mod m)

и затем вычисляется НОД(n, m).