Ро-алгоритм Полларда

Эта статья — об алгоритме факторизации чисел. О методе дискретного логарифмирования см. ро-метод Полларда для дискретного логарифмирования.

Числовая последовательность зацикливается, начиная с некоторого n. Цикл может быть представлен в виде греческой буквы ρ.

Ро-алгоритм ([math]\displaystyle{ \rho }[/math]-алгоритм) — предложенный Джоном Поллардом в 1975 году алгоритм, служащий для факторизации (разложения на множители) целых чисел. Данный алгоритм основывается на алгоритме Флойда поиска длины цикла в последовательности и некоторых следствиях из парадокса дней рождения. Алгоритм наиболее эффективен при факторизации составных чисел с достаточно малыми множителями в разложении. Сложность алгоритма оценивается как [math]\displaystyle{ O(N^{1/4}) }[/math].

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

История алгоритма

В конце 60-х годов XX века Роберт Флойд придумал достаточно эффективный метод решения задачи нахождения цикла, также известный, как алгоритм «черепаха и заяц». Джон Поллард, Дональд Кнут и другие математики проанализировали поведение этого алгоритма в среднем случае. Было предложено несколько модификаций и улучшений алгоритма.

В 1975 году Поллард опубликовал статью, в которой он, основываясь на алгоритме Флойда обнаружения циклов, изложил идею алгоритма факторизации чисел, работающего за время, пропорциональное [math]\displaystyle{ N^{1/4} }[/math]. Автор алгоритма назвал его методом факторизации Монте-Карло, отражая кажущуюся случайность чисел, генерируемых в процессе вычисления. Однако позже метод всё-таки получил своё современное название — ρ-aлгоритм Полларда.

В 1981 году Ричард Брент и Джон Поллард с помощью алгоритма нашли наименьшие делители чисел Ферма [math]\displaystyle{ F_{n} = 2^{2^n}+1 }[/math] при [math]\displaystyle{ 5 \leq n \leq 13 }[/math]. Скорость алгоритма сильно зависит лишь от величины наименьшего делителя исходного числа, но не от самого числа. Так, поиск наименьшего делителя седьмого числа Ферма — [math]\displaystyle{ \begin{array}{lll}F_7= 340282366920938463463374607431768211457 = 59\,649\,589\,127\,497\,217 \cdot 5\,704\,689\,200\,685\,129\,054\,721;\end{array} }[/math], занимает гораздо больше времени, чем поиск делителя двенадцатого числа Ферма (т.к. его делитель 114689 значительно меньше, хотя само число состоит более чем из 1200 десятичных цифр).

В рамках проекта «Cunningham project» алгоритм Полларда помог найти делитель длиной 19 цифр числа [math]\displaystyle{ 2^{2386}+1 }[/math]. Большие делители также могли бы быть найдены, однако открытие метода факторизации с помощью эллиптических кривых сделало алгоритм Полларда неконкурентоспособным.

Описание алгоритма

Оригинальная версия

Рассматривается последовательность целых чисел [math]\displaystyle{ {x_{n}} }[/math], такая что [math]\displaystyle{ x_{0}=2 }[/math] и [math]\displaystyle{ x_{i+1}=(x_{i}^{2}-1\,) (\mathrm{mod}\, N) }[/math], где [math]\displaystyle{ N }[/math] — число, которое нужно факторизовать. Оригинальный алгоритм выглядит следующим образом:

1. Вычисляются тройки чисел

[math]\displaystyle{ (x_{i}, x_{2i}, Q_{i}), i=1,2,... }[/math], где [math]\displaystyle{ Q_{i} \equiv \prod_{j=1}^{i}(x_{2j}-x_{j})\, (\mathrm{mod}\, N) }[/math].

Причём каждая такая тройка получается из предыдущей.

2. Каждый раз, когда число [math]\displaystyle{ i }[/math] кратно числу [math]\displaystyle{ m }[/math] (скажем, [math]\displaystyle{ m=100 }[/math]), вычисляется наибольший общий делитель [math]\displaystyle{ d_{i}=\mathrm{GCD}(Q_{i},N) }[/math] любым известным методом.

3. Если [math]\displaystyle{ 1 \lt d_{i} \lt N }[/math], то частичное разложение числа [math]\displaystyle{ N }[/math] найдено, причём [math]\displaystyle{ N = d_{i} \times (N/d_{i}) }[/math].

Найденный делитель [math]\displaystyle{ d_{i} }[/math] может быть составным, поэтому его также необходимо факторизовать. Если число [math]\displaystyle{ N/d_{i} }[/math] составное, то продолжаем алгоритм с модулем [math]\displaystyle{ N' = N/d_{i} }[/math].

4. Вычисления повторяются [math]\displaystyle{ S }[/math] раз. Если при этом число не было до конца факторизовано, выбирается, например, другое начальное число [math]\displaystyle{ x_{0} }[/math].

Современная версия

Пусть [math]\displaystyle{ N }[/math] составное целое положительное число, которое требуется разложить на множители. Алгоритм выглядит следующим образом:

1. Случайным образом выбирается небольшое число [math]\displaystyle{ x_{0} }[/math] и строится последовательность [math]\displaystyle{ \{x_{n}\}, n = 0, 1, 2, ... }[/math], определяя каждое следующее как [math]\displaystyle{ x_{n+1} = F(x_{n})\, (\mathrm{mod}\,\, N) }[/math].
2. Одновременно на каждом i-ом шаге вычисляется [math]\displaystyle{ d = \mathrm{GCD}(N,|x_{i}-x_{j}|) }[/math] для каких-либо [math]\displaystyle{ i }[/math], [math]\displaystyle{ j }[/math] таких, что [math]\displaystyle{ j\lt i }[/math], например, [math]\displaystyle{ i = 2j }[/math].
3. Если [math]\displaystyle{ d\gt 1 }[/math], то вычисление заканчивается, и найденное на предыдущем шаге число [math]\displaystyle{ d }[/math] является делителем [math]\displaystyle{ N }[/math]. Если [math]\displaystyle{ N/d }[/math] не является простым числом, то процедуру поиска делителей продолжается, взяв в качестве [math]\displaystyle{ N }[/math] число [math]\displaystyle{ N'=N/d }[/math].

На практике функция [math]\displaystyle{ F(x) }[/math] выбирается не слишком сложной для вычисления (но в то же время не линейным многочленом), при условии того, что она не должна порождать взаимно однозначное отображение. Обычно в качестве [math]\displaystyle{ F(x) }[/math] выбираются функции [math]\displaystyle{ F(x) = x^2 \pm 1 (\mathrm{mod}\, N) }[/math] или [math]\displaystyle{ F(x) = x^2 \pm a (\mathrm{mod}\, N) }[/math]. Однако функции [math]\displaystyle{ x^2-2 }[/math] и [math]\displaystyle{ x^2 }[/math] не подходят.

Если известно, что для делителя [math]\displaystyle{ p }[/math] числа [math]\displaystyle{ N }[/math] справедливо [math]\displaystyle{ p \equiv 1\, (\mathrm{mod}\, k) }[/math] при некотором [math]\displaystyle{ k \gt 2 }[/math], то имеет смысл использовать [math]\displaystyle{ F(x) = x^k + b }[/math].

Существенным недостатком алгоритма в такой реализации является необходимость хранить большое число предыдущих значений [math]\displaystyle{ x_{j} }[/math].

Улучшения алгоритма

Изначальная версия алгоритма обладает рядом недостатков. В настоящий момент существует несколько подходов к улучшению оригинального алгоритма.

Пусть [math]\displaystyle{ F(x)=(x^2-1)\bmod N }[/math]. Тогда, если [math]\displaystyle{ (x_j-x_i)\equiv0\pmod p }[/math], то [math]\displaystyle{ (F(x_j)-F(x_i))\equiv0\pmod p }[/math], поэтому, если пара [math]\displaystyle{ (x_i,x_j) }[/math] даёт решение, то решение даст любая пара [math]\displaystyle{ (x_{i+k},x_{j+k}) }[/math].

Поэтому нет необходимости проверять все пары [math]\displaystyle{ (x_i,x_j) }[/math], а можно ограничиться парами вида [math]\displaystyle{ (x_i,x_j) }[/math], где [math]\displaystyle{ j=2^k }[/math], и [math]\displaystyle{ k }[/math] пробегает набор последовательных значений 1, 2, 3, …, а [math]\displaystyle{ i }[/math] принимает значения из интервала [math]\displaystyle{ [2^k+1;2^{k+1}] }[/math]. Например, [math]\displaystyle{ k=3 }[/math], [math]\displaystyle{ j=2^3=8 }[/math], а [math]\displaystyle{ i\in[9;16] }[/math].

Эта идея была предложена Ричардом Брентом в 1980 году и позволяет уменьшить количество выполняемых операций приблизительно на 25 %.

Ещё одна вариация ρ-алгоритма Полларда была разработана Флойдом. Согласно Флойду, значение [math]\displaystyle{ y }[/math] обновляется на каждом шаге по формуле [math]\displaystyle{ y=F^2(y)=F(F(y)) }[/math], поэтому на шаге [math]\displaystyle{ i }[/math] будут получены значения [math]\displaystyle{ x_i=F^i(x_0) }[/math], [math]\displaystyle{ y_i=x_{2i}=F^{2i}(x_0) }[/math], и НОД на этом шаге вычисляется для [math]\displaystyle{ N }[/math] и [math]\displaystyle{ y-x }[/math].

Пример факторизации числа

Данный пример наглядно демонстрирует ρ-алгоритм факторизации (версия алгоритма, с улучшением Флойда), для числа N = 8051:

Таблица: факторизация числа 8051

n = 8051, F(x) = (x2 + 1) mod n , x0 = y0 = 2
i	xi=F(xi-1)	yi=F(F(yi-1))	НОД(|xi − yi|, 8051)
1	5	26	1
2	26	7474	1
3	677	871	97

Используя другие варианты полинома [math]\displaystyle{ F(x) }[/math], можно также получить делитель 83:

Таблица: факторизация числа 8051

n = 8051, F(x) = (x2 + 3) mod n , x0 = y0 = 2
i	xi=F(xi-1)	yi=F(F(yi-1))	НОД(|xi − yi|, 8051)
1	7	52	1
2	52	1442	1
3	2707	778	1
4	1442	3932	83

Таким образом, d₁ = 97, d₂ = 83 — нетривиальные делители числа 8051.

После нахождения делителя числа, в ρ-алгоритме предлагается продолжать вычисления и искать делители числа [math]\displaystyle{ N/d }[/math], если [math]\displaystyle{ N/d }[/math] не является простым. В этом простом примере данного шага совершать не потребовалось.

Обоснование ρ-алгоритма Полларда

Алгоритм основывается на известном парадоксе дней рождения.

Парадокс дней рождений, кратко: Пусть [math]\displaystyle{ \lambda \gt 0 }[/math]. Для случайной выборки из [math]\displaystyle{ l+1 }[/math] элементов, каждый из которых меньше [math]\displaystyle{ q }[/math], где [math]\displaystyle{ l = \sqrt{2 \lambda q} }[/math], вероятность того, что два элемента окажутся одинаковыми [math]\displaystyle{ p \gt 1 - e^{-\lambda} }[/math].

Следует отметить, что вероятность [math]\displaystyle{ p = 0.5 }[/math] в парадоксе дней рождения достигается при [math]\displaystyle{ \lambda \approx 0.69 }[/math].

Пусть последовательность [math]\displaystyle{ \{u_{n}\} }[/math] состоит из разностей [math]\displaystyle{ x_{i} - x_{j} }[/math], проверяемых в ходе работы алгоритма. Определяется новая последовательность [math]\displaystyle{ \{z_{n}\} }[/math], где [math]\displaystyle{ z_{n} = u_{n}\, \mathrm{mod}\, q }[/math], [math]\displaystyle{ q }[/math] — меньший из делителей числа [math]\displaystyle{ N }[/math].

Все члены последовательности [math]\displaystyle{ \{z_{n}\} }[/math] меньше [math]\displaystyle{ \sqrt{N} }[/math]. Если рассматривать её как случайную последовательность целых чисел, меньших [math]\displaystyle{ q }[/math], то, согласно парадоксу дней рождения, вероятность того, что среди [math]\displaystyle{ l+1 }[/math] её членов попадутся два одинаковых, превысит [math]\displaystyle{ 1/2 }[/math] при [math]\displaystyle{ \lambda \approx 0.69 }[/math], тогда [math]\displaystyle{ l }[/math] должно быть не меньше [math]\displaystyle{ \sqrt{2\lambda q} \approx \sqrt{1.4 q} \approx 1.18\sqrt{q} }[/math].

Если [math]\displaystyle{ z_{i}=z_{j} }[/math], тогда [math]\displaystyle{ x_{i}-x_{j} \equiv 0\, \mathrm{mod}\, q }[/math], то есть, [math]\displaystyle{ x_{i}-x_{j}=kq }[/math] для некоторого целого [math]\displaystyle{ k }[/math]. Если [math]\displaystyle{ x_{i} \neq x_{j} }[/math], что выполняется с большой вероятностью, то искомый делитель [math]\displaystyle{ q }[/math] числа [math]\displaystyle{ N }[/math] будет найден как [math]\displaystyle{ \mathrm{GCD}(N, |x_{i}-x_{j}|) }[/math]. Поскольку [math]\displaystyle{ \sqrt{q} \leq n^{1/4} }[/math], то с вероятностью, превышающей [math]\displaystyle{ 1/2 }[/math] , делитель [math]\displaystyle{ N }[/math] будет найден за [math]\displaystyle{ 1.18 \times N^{1/4} }[/math] итераций.

Сложность алгоритма

Чтобы оценить сложность алгоритма, рассматривается последовательность, строящаяся в процессе вычислений, как случайная (разумеется, ни о какой строгости при этом говорить нельзя). Чтобы полностью факторизовать число [math]\displaystyle{ N }[/math] длиной [math]\displaystyle{ \beta }[/math] бит, достаточно найти все его делители, не превосходящие [math]\displaystyle{ \sqrt{N} }[/math], что требует максимум порядка [math]\displaystyle{ \sqrt{N} }[/math] арифметических операций, или [math]\displaystyle{ N^{1/4}\beta^2 = 2^{\beta/4}\beta^2 }[/math] битовых операций.

Поэтому сложность алгоритма оценивается, как [math]\displaystyle{ O(N^{1/4}) }[/math]. Однако в этой оценке не учитываются накладные расходы по вычислению наибольшего общего делителя. Полученная сложность алгоритма, хотя и не является точной, достаточно хорошо согласуется с практикой.

Справедливо следующее утверждение: пусть [math]\displaystyle{ N }[/math] — составное число. Тогда существует такая константа [math]\displaystyle{ C }[/math], что для любого положительного числа [math]\displaystyle{ \lambda }[/math] вероятность события, состоящего в том, что ρ-алгоритм Полларда не найдет нетривиального делителя [math]\displaystyle{ N }[/math] за время [math]\displaystyle{ C\sqrt{\lambda \sqrt N}(\log N)^2 }[/math], не превосходит величины [math]\displaystyle{ e^{-\lambda} }[/math]. Данное утверждение следует из парадокса дней рождения.

Особенности реализации

Объём памяти, используемый алгоритмом, можно значительно уменьшить.

 int Rho-Поллард (int N)
 { 
   int x = random(1, N-2);
   int y = 1; int i = 0; int stage = 2;
   while (Н.О.Д.(N, abs(x - y)) == 1)
   {
     if (i == stage){
       y = x;
       stage = stage*2; 
     }
     x = (x*x + 1) (mod N);
     i = i + 1;
   }
   return Н.О.Д(N, abs(x-y));
 }

В этом варианте вычисление требует хранить в памяти всего три переменные [math]\displaystyle{ N }[/math], [math]\displaystyle{ x }[/math], и [math]\displaystyle{ y }[/math], что выгодно отличает алгоритм в такой реализации от других методов факторизации чисел.

Распараллеливание алгоритма

Алгоритм Полларда допускает распараллеливание с использованием как систем с разделяемой памятью, так и систем с распределенной памятью (передача сообщений), однако второй случай является наиболее интересным с практической точки зрения.

Система с распределенной памятью

Существующий метод распараллеливания заключается в том, что каждый вычислительный узел исполняет один и тот же последовательный алгоритм, однако, исходное число [math]\displaystyle{ x_0 }[/math] и/или полином [math]\displaystyle{ F(x) }[/math] берутся различными. Для упрощения распараллеливания, предлагается получать их из генератора случайных чисел. Однако такая параллельная реализация не даёт линейного ускорения.

Предположим что есть [math]\displaystyle{ P }[/math] одинаковых исполнителей. Если мы используем [math]\displaystyle{ P }[/math] различных последовательностей (то есть различных полиномов [math]\displaystyle{ F(x) }[/math]), то вероятность того, что первые [math]\displaystyle{ k }[/math] чисел в этих последовательностях будут различными по модулю [math]\displaystyle{ p }[/math], будет примерно равна [math]\displaystyle{ \exp({-k^2 P}/2p) }[/math]. Таким образом, максимальное ускорение можно оценить как [math]\displaystyle{ P^{1/2} }[/math].

Ричард Крэндалл предположил, что достижимо ускорение [math]\displaystyle{ O(P/(\log P)^2) }[/math], однако данное утверждение пока не проверено.

Система с общей памятью

Предыдущий метод, очевидно, можно использовать и на системах с общей памятью, однако, гораздо разумнее использовать единый генератор [math]\displaystyle{ F(x) }[/math].