Тест Миллера — Рабина

Число 18 446 744 073 709 551 557 простое? Перебор делителей до корня — это два миллиарда делений, и Python по нашему замеру потратит на них около минуты. Тест Миллера — Рабина отвечает за десятую долю миллисекунды и не находит при этом ни одного делителя. Он вызывает свидетелей — числа-основания — и сверяет их показания. Простое число запутать их не может, а составное почти наверняка проговорится. Ниже суд идёт в вашем браузере. Выберите подсудимого — среди них есть рекордсмены, которые обманывают тринадцать свидетелей подряд, — или введите своё число.

В каждой строке одно основание a даёт показания двумя способами. Колонка «Ферма» — старый тест: остаток aⁿ⁻¹ по модулю n. Для простого n это всегда 1 — так говорит малая теорема Ферма. Колонка «цепочка» — та же степень, но посчитанная с остановками: сначала a^d, потом её квадрат, квадрат квадрата и так далее до aⁿ⁻¹. Посмотрите на 561: тест Ферма на большинстве оснований молчит, а цепочка выдаёт число с первого же основания. На ψ₁₃ лгут тринадцать оснований подряд, и только 43 говорит правду. А на числе Кармайкла из 61 знака тест Ферма не заметит подвоха ни при каком основании, взаимно простом с n, — зато цепочка не только ловит его, но и выдаёт делитель.

Идея

По простому модулю у единицы ровно два квадратных корня: 1 и −1. Запишем n − 1 = 2^s·d с нечётным d и будем возводить a^d в квадрат, пока не дойдём до aⁿ⁻¹. Если n простое, в конце цепочки стоит 1, и появиться она может только сразу или сразу после −1. Если в цепочке единица возникла из чего-то другого — или до единицы цепочка вообще не дошла, — n составное, и основание a это доказало.

Тест Ферма смотрит только на последнее звено цепочки. Его обманывает любое составное n, для которого aⁿ⁻¹ ≡ 1, и такие числа есть для любого основания. Цепочка требует больше: не просто прийти к единице, а прийти к ней правильной дорогой. Обмануть её тоже можно, но гораздо труднее.

Назовём основание a свидетелем, если цепочка его выдаёт: a^d ≢ 1 и ни одно из чисел a^d, a^2d, …, a^2s−1d не равно −1 по модулю n. Иначе говорят, что n — сильно вероятно простое по основанию a. Составное число, которое прошло проверку, называют сильным псевдопростым по основанию a, а само основание — сильным лжецом.

Допрос по шагам

Возьмём простое 97. Тогда 96 = 2⁵ · 3: s = 5, d = 3. Основание 5: 5³ = 125 ≡ 28, дальше квадраты — 28² ≡ 8, 8² = 64, 64² ≡ 22, 22² ≡ 96 ≡ −1, и последний квадрат (−1)² = 1. Цепочка 28, 8, 64, 22, −1, 1 пришла к единице через −1, как и положено простому.

Теперь 2047 = 23 · 89. Здесь 2046 = 2 · 1023, звеньев два. Основание 3: 3¹⁰²³ ≡ 1565 — ни 1, ни −1, а квадратов больше делать нельзя: следующее звено уже aⁿ⁻¹. Свидетель. Основание 2: 2¹⁰²³ ≡ 1, потому что 2¹¹ = 2048 ≡ 1, а 1023 делится на 11. Двойка лжёт. Число 2047 — наименьшее сильное псевдопростое по основанию 2, и всего у него 242 лжеца среди 2046 оснований, меньше 12 %.

На каждое основание уходит одно быстрое возведение в степень — около log₂ n умножений по модулю — и не больше s возведений в квадрат. Основания перебирают по очереди, и первый же свидетель заканчивает суд: «составное» — это доказательство, а не догадка. Если ни одно основание не возразило, ответ «вероятно простое», и насколько «вероятно», решает следующий раздел.

Почему простое не проваливается, а составное проваливается часто

Обратное утверждение и делает тест полезным: у составного числа лжецов мало. Докажем, что их не больше половины. Нам понадобится одно наблюдение про остатки, обратимые по модулю n, — те, что взаимно просты с n; их φ(n) штук.

Отсюда главное свойство теста. Случайное основание ловит составное n с вероятностью больше ½, и k независимых оснований ошибаются все сразу с вероятностью меньше 2^−k. Ошибка односторонняя: «составное» всегда верно, сомнение остаётся только в ответе «простое». Это алгоритм Монте-Карло, но с ошибкой, которую можно сделать меньше вероятности сбоя железа.

На самом деле граница вдвое строже. В 1980 году Майкл Рабин и независимо Луи Монье доказали, что у нечётного составного n > 9 лжецов не больше φ(n)/4, и ошибка k раундов меньше 4^−k. Монье нашёл и точную формулу для их числа. Пусть у n ровно r различных простых делителей, p − 1 = 2^s_p·d_p с нечётным d_p, и ν — наименьшее из s_p. Тогда сильных лжецов

а лжецов Ферма F(n) = ∏ gcd(n − 1, p − 1). Для 2047 = 23 · 89 формула даёт (1 + 3/3) · 11 · 11 = 242 — столько мы и насчитали. Вывод длиннее нашего доказательства: он опирается на то, что обратимые остатки по модулю степени нечётного простого получаются степенями одного остатка, — его можно найти в статье Монье 1980 года и в книге Крэндалла и Померанса «Простые числа». Мы его не повторяем, а проверяем: листинг ниже сравнивает формулу с прямым перебором всех оснований для всех нечётных составных до тысячи.

Атлас лжецов

Четверть — это худший случай. Как лгут настоящие числа? Сцена ниже раскладывает на множители каждое нечётное составное число до миллиона — их 421 502 — и по формуле Монье считает, какая доля оснований ему лжёт. Каждая точка — одно число, высота — доля лжецов, шкала логарифмическая. Нажмите на точку, чтобы узнать, что это за число, и отправьте его на суд в первую сцену.

Нижняя кромка облака — плотная кривая 2/(n − 1): числа, которым лгут только 1 и n − 1. Эти два основания проходят проверку у любого нечётного n, и брать их бессмысленно. Таких чисел до миллиона 261 279, 62 %. Каждая следующая кривая выше — свой небольшой набор лжецов: 6, 10, 14, 18 оснований. Средняя доля лжецов по всем 421 502 числам — 0,0235 %, в тысячу с лишним раз меньше четверти.

До самой границы ¼ поднимаются 20 чисел. Девятнадцать из них устроены одинаково — p·(2p − 1), где p и 2p − 1 простые, а p при делении на 4 даёт остаток 3: 15 = 3 · 5, 91 = 7 · 13, 703 = 19 · 37 и дальше до 954 271 = 691 · 1381. Двадцатое — 8911 = 7 · 19 · 67, число Кармайкла из трёх простых вида 4k + 3. Из того же семейства p(2p − 1) и рекордсмены первой сцены: ψ₁₂ = 399 165 290 221 · 798 330 580 441 и ψ₁₃ = 1 287 836 182 261 · 2 575 672 364 521. Это не совпадение: чем больше у числа лжецов, тем больше шансов, что солгут все основания подряд.

Переключите сцену на тест Ферма. Облако поднимется, а над линией ½ окажутся одни числа Кармайкла. Это та же лемма о половине: если у n есть хоть один свидетель Ферма, взаимно простой с n, то лжецы Ферма — набор, замкнутый относительно умножения, — занимают не больше половины. Исключение — числа, у которых таких свидетелей нет вовсе.

Разница между худшим случаем и средним объясняет, почему стандарты по-разному считают раунды. Для числа, которое подбирал противник, честная оценка — 4^−k: OpenSSL делает не меньше 64 раундов для чисел до 2048 бит, это 2⁻¹²⁸, и 128 раундов для чисел длиннее. Для случайного кандидата при генерации ключа ошибка несравнимо меньше, и документация OpenSSL оценивает шанс получить составное вместо 2048-битного простого как 2⁻⁶⁴⁸.

Числа Кармайкла: почему одного Ферма мало

Составное n, для которого aⁿ⁻¹ ≡ 1 (mod n) при всех a, взаимно простых с n, называют числом Кармайкла. Тест Ферма ловит его, только если случайно выбрать основание с общим делителем, — а у числа из трёх двадцатизначных простых это почти невероятно. Распознать такие числа по разложению помогает признак Корсельта.

Обратное тоже верно — у числа Кармайкла нет квадратных делителей и p − 1 делит n − 1, — но для доказательства нужно существование первообразного корня, и мы его опускаем. Альвин Корсельт сформулировал признак в 1899 году, не зная ни одного примера. Первые семь чисел — 561, 1105, 1729, 2465, 2821, 6601 и 8911 — нашёл ещё в 1885 году Вацлав Шимерка, но его работу не заметили, и имя числам дал Роберт Кармайкл, описавший 561 в 1910 году. Из признака сразу видно, что простых делителей у числа Кармайкла хотя бы три: если n = pq и p < q, то n − 1 = p(q − 1) + (p − 1), и q − 1 должно делить меньшее его число p − 1.

Признак позволяет строить числа Кармайкла любого размера. Джек Черник в 1939 году заметил: если 6k + 1, 12k + 1 и 18k + 1 простые, их произведение — число Кармайкла. Действительно, произведение равно 1296k³ + 396k² + 36k + 1, и n − 1 = 36k(36k² + 11k + 1) делится на 6k, 12k и 18k. При k = 1 выходит 7 · 13 · 19 = 1729. Кнопка «новый Кармайкл» в первой сцене ищет такое k из двадцати цифр — обычно за несколько тысяч или десятков тысяч попыток — и получает число из шестидесяти с лишним знаков, на котором тест Ферма бессилен.

Таблица ниже — сколько чисел каждого сорта меньше 10^k. Мы пересчитали их сами до миллиарда — счёт совпал с опубликованными таблицами (Померанс, Селфридж и Вагстафф, 1980; Пинч для чисел Кармайкла).

Одно основание 2 на миллиарде ошибается 1282 раза на пятьдесят миллионов простых. Из 646 чисел Кармайкла до миллиарда сильную проверку по основанию 2 проходят только 43. Чисел Кармайкла бесконечно много — это доказали Алфорд, Гранвилл и Померанс в 1994 году, — но для сильного теста они не особенные: лжецов у них не больше четверти, как у всех.

Сколько оснований нужно на самом деле

Для чисел ограниченной длины случайность можно убрать совсем. Обозначим ψ_k наименьшее составное число, которое проходит сильный тест по всем первым k простым основаниям. Если n < ψ_k и прошло их все — оно простое, без всяких вероятностей. Значения ψ_k искали перебором с хитрыми отсечениями несколько десятилетий:

Число ψ₁₂ ≈ 3,2·10²³ больше 2⁶⁴ ≈ 1,8·10¹⁹, поэтому первых двенадцати простых, от 2 до 37, хватает для любого 64-битного числа. Короче набор Джима Синклера 2011 года: семь оснований 2, 325, 9375, 28178, 450775, 9780504, 1795265022 проверяют все n < 2⁶⁴. Такие наборы проверяют по списку всех псевдопростых по основанию 2 меньше 2⁶⁴, который вычислил Ян Фейтсма. Михал Форишек и Якуб Янчина в 2015 году пошли дальше: для 32-битных чисел хватает одного раунда, если основание выбирать из заранее составленной таблицы по хешу самого числа, а для 64-битных — двух или трёх.

Для чисел любой длины детерминированный вариант есть только в предположении. Гэри Миллер в 1976 году показал, что если верна расширенная гипотеза Римана, то достаточно проверить все основания до O(ln² n); Эрик Бах в 1990 году уточнил границу до 2 ln² n. Гипотеза не доказана, и на практике полагаются на случайные основания.

Код

Первая программа — суд из сцены в коде. Функция strong строит цепочку целиком, чтобы её можно было показать; в рабочем тесте цикл прерывают, как только встретили −1. Числа длинные, поэтому в Go — math/big, в Python — обычный int и встроенная pow(a, d, n), то есть быстрое возведение в степень по модулю.

number       digits  fermat  strong  bases 2..41 say            with 43 too
2^64 - 59        20   13/13   13/13  probably prime             probably prime
561               3   10/13    0/13  composite (witness 2)      composite (witness 2)
3215031751       10   13/13    6/13  composite (witness 11)     composite (witness 11)
psi12            24   12/13   12/13  composite (witness 41)     composite (witness 41)
psi13            25   13/13   13/13  probably prime             composite (witness 43)
chernick         61   13/13    2/13  composite (witness 2)      composite (witness 2)

Тест Ферма обманут числом Кармайкла из 61 знака на всех тринадцати основаниях, сильный тест — только на двух. А дальше та самая улика: цепочка для основания 2 — x x 1 1 1, единица появилась из числа, которое не равно ±1. Это число x — квадратный корень из единицы, и n делит (x − 1)(x + 1), не деля ни одну из скобок. Значит, НОД(x − 1, n) и НОД(x + 1, n) — нетривиальные делители. Программа получает 12k + 1 = 120 000 000 000 000 138 301 и произведение двух других простых. Тест, который не искал делителей, нашёл делитель даром.

Вторая программа сверяет формулу Монье с перебором и строит маленький атлас — до 300 000, чтобы уложиться в доли секунды. Доля чисел, которым лгут только 1 и n − 1, там 63,1 %, а в сцене до миллиона — 62,0 %: она падает, но очень медленно. Эрик Бах и Эндрю Шаллю в 2013 году показали, что доля таких чисел стремится к нулю — со скоростью, обратной тройному логарифму.

Третья — рабочая версия для 64-битных чисел с основаниями Синклера. В Python здесь нечего бояться, а в Go произведение двух 64-битных остатков не помещается в машинное слово. Функция bits.Mul64 возвращает все 128 бит произведения, bits.Div64 делит их на n. Перед возведением в степень стоит пробное деление на двенадцать маленьких простых: у 85 % всех чисел есть делитель не больше 37, и их отсеивают за малую долю того, что стоит один раунд. Программа сверяется с решетом до ста тысяч, ловит ψ₉, перед которым пасуют основания 2 … 23, находит восемь самых больших простых, помещающихся в 64 бита, и считает простые в окне после 10¹⁸: их 2398 при ожидаемых по закону распределения простых 2413.

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

Один раунд — одно возведение в степень: ⌊log₂ n⌋ возведений в квадрат и не больше стольких же домножений по модулю, плюс не больше s квадратов цепочки. Если числа длиной L = log n умножать столбиком, умножение стоит O(L²), и раунд — O(L³). Значит, k раундов — O(k log³ n). С умножением Карацубы показатель опускается до 2,585, с быстрым преобразованием Фурье — почти до двух. Перебор делителей до √n стоит √n = 10^L/2 делений: каждые две новые цифры — вдесятеро дольше. Замер: простое число из d цифр, перебор против двадцати раундов теста.

Перебор идёт по прямой вверх: на Python каждые две цифры умножают время в десять раз, и 16-значное простое проверяется уже почти шесть секунд; 20-значное заняло бы около десяти минут, 30-значное — почти два года. Миллер — Рабин для 20-значного числа тратит меньше миллисекунды, для числа из 512 цифр — полсекунды на Python и 166 миллисекунд в Go. На длинных числах его линия идёт с наклоном между двумя и тремя: удвоили длину — время выросло в 5–8 раз в Python, а в Go с большим разбросом, от 3,5 до 11 раз. Теория обещает наклон 3 — двойка от умножения столбиком и единица от числа умножений; на таких длинах заметную долю времени съедают накладные расходы, которые растут медленнее.

Где тест ломается

Фиксированные основания против того, кто подбирает число. Проверка «по основаниям 2…41» надёжна для чисел меньше ψ₁₃ и совершенно ненадёжна для чисел, которые подобраны против неё: ψ₁₃ — уже такое число. Строить их можно и больше. Франсуа Арно в 1995 году опубликовал число Кармайкла из 397 цифр, которое проходит сильный тест по всем простым основаниям меньше 307. В 2018 году Мартин Альбрехт, Джейк Массимо, Кеннет Патерсон и Юрай Соморовски проверили криптографические библиотеки и построили составные числа, которые библиотеки принимали за простые: для OpenSSL — 2048-битные числа, которые в настройках по умолчанию проходили проверку с вероятностью 1/16 вместо обещанных 2⁻⁸⁰; для GNU GMP — 1024-битные, которые всегда проходили с минимальным рекомендованным числом раундов; в Cryptlib, LibTomCrypt, JavaScript Big Number и WolfSSL — числа, которые проходили всегда. Опасность не теоретическая: противник подсовывает «простое» в параметрах Диффи — Хеллмана, и дискретный логарифм становится лёгким. В 2019 году OpenSSL добавил функцию проверки, у которой нельзя попросить слишком мало раундов: не меньше 64 для чисел до 2048 бит и 128 для более длинных; изменение ссылается на исследования Массимо и Патерсона.

Предсказуемые «случайные» основания. Оценка 4^−k верна, только если основания выбраны случайно и независимо от n. Основания из генератора с известным зерном — те же фиксированные основания, просто их сложнее увидеть.

Защита. Авторы той работы советуют тест Бейли — PSW: один сильный раунд по основанию 2 плюс сильный тест Люка, проверка совсем другой природы. Числа, которое обмануло бы обе части сразу, не нашли до сих пор, а для всех чисел меньше 2⁶⁴ его отсутствие проверено. В документации Go прямо сказано, что ProbablyPrime точен на 100 % для чисел меньше 2⁶⁴, но не годится для чисел, которые мог подобрать противник.

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

• Ключи. RSA нужны два случайных простых по 1024 бита и больше. Программа берёт случайное нечётное число, отсеивает его малыми простыми и проверяет тестом Миллера — Рабина, пока не встретит простое; как устроено отсеивание, рассказано в статье о решете Эратосфена. OpenSSL проверяет кандидатов не меньше чем 64 раундами.
• Библиотеки. В Go big.Int.ProbablyPrime(n) делает n раундов с псевдослучайными основаниями и тест Бейли — PSW. В GMP mpz_probab_prime_p(n, reps) — пробные деления, Бейли — PSW и ещё reps − 24 раундов. В стандартной библиотеке Python проверки простоты нет, но тест пишется в десять строк на встроенной pow(a, d, n).
• Разложение на множители. Вместе с ρ-алгоритмом Полларда тест раскладывает числа: кусок, который он назвал простым, дальше не делят.

История

Малая теорема Ферма известна с 1640 года, а в 1819 году Пьер Саррюс показал, что просто обратить её нельзя: 2³⁴¹ ≡ 2 (mod 341), хотя 341 = 11 · 31. Тест с цепочкой квадратов предложил Гэри Миллер в 1976 году — детерминированный, но доказанный только при условии расширенной гипотезы Римана. В 1980 году Майкл Рабин превратил его в вероятностный тест, которому никакие гипотезы не нужны, и доказал оценку в четверть; тогда же и независимо её получил Луи Монье. В том же году Померанс, Селфридж и Вагстафф пересчитали все псевдопростые до 25·10⁹ и нашли первые ψ_k.

В 2002 году Маниндра Агравал, Нирадж Каял и Нитин Саксена построили детерминированный тест простоты с полиномиальным временем без всяких гипотез, и вопрос «можно ли проверить простоту быстро» был закрыт в теории. На практике по-прежнему бросают кости: тест Миллера — Рабина несравнимо быстрее, а его ошибку можно сделать сколь угодно малой.

Задачи

В первой задаче нужно найти все числа Кармайкла в окне далеко от нуля — и придумать, как не раскладывать на множители каждое число окна. Во второй — сосчитать лжецов числа по его разложению, то есть вывести формулу Монье самостоятельно, хотя бы для одного случая.