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

ρ-алгоритм Полларда — алгоритм факторизації цілих чисел, що ґрунтується на пошуку циклу в послідовності і деяких наслідках із парадоксу днів народжень. Його запропонував Джон Поллард^[en] (1975). Алгоритм найбільш ефективний для факторизації складених чисел із досить малими множниками в розкладі. Обчислювальна складність оцінюється як ${\displaystyle O(N^{1/4})}$.

Не плутати з алгоритмом (p-1), який також розробив Джон Поллард.

Числова послідовність зациклюється, починаючи з деякого n. Цикл виглядає як грецька літера ρ.

У всіх варіантах ρ-алгоритму Полларда будується числова послідовність, елементи якої, починаючи з деякого номера n, утворюють цикл, що можна проілюструвати розташуванням членів послідовності у вигляді грецької літери ρ. Це й послужило назвою для сімейства методів.

Історія алгоритму

Наприкінці 60-х років XX століття Дональд Кнут опублікував досить ефективний алгоритм пошуку циклу в послідовності, також відомий, як алгоритм «черепаха та заєць», який він пов'язував з ім'ям Флойда^[1]. Джон Поллард^[en], Дональд Кнут та інші математики проаналізували поведінку цього алгоритму в середньому випадку. Було запропоновано кілька модифікацій та поліпшень алгоритму.

У 1975 році Поллард опублікував статтю, в якій він, ґрунтуючись на алгоритмі Флойда виявлення циклів, виклав ідею алгоритму факторизації чисел, який виконується за час, пропорційний ${\displaystyle N^{1/4}}$^[2]. Автор назвав його методом факторизації Монте-Карло, тому, що в процесі обчислення генерується псевдовипадкова послідовність чисел. Проте пізніше метод все-таки назвали ρ-алгоритмом Полларда^[3].

У 1981 році Річард Брент^[ru] і Джон Поллард за допомогою цього алгоритму знайшли менший дільник восьмого числа Ферма ${\displaystyle F_{8}=2^{2^{8}}+1}$^[4].

Так, ${\displaystyle F_{8}=1238926361552897\cdot p_{62}}$, де ${\displaystyle p_{62}}$ — просте число, що складається з 62 десяткових цифр.

У межах проекту «Cunningham project^[en]» алгоритм Полларда допоміг знайти дільник числа ${\displaystyle 2^{2386}+1}$ довжиною 19 цифр. Більші дільники також можна б знайти, але відкриття методу факторизації за допомогою еліптичних кривих^[ru] зробило алгоритм Полларда неконкурентоспроможним^[5].

Опис алгоритму

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

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

1. Будемо обчислювати трійки чисел

${\displaystyle (x_{i},x_{2i},Q_{i}),i=1,2,...}$, де ${\displaystyle Q_{i}\equiv \prod _{j=1}^{i}(x_{2j}-x_{j})\,(\mathrm {mod} \,N)}$.

Причому кожна така трійка отримується з попередньої.

2. Щоразу, коли число ${\displaystyle i}$ кратне числу ${\displaystyle m}$ (скажімо, ${\displaystyle m=100}$), будемо обчислювати найбільший спільний дільник ${\displaystyle d_{i}=\mathrm {GCD} (Q_{i},N)}$ будь-яким відомим методом.

3. Якщо ${\displaystyle 1<d_{i}<N}$, то знайдено часткове розкладання числа ${\displaystyle N}$, причому ${\displaystyle N=d_{i}\times (N/d_{i})}$.

Знайдений дільник ${\displaystyle d_{i}}$ може бути складовим, тому його також необхідно факторизувати. Якщо число ${\displaystyle N/d_{i}}$ складене, то продовжуємо алгоритм з модулем ${\displaystyle N'=N/d_{i}}$.

4. Обчислення повторюються ${\displaystyle S}$ раз. Наприклад, можна зупинити алгоритм при ${\displaystyle i=S=10^{5}}$. Якщо при цьому число не було до кінця факторизовано, можна вибрати, наприклад, інше початкове число ${\displaystyle x_{0}}$.

Сучасна версія

Нехай ${\displaystyle N}$ складене ціле додатне число, яке потрібно розкласти на множники. Алгоритм виглядає таким чином:^[7]

1. Вибираємо невелике число ${\displaystyle x_{0}}$ та будуємо послідовність ${\displaystyle \{x_{n}\},n=0,1,2,...}$, визначаючи кожне наступне як ${\displaystyle x_{n+1}=F(x_{n})\,(\mathrm {mod} \,\,N)}$.
2. Одночасно на кожному i-ому кроці обчислюємо ${\displaystyle d=\mathrm {GCD} (N,|x_{i}-x_{j}|)}$ для будь-яких ${\displaystyle i}$, ${\displaystyle j}$ таких, що ${\displaystyle j<i}$, наприклад, ${\displaystyle i=2j}$.
3. Якщо виявили, що ${\displaystyle d>1}$, то обчислення закінчується, і знайдене на попередньому кроці число ${\displaystyle d}$ є дільником ${\displaystyle N}$. Якщо ${\displaystyle N/d}$ не є простим числом, то процедуру пошуку дільників можна продовжити, узявши як ${\displaystyle N}$ число ${\displaystyle N'=N/d}$.

Як на практиці обирати функцію ${\displaystyle F(x)}$? Функція має бути не надто складною для обчислення, але в той же час не має бути лінійним многочленом, а також не повинна породжувати взаємно однозначне відображення. Зазвичай за ${\displaystyle F(x)}$ беруть функцію ${\displaystyle F(x)=x^{2}\pm 1(\mathrm {mod} \,N)}$ або ${\displaystyle F(x)=x^{2}\pm a(\mathrm {mod} \,N)}$^[8]. Однак не слід застосовувати функції ${\displaystyle x^{2}-2}$ та ${\displaystyle x^{2}}$^[6].

Якщо відомо, що для дільника ${\displaystyle p}$ числа ${\displaystyle N}$ справедливо ${\displaystyle p\equiv 1\,(\mathrm {mod} \,k)}$ при деякому ${\displaystyle k>2}$, то має сенс застосувати ${\displaystyle F(x)=x^{k}+b}$^[6].

Істотним недоліком алгоритму в такий реалізації є необхідність зберігати велику кількість попередніх значень ${\displaystyle x_{j}}$.

Покращення алгоритму

Початкова версія алгоритму має низку недоліків. На даний момент^[коли?] існує кілька підходів до поліпшення оригінального методу.

Нехай ${\displaystyle F(x)=(x^{2}-1)\mathrm {mod} \,N}$. Зауважимо, що й ${\displaystyle (x_{j}-x_{i})\equiv 0(\mathrm {mod} \,p)}$, то ${\displaystyle (f(x_{j})-f(x_{i}))\equiv 0(\mathrm {mod} \,p)}$, тому, якщо пара ${\displaystyle (x_{i},x_{j})}$ дає нам розв'язок, то розв'язок дасть будь-яка пара ${\displaystyle (x_{i+k},x_{j+k})}$.

Тому, немає потреби перевіряти всі пари ${\displaystyle (x_{i},x_{j})}$, а можна обмежитися парами виду ${\displaystyle (x_{i},x_{j})}$, де ${\displaystyle j=2^{k}}$, і ${\displaystyle k}$ пробігає набір послідовних значень 1, 2, 3,…, а ${\displaystyle i}$ набуває значення з інтервалу ${\displaystyle [2^{k}+1;2^{k+1}]}$. Наприклад, ${\displaystyle k=3}$, ${\displaystyle j=2^{3}=8}$, а ${\displaystyle i\in [9;16]}$^[9].

Цю ідею запропонував Річард Брент^[ru] у 1980 році^[10] і вона дозволяє зменшити кількість виконуваних операцій приблизно на чверть (25%)^[11].

Ще одну варіацію ρ-методу Поларда розробив Флойд. За Флойдом, значення ${\displaystyle y}$ оновлюється на кожному кроці за формулою ${\displaystyle y=F^{2}(y)=F(F(y))}$, тому на кроці i будуть отримані значення ${\displaystyle x_{i}=F^{i}(x_{0})}$, ${\displaystyle y_{i}=x_{2i}=F^{2i}(x_{0})}$, і НСД на цьому кроці обчислюється для ${\displaystyle N}$ та ${\displaystyle y-x}$^[7].

Приклад факторизації числа

Нехай ${\displaystyle N=8051}$, ${\displaystyle F(x)=(x^{2}+1)\,\mathrm {mod} \,8051}$, ${\displaystyle x_{0}=y_{0}=2}$, ${\displaystyle y_{i+1}=F(F(y_{i}))}$.

i	xi	yi	НСД (|xi −yi|, 8051)
1	5	26	1
2	26	7474	1
3	677	871	97

Таким чином, 97 — нетривіальний дільник числа 8051. Використовуючи інші варіанти поліному ${\displaystyle F(x)}$, можна також отримати дільник 83.

Обґрунтування ρ-методу Полларда

Алгоритм ґрунтується на відомому парадоксі днів народження.

Теорема (Парадокс днів народження)

Нехай ${\displaystyle \lambda >0}$. Для випадкової вибірки з ${\displaystyle l+1}$ елементів, кожний з яких менший від ${\displaystyle q}$, де ${\displaystyle l={\sqrt {2\lambda q}}}$, ймовірність того, що два елементи виявляться однаковими ${\displaystyle p>1-\exp ^{-\lambda }}$.

Слід зазначити, що ймовірність ${\displaystyle p=0.5}$ в парадоксі днів народження досягається при ${\displaystyle \lambda \approx 0.69}$.

Нехай послідовність ${\displaystyle \{u_{n}\}}$ складається з різниць ${\displaystyle x_{i}-x_{j}}$, що перевіряються під час роботи алгоритму. Визначимо нову послідовність ${\displaystyle \{z_{n}\}}$, де ${\displaystyle z_{n}=u_{n}\,\mathrm {mod} \,q}$, ${\displaystyle q}$ — менший з дільників числа ${\displaystyle N}$.

Усі члени послідовності ${\displaystyle \{z_{n}\}}$ менші ${\displaystyle {\sqrt {N}}}$. Якщо розглядати її як випадкову послідовність цілих чисел, менших ${\displaystyle q}$, то, згідно з парадоксом днів народження, імовірність того, що серед ${\displaystyle l+1}$ її членів трапляться два однакових, перевищить ${\displaystyle 1/2}$ при ${\displaystyle \lambda \approx 0.69}$, тоді ${\displaystyle l}$ має бути не менше ${\displaystyle {\sqrt {2\lambda q}}\approx {\sqrt {1.4q}}\approx 1.18{\sqrt {q}}}$.

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

Складність алгоритму

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

Тому складність алгоритму оцінюється, як ${\displaystyle O(N^{1/4})}$^[12]. Однак у цій оцінці не враховуються накладні витрати з обчислення найбільшого спільного дільника. Отримана складність алгоритму, хоча і не є точною, проте достатньо добре узгоджується з практикою.

Виконується така теорема. Нехай ${\displaystyle N}$ — складене число. Тоді існує така стала ${\displaystyle C}$, що для будь-якого додатного числа ${\displaystyle \lambda }$ ймовірність події, що полягає в тому, що ρ-метод Поларда не знайде нетривіального дільника ${\displaystyle N}$ за час ${\displaystyle C{\sqrt {\lambda {\sqrt {N}}}}(\log N)^{2}}$, не перевершує величини ${\displaystyle e^{-\lambda }}$. Ця теорема випливає з парадоксу днів народження.

Особливості реалізації

Обсяг пам'яті, використовуваний алгоритмом, можна значно зменшити.

 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));
 }

у цьому варіанті обчислення потребує зберігати в пам'яті всього три змінні ${\displaystyle N}$, ${\displaystyle x}$, і ${\displaystyle y}$, що вигідно відрізняє метод в такій реалізації від інших методів факторизації чисел^[7].

Розпаралелювання алгоритму

Алгоритм Полларда дозволяє розпаралелювання з використанням будь-якого стандарту паралельних обчислень (наприклад, OpenMP та ін.).

Існує декілька варіантів розпаралелювання, але їх спільна ідея полягає в тому, що кожний процесор виконує послідовний алгоритм, причому початкове число ${\displaystyle x_{0}}$ та/або поліном ${\displaystyle F(x)}$ мають бути різними для кожного процесора. Очікується, що на якомусь процесорі початкові параметри (випадково) виявляться більш вдалими і дільник буде знайдено швидше, однак цей випадок недетермінований (імовірнісний). Прискорення від такої паралельної реалізації значно менше лінійного.

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

Річард Крандалл припустив, що можна досягти прискорення ${\displaystyle O(P/(\log {P})^{2})}$, однак на 2000-й рік це твердження не було перевірено^[13].

Див. також

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