Решето Эратосфена

Чтобы выпустить ключ RSA на 2048 бит, программе нужны два случайных простых числа по 1024 бита — в каждом 309 десятичных цифр. Готового списка таких чисел нет и быть не может, поэтому программа берёт случайное нечётное число и идёт от него вверх, пока не встретит простое. Простые среди нечётных чисел такого размера попадаются в среднем одно на 355, а проверка кандидата тестом Миллера — Рабина — это возведение в степень с показателем из 309 цифр по модулю 309-значного числа, больше тысячи умножений огромных чисел. Можно проверять всех подряд. А можно сначала просеять: вычеркнуть кандидатов, у которых есть маленький делитель, — их больше восьмидесяти процентов, — и тратить дорогой тест только на выживших. Ниже это происходит в вашем браузере, на настоящих числах.

Каждая клетка — нечётный кандидат: x, x + 2, x + 4 и так далее. Сначала решето вычёркивает каждую третью клетку — кратные трёх, потом каждую пятую, седьмую, и через несколько простых окно уже в дырах. Потом по выжившим идёт дорогой тест, и первое число, которое он не опровергнет, — простое. Тест вероятностный, но после 64 раундов шанс ошибиться меньше 4⁻⁶⁴. Подвигайте ползунок «решето до». Без решета тест достаётся каждому кандидату. Решето до тысячи оставляет около 16 % кандидатов, и дорогих проверок нужно примерно вшестеро меньше. Решето до миллиона оставляет 8 % — выигрыш растёт, но всё медленнее, а само решето заметно дорожает: это 78 тысяч делений огромного числа на маленькие.

Идея

Не проверять числа на делимость по одному, а идти в обратную сторону — от делителей к кратным. Каждое простое p само показывает, где стоят его кратные: через каждые p чисел. Шагать по массиву с шагом p и ставить отметки — ни одного деления, только сложения.

Классическое решето ищет все простые до n. Выпишем числа от 2 до n. Самое маленькое невычеркнутое число — простое: если бы у него был меньший делитель, оно было бы вычеркнуто как кратное этого делителя. Объявляем его простым, вычёркиваем его кратные и повторяем. Две детали делают решето быстрым:

• кратные p вычёркиваем начиная с p², а не с 2p: числа 2p, 3p, …, (p − 1)p имеют делитель меньше p и уже вычеркнуты;
• останавливаемся, когда p² превысит n: вычёркивать дальше нечего, и всё, что осталось, — простые.

Подробный школьный разбор с картинкой и доказательством — в главе «Простые числа» «Царицы наук». Здесь нас интересует другое: сколько работы делает решето, где оно делает её зря и как его ускоряют на практике.

Сколько раз вычёркивают одно число

Двойка вычёркивает 4, 6, 8, …, тройка — 9, 12, 15, … Число 12 попадёт под обе: его вычеркнут дважды. Число 30 — трижды: двойкой, тройкой и пятёркой. Решето не помнит, что число уже вычеркнуто, и тратит время на повторы. Пройдите решето по шагам и посмотрите на точки: столько раз вычеркнули число.

До ста решето делает 104 вычёркивания, а составных чисел там 74: тридцать вычёркиваний лишние. Переключитесь на «Линейное» — это другое решето, о нём ниже: оно вычёркивает каждое составное ровно один раз, его наименьшим простым делителем, и делает ровно 74 записи. Подпись под числом показывает, как его разложило линейное решето: 60 = 2 · 30.

Почему остаются ровно простые

Та же мысль даёт и остановку на √n: у составного числа до n всегда есть простой делитель не больше √n. Перебор делителей одного числа останавливается на корне по той же причине.

Код

Массив флагов «простое ли число», по байту на число. В Python вычёркивание кратных одного простого — одно присваивание срезу is_prime[p*p::p] = bytes(count): цикл крутится внутри интерпретатора, на C, и это в десятки раз быстрее питоновского for. В Go пишем обычный цикл: компилятор сделает из него то же самое. Программа печатает простые до ста и считает, сколько вычёркиваний решето делает на разных n.

         n   primes  crossings  per number  ln ln n
      1000      168       1411       1.411    1.933
   1000000    78498    2122048       2.122    2.626
  10000000   664579   22850051       2.285    2.780

Посмотрите на столбец «на число»: от тысячи до десяти миллионов n выросло в десять тысяч раз, а среднее число вычёркиваний — с 1,4 до 2,3. Решето почти линейно, и следующий раздел объясняет, откуда это «почти».

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

Простое p вычёркивает примерно n/p чисел. Всего вычёркиваний меньше, чем

Сумма обратных простых растёт как двойной логарифм — это теорема Мертенса 1874 года, константа 0,2615 тоже его. Двойной логарифм растёт так медленно, что его трудно отличить от константы: для n = 10⁷ он равен 2,78, для n = 10¹⁸ — 3,72. Отсюда O(n log log n) по времени. Память — O(n): байт на число, а если хранить биты и только нечётные числа — шестнадцать чисел на байт.

Перебор делителей для каждого числа до n гораздо дороже: простое k требует около √k делений, и всего выходит порядка n√n/ln n. Ниже замер: четыре способа на одних и тех же n, оси логарифмические, поэтому наклон прямой показывает степень. У перебора наклон около 1,4 — полторы степени, слегка придавленные логарифмом, — у решёт около 1.

Две вещи на графике неочевидны. Линейное решето, у которого меньше всего записей, не быстрее обычного: в Go они идут вровень, а в Python линейное отстаёт раз в сорок — его внутренний цикл написан на Python, а у обычного решета цикл спрятан в срезе. И на самом большом n простое решето в Go дорожает быстрее сегментного. Это не алгоритм, а память. Массив в 100 МБ не помещается ни в какой кэш процессора, и шаг с большим p каждый раз попадает в новую строку памяти. Сегментное решето делает те же вычёркивания по кускам, которые помещаются в кэш, и на 10⁸ выигрывает почти втрое.

А теперь то же в вашем браузере. Выберите n и запустите: способы пойдут по очереди, в таблице — память, число операций и время.

На миллиарде остаётся одно сегментное решето: байт на число — это гигабайт, а сегментному хватает окна в 32 КБ и трёх с небольшим тысяч простых до √10⁹ ≈ 31 623. Ответ 50 847 534 сходится с таблицами.

Линейное решето

Можно ли вычёркивать каждое составное ровно один раз? Да, если поручить его одному определённому простому — наименьшему делителю. Каждое составное m единственным образом записывается как m = p · i, где p — наименьший простой делитель m; тогда у i все простые делители не меньше p. Значит, достаточно для каждого i перебрать простые p, не большие наименьшего делителя i, и записать p как наименьший делитель числа p · i. Каждая пара (p, i) даёт своё число, и каждое составное получается ровно из одной пары: всего записей столько же, сколько составных. Это линейное решето, его время работы O(n) показали Дэвид Грис и Джаядев Мисра в 1978 году.

Программа подтверждает: до 400 000 составных 366 139, и записей ровно столько же. Но главная ценность линейного решета не в скорости — выше видно, что на практике оно её не даёт, — а в побочном продукте: таблице наименьших простых делителей spf. С ней любое число до n раскладывается на множители за O(log n) делений: делим на spf[m], пока не останется единица. Так программа нашла число делителей каждого числа до 400 000: больше всех их у 332 640, целых 192. Платить за это приходится памятью: четыре байта на число вместо одного бита.

Таблицу наименьших делителей можно получить и обычным решетом: при вычёркивании записывать в число p, если там ещё пусто. Повторные вычёркивания останутся, но цикл проще, и на практике разница невелика.

Сегментное решето: окно где угодно

Чтобы вычеркнуть составные в окне [L, R), не нужен массив от нуля: достаточно простых до √R и массива размером с само окно. Для каждого простого p находим первое кратное в окне — ⌈L/p⌉ · p, но не меньше p², чтобы не вычеркнуть само p, — и шагаем по окну с шагом p. Это сегментное решето. Разрезав отрезок от 0 до n на окна по 32 КБ, получаем то же решето с памятью O(√n), и каждое окно целиком живёт в кэше — отсюда выигрыш в замере выше.

Окно может лежать далеко от нуля. Программа ниже считает простые до десяти миллионов окнами по 32 768 чисел, а потом просеивает миллион чисел после 10¹²: для этого хватает 78 498 простых до миллиона, а весь массив от нуля занял бы терабайт.

[10^12, 10^12 + 10^6): 36249 primes, using 78498 primes up to 10^6
first 1000000000039, last 1000000999999
widest gap 236 after 1000000494263

Простых в этом миллионе 36 249 — примерно одно на 27,6 числа, как и обещает закон распределения простых: ln 10¹² ≈ 27,6. А самый широкий промежуток без простых — 236 чисел подряд.

Первая сцена статьи — тоже сегментное решето, только неполное. До √x для 309-значного x никакого решета не хватит: это число из 155 цифр. Поэтому окно просеивают лишь малыми простыми, до тысячи или до миллиона. Выжившие — не обязательно простые, это числа без маленьких делителей; остальное доделывает тест Миллера — Рабина. Доля выживших тоже предсказуема — по третьей теореме Мертенса:

где γ ≈ 0,5772 — постоянная Эйлера. Для B = 719 формула даёт 17,1 %, точное произведение — 17,0 %, а в окне из листинга ниже выжило 344 кандидата из 2000, 17,2 %. Знаменатель — логарифм, поэтому каждый следующий десятикратный рост B срезает всё меньше: с тысячи до миллиона доля падает только вдвое, с 16 % до 8 %, а делений огромного числа на малые простые становится в 470 раз больше.

Число здесь «случайное» понарошку — блоки SHA-256 от фразы, чтобы Python и Go нашли одно и то же. От старта до простого 457 нечётных кандидатов. Решето по 127 простым до 719 сократило дорогие проверки до 76, решето до ста тысяч — до 44.

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

• Ключи. OpenSSL перед тестом Миллера — Рабина отсеивает кандидатов первыми простыми: 127 нечётными до 719 для 1024-битных чисел и 383 до 2657 для 2048-битных. В функции BN_generate_prime_ex остатки случайного старта от деления на эти простые считаются один раз, а дальше кандидат сдвигается на 2, и отбрасываются те сдвиги, при которых остаток обнуляется, — то же решето на окне, только ленивое.
• primesieve. Открытая библиотека Кима Валиша для быстрой генерации простых: сегментное решето Эратосфена с колёсной факторизацией — кратные самых маленьких простых выброшены заранее — и с размером окна под кэши L1 и L2 конкретного процессора.
• Гипотеза Гольдбаха. Томаш Оливейра-и-Силва, Зигфрид Херцог и Сильвио Парди проверили, что каждое чётное число до 4·10¹⁸ — сумма двух простых. Простые для этого давало сегментное решето, дружественное к кэшу; ушло около 782 процессорных лет. Работа опубликована в 2014 году.
• Олимпиады. Таблица наименьших делителей до 10⁶–10⁷ — стандартная заготовка: разложение чисел, функция Эйлера, число делителей для многих запросов сразу.

Где решето не нужно и где на нём спотыкаются

• Одно число. Чтобы проверить на простоту одно число, решето не нужно. Для 12-значного хватит перебора делителей до √n — миллиона делений, а для 20-значного и больших — теста Миллера — Рабина, нескольких возведений в степень. Решето окупается, когда нужны все простые подряд или проверка многих близких чисел.
• Память. Байт на число до 10⁹ — гигабайт. Биты и одни нечётные числа дают 62,5 МБ, сегменты — десятки килобайт.
• Переполнение. p * p при p ≥ 46 341 не помещается в 32-битное целое. В Go int 64-битный, в Python целые не переполняются, но в C и Java на этом спотыкаются постоянно.
• Окно с единицей. Если окно начинается с 0 или 1, их нужно вычеркнуть отдельно, а простые из самого окна не вычёркивать — отсюда «не меньше p²» в сегментном решете.

История

Сочинения Эратосфена, жившего в III веке до нашей эры, почти не сохранились. О решете мы знаем из «Введения в арифметику» Никомаха, написанного около 100 года нашей эры: Никомах приписывает способ Эратосфену и описывает его на одних нечётных числах — оптимизация «только нечётные» старше самого слова «алгоритм». Линейное решето с доказательством линейного времени опубликовали Дэвид Грис и Джаядев Мисра в 1978 году. Сегментное решето известно с 1970-х; одна из первых работ — статья Картера Бейса и Ричарда Хадсона 1977 года, где они просеивали простые до 10¹². Пол Притчард в 1981 году построил колёсное решето, которому нужно даже меньше O(n) операций, но на практике хорошо написанное решето Эратосфена с колесом не уступает ему.

Задачи

Первая — сегментное решето в чистом виде, на окнах далеко от нуля. Во второй решето нужно не для простых: сначала придётся понять, что считать.