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

В 1903 году Фрэнк Нельсон Коул вышел к доске на заседании математического общества, молча посчитал 2⁶⁷ − 1, потом перемножил 193 707 721 и 761 838 257 287 — результаты совпали. Позже он признавался, что на поиск делителей ушли «три года воскресений». Перебору делителей пришлось бы дойти до 193 707 721 — это почти сто миллионов делений на одни нечётные числа. ρ-алгоритм Полларда находит тот же делитель за четырнадцать тысяч шагов. Он вообще не ищет делитель: он ждёт, пока псевдослучайная последовательность сама его выдаст, повторившись по модулю неизвестного нам простого числа. Ниже он раскладывает настоящие числа — от числа Ферма, которое Эйлер разложил в 1732 году, до 78-значного числа Ферма F₈.

Шаг здесь — одно вычисление x ↦ x² + c по модулю n, НОД — сколько раз считался наибольший общий делитель. Пока алгоритм работает, строка под ответом показывает, докуда он «дотянулся»: после s шагов делитель меньше (s/1,25)², скорее всего, уже был бы найден. После миллиона шагов это 6·10¹¹, а перебор за миллион делений проверяет делители только до двух миллионов. На 30-значном числе с 13-значным делителем ρ тратит около двух миллионов шагов, на 40-значном с 16-значным — около тридцати миллионов. Время перебора в тех же случаях — часы и месяцы. Число Ферма F₈ с 16-значным делителем браузер раскладывает за несколько секунд, а с приёмом, о котором ниже, — в несколько раз быстрее.

Идея

Возьмём последовательность x_i+1 = x_i² + c mod n и посмотрим на неё по модулю неизвестного простого делителя p. Там всего p возможных значений, и по парадоксу дней рождения повтор x_i ≡ x_j (mod p) появится примерно через √p шагов. Заметить его можно, не зная p: при повторе p делит x_i − x_j, и НОД(x_i − x_j, n) больше единицы.

Повтор по модулю p наступает гораздо раньше, чем по модулю n: значений по модулю n в n/p раз больше. Поэтому НОД обычно равен p или его кратному, но не самому n. Осталось не перебирать все пары (i, j) — их квадрат от числа шагов, — а заметить повтор, храня всего пару чисел. Для этого есть приёмы поиска цикла.

Буква ρ

По модулю простого p последовательность детерминирована: следующее значение зависит только от текущего. Значит, после первого повтора она ходит по кругу. Нарисуйте её путь — сначала хвост из μ значений, потом цикл длины λ, — и получится греческая буква ρ; отсюда название. Ниже настоящие последовательности по модулю 101, 1009 и 10 007. Меняйте c и стартовое значение и смотрите, как бегают черепаха и заяц.

Флойд. Черепаха делает шаг, заяц — два, и на каждом ходу их значения сравнивают. Пока оба на хвосте, заяц убегает вперёд, но, попав на цикл, начинает нагонять черепаху на одну точку за ход и встречает её, не пройдя и круга. В алгоритме Полларда «сравнить» значит посчитать НОД(x_i − x_2i, n): равенства по модулю p мы не видим, зато видим общий делитель.

Брент. Заяц бежит по одному шагу, а черепаха стоит и перепрыгивает к зайцу, когда тот пробежал 1, 2, 4, 8, … шагов. Отрезок между прыжками рано или поздно станет длиннее цикла, и заяц придёт к черепахе. Каждое значение вычисляется один раз, а не три, поэтому вычислений f меньше.

Гистограмма под буквой — длина ρ для каждого из p возможных стартов при выбранном c. Линия — то, что предсказывает случайное отображение на p точках, со средним √(πp/2) ≈ 1,25√p. Отображение x² + c вовсе не случайное, но ведёт себя похоже, и на этом держится весь алгоритм.

Почему √p, а не √n

Алгоритм держится на трёх утверждениях и одной гипотезе. Первые два утверждения — о самой последовательности, третье — о настоящих случайных числах, а гипотеза — в том, что x² + c ведёт себя как случайное.

Отсюда буква ρ: по модулю p всего p значений, следующее определяется текущим, поэтому не позже чем через p шагов значение повторится, и дальше последовательность пойдёт по кругу. Хвост μ и цикл λ вместе не длиннее p.

Остаётся оценить μ + λ — номер первого повтора. Если бы значения выбирались случайно и независимо из p возможных, это была бы задача о днях рождения.

При k = 1,2√p правая часть меньше половины, при k = 3√p — около одного процента: повтор почти наверняка случится за несколько √p шагов. Среднее можно посчитать точно, оно около √(πp/2) ≈ 1,25√p. А по модулю n значений в n/p раз больше, повтор там ждать в √(n/p) раз дольше — поэтому НОД почти всегда даёт p, а не n.

        p  tries  rho/sqrt p  Floyd/sqrt p  Brent/sqrt p
    10007   1000       1.224         3.042         1.927
   100003    600       1.299         3.236         2.049
  1000003    300       1.249         3.098         2.028
 10000019    100       1.318         3.187         2.105

Флойду в среднем нужно около 3,1√p вычислений f — черепаха делает одно, заяц два на каждом ходу, — Бренту около 2,0√p, на треть меньше. Сам Брент в статье 1980 года оценил выигрыш в поиске цикла в 36 %, а в разложении на множители — в 24 %.

Код

Разложение целиком: маленькие делители и простота кусков — на совести теста Миллера — Рабина, а составной кусок делит ρ. В варианте Брента два улучшения. Черепаха ждёт на степенях двойки, как в сцене. И НОД считается не на каждом шаге: разности x − y перемножаются по модулю n, и раз в 128 шагов берётся один НОД произведения. Если p делит хоть одну разность, оно делит и произведение. Редкий риск — что в одну пачку попадут сразу все делители и НОД окажется равен n, — снимается повтором последней пачки по одному шагу.

2^64 + 1 = 274177 * 67280421310721
    Floyd: 2424 steps, 808 gcds; Brent: 1918 steps, 13 gcds; trial division: 4101232 divisions
2^67 - 1 = 193707721 * 761838257287
    Floyd: 16584 steps, 5528 gcds; Brent: 13822 steps, 50 gcds; trial division: 96853860 divisions
p * q = 999999937 * 1000000007
    Floyd: 36777 steps, 12259 gcds; Brent: 28670 steps, 102 gcds; trial division: 499999968 divisions

Шагов у Брента примерно на пятую часть меньше, а НОД — в сто с лишним раз: на 2⁶⁷ − 1 их 50 вместо 5528. НОД двух длинных чисел — это целый алгоритм Евклида, заметно дороже одного умножения, так что пачки экономят больше, чем сам поиск цикла. Число 2⁶⁴ + 1 = 274 177 · 67 280 421 310 721 разложил Фортюне Ландри в 1880 году, без всяких вычислительных машин.

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

Шагов — порядка √p, где p — наименьший простой делитель, и на шаг уходит умножение и взятие остатка по модулю n: O(log² n) столбиком. Всего O(√p · log² n), памяти O(1). Худший случай — два множителя одной длины, p ≈ √n, тогда шагов порядка n^1/4. Для сравнения, перебор делает до p делений, если нужно найти p, и до √n, если нужно доказать простоту.

Замер показывает главное свойство. Числа n из 64, 128 и 256 бит, у каждого наименьший делитель p заданной длины, время до первого делителя:

Все три прямые идут с наклоном около ½: в шестнадцать миллионов раз больший p обходится в несколько тысяч раз дороже. Длина n на наклон не влияет — она только сдвигает прямую вверх, потому что умножать длинные числа дороже: при p около 2³⁰ 256-битное число раскладывается в 1,7 раза дольше 64-битного на Python и почти в 5 раз — на Go. Ещё одна неожиданность: Go с math/big здесь не быстрее Python. Числа короткие, и время уходит не на арифметику, а на вызовы и выделения памяти; на uint64 с bits.Mul64 Go обогнал бы Python в разы.

Приём для чисел Ферма

Иногда о делителях что-то известно заранее. Эйлер доказал, что простые делители числа Ферма F_m = 2^2m + 1 имеют вид k·2^m+1 + 1, а Эдуард Люка в 1878 году усилил это до k·2^m+2 + 1. У делителя F₈ это значит, что p − 1 делится на 2¹⁰. Возьмём вместо x² + c отображение x¹⁰²⁴ + c. По модулю такого p степень x¹⁰²⁴ принимает всего около p/1024 различных значений, последовательность бегает по множеству в тысячу раз меньшему, и цикл в √1024 = 32 раза короче. Шаг стоит десять возведений в квадрат вместо одного, но выигрыш больше. В первой сцене выключите этот приём флажком и сравните: без него F₈ требует около тридцати миллионов шагов, с ним — меньше миллиона.

Ровно на этом числе Ричард Брент и Джон Поллард в 1980 году показали силу метода. F₈ было известно составным, но ни одного делителя не знали. Модификация ρ-алгоритма нашла 1 238 926 361 552 897 за два часа работы машины Univac 1100/42 — около 23 миллионов умножений по модулю F₈. Простоту оставшегося 62-значного множителя первым доказал Хью Уильямс.

Что делает разложение трудным

ρ-алгоритм быстр ровно настолько, насколько мал наименьший делитель. У случайного большого числа почти всегда есть маленькие делители, и ρ вместе с тестом Миллера — Рабина раскладывает его до последнего, обычно простого, куска. Трудны числа, у которых маленьких делителей нет вовсе: произведение двух простых одной длины. Для 30-значного произведения двух 15-значных простых ρ нужно порядка 6·10⁷ шагов — секунды в браузере; для 60-значного — порядка 10¹⁵, для 600-значного — 10¹⁵⁰. Именно такие числа и берут в RSA, и подробнее о стойкости — в главе «Секрет на виду у всех» базового курса и в разделе «Взлом» «Царицы наук».

Другие методы делят ту же задачу иначе. Метод p − 1, тоже Полларда (1974), находит p мгновенно, если у p − 1 все простые делители маленькие, и бессилен в остальных случаях. Метод эллиптических кривых Хендрика Ленстры (1985) — тоже «ищет маленький делитель», но его время растёт медленнее любой степени p; им находят делители из десятков цифр. Квадратичное решето и решето числового поля работают иначе: их время зависит от размера n, а не p, и для чисел из сотен цифр без маленьких делителей они вне конкуренции. Решето числового поля Поллард тоже предложил первым, в 1988 году; в 1990-м им разложили девятое число Ферма из 155 цифр. На практике методы выстраиваются в конвейер: пробное деление, ρ, эллиптические кривые, решето.

Тот же приём для логарифмов

В 1978 году Поллард применил ту же идею к дискретному логарифму: найти x, при котором g^x = h в группе из N элементов. Псевдослучайное блуждание по группе снова зацикливается за √N шагов, а из совпадения двух точек пути извлекается x. Для эллиптических кривых лучших общих методов не известно, поэтому стойкость кривой над полем из 2²⁵⁶ элементов и оценивают в 2¹²⁸ шагов — об этом в главе «Кривые в браузере».

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

• Утилита factor из GNU coreutils раскладывает числа так же, как программа выше: пробное деление по таблице малых простых, проверка простоты Миллером — Рабином и тестом Люка, затем ρ Полларда — Брента или метод SQUFOF для оставшегося составного куска.
• PARI/GP: функция factorint сочетает ρ-алгоритм, SQUFOF Шенкса, метод эллиптических кривых и квадратичное решето.
• Олимпиады: разложить на множители числа до 10¹⁸ — стандартная задача на связку Миллера — Рабина и ρ с умножением через 128-битное произведение.
• Криптоанализ эллиптических кривых: ρ для логарифмов, распараллеленный на много машин, — основной способ рекордных вычислений дискретного логарифма на кривых.

История

Джон Поллард опубликовал метод в 1975 году в короткой статье «Метод Монте-Карло для разложения на множители». Цикл в ней ловился способом, который Дональд Кнут в 1969 году без ссылки приписал Роберту Флойду; в работах самого Флойда этого алгоритма нет, и, возможно, он фольклорный. В 1980 году Ричард Брент заменил его своим поиском цикла с черепахой на степенях двойки и в том же году вместе с Поллардом разложил F₈. Делители чисел Ферма и Мерсенна, раньше стоившие «лет воскресений», стали задачей на часы машинного времени, а сейчас — на доли секунды в браузере.

Задачи

Восстановите начало таблицы Каннингема: разложите 2^k ± 1 для всех k до ста. Перебора до миллиона хватит для большинства строк, но не для всех.