Решето Аткина

Решето Аткина — алгоритм нахождения всех простых чисел до заданного целого числа N. Алгоритм был создан и . Заявленная авторами асимптотическая скорость работы алгоритма соответствует скорости лучших ранее известных алгоритмов просеивания, но в сравнении с ними решето Аткина требует меньше памяти.

Описание


 Основная идея алгоритма состоит в использовании неприводимых квадратичных форм (представление чисел в виде ax²+by²). Предыдущие алгоритмы в основном представляли собой различные модификации решета Эратосфена, где использовалось представление чисел в виде редуцированных форм (как правило — в виде произведения xy).
 В упрощённом виде алгоритм может быть представлен следующим образом:

  • Все числа, равные (по модулю 60) 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56 или 58, делятся на два и заведомо не простые. Все числа, равные (по модулю 60) 3, 9, 15, 21, 27, 33, 39, 45, 51 или 57, делятся на три и тоже не являются простыми. Все числа, равные (по модулю 60) 5, 25, 35 или 55, делятся на пять и также не простые. Все эти остатки (по модулю 60) игнорируются.
    • Все числа, равные (по модулю 60) 1, 13, 17, 29, 37, 41, 49 или 53, имеют остаток от деления на 4, равный 1. Эти числа являются простыми тогда и только тогда, когда количество решений уравнения 4x² + y² = n нечётно и само число не кратно никакому квадрату простого числа (:en:square-free integer).
    • Числа, равные (по модулю 60) 7, 19, 31, или 43, имеют остаток от деления на 6, равный 1. Эти числа являются простыми тогда и только тогда, когда количество решений уравнения 3x² + y² = n нечётно и само число не кратно никакому квадрату простого.
    • Числа, равные (по модулю 60) 11, 23, 47, или 59, имеют остаток от деления на 12, равный 11. Эти числа являются простыми тогда и только тогда, когда количество решений уравнения 3x² − y² = n (для x \textgreater y) нечётно и само число n не кратно никакому квадрату простого.
    • Отдельный шаг алгоритма вычеркивает числа, кратные квадратам простых чисел. Так как ни одно из рассматриваемых чисел не делится на 2, 3, или 5, то, соответственно, они не делятся и на их квадраты. Поэтому проверка, что число не кратно квадрату простого числа, не включает 2², 3², и 5².

Сегментация


 Для уменьшения требований к памяти «просеивание» производится порциями (сегментами, блоками), размер которых составляет примерно N.

Предпросеивание


 Для ускорения работы алгоритм игнорирует все числа, которые кратны одному из нескольких первых простых чисел (2, 3, 5, 7, \ldots). Это делается путём использования стандартных структур данных и алгоритмов их обработки, предложенных ранее Полом Притчардом . Они известны под названием . Количество первых простых чисел выбирается в зависимости от заданного числа N. Теоретически предлагается брать первые простые примерно до logN. Это позволяет улучшить асимптотическую оценку скорости алгоритма на множитель O(1loglogN). При этом требуется дополнительная память, которая с ростом N ограничена как explogN. Увеличение требований к памяти оценивается как O(No(1)).
 Версия, представленная на сайте одного из авторов, оптимизирована для поиска всех простых чисел до миллиарда (log109log230=305,5), в ней исключаются из вычислений числа, кратные двум, трём, пяти и семи (2 × 3 × 5 × 7 = 210).

Оценка сложности


 По оценке авторов алгоритм имеет асимптотическую сложность O(NloglogN) и требует O(N12+o(1)) бит памяти. Ранее были известны алгоритмы столь же асимптотически быстрые, но требующие существенно больше памяти. Теоретически в данном алгоритме сочетается максимальная скорость работы при наименьших требованиях к памяти. Реализация алгоритма, выполненная одним из авторов, показывает достаточно высокую практическую скорость.
 В алгоритме используется два вида оптимизации, которые существенно повышают его эффективность (по сравнению с упрощённой версией).
 Ниже представлена реализация упрощённой версии на языке программирования C, иллюстрирующая основную идею алгоритма — использование квадратичных форм:
 \beginShaded \beginHighlighting[] \DataTypeTokint \NormalToklimit = \DecValTok1000\NormalTok; \DataTypeTokint \NormalToksqr_lim; \DataTypeTokbool \NormalTokis_prime[\DecValTok1001\NormalTok]; \DataTypeTokint \NormalTokx2, y2; \DataTypeTokint \NormalToki, j; \DataTypeTokint \NormalTokn;
 \CommentTok// Инициализация решета \NormalToksqr_lim = (\DataTypeTokint\NormalTok)sqrt((\DataTypeToklong \DataTypeTokdouble\NormalTok)limit); \KeywordTokfor \NormalTok(i = \DecValTok0\NormalTok; i <= limit; i++) is_prime[i] = \KeywordTokfalse\NormalTok; \NormalTokis_prime[\DecValTok2\NormalTok] = \KeywordToktrue\NormalTok; \NormalTokis_prime[\DecValTok3\NormalTok] = \KeywordToktrue\NormalTok;
 \CommentTok// Предположительно простые — это целые с нечётным числом \CommentTok// представлений в данных квадратных формах. \CommentTok// x2 и y2 — это квадраты i и j (оптимизация). \NormalTokx2 = \DecValTok0\NormalTok; \KeywordTokfor \NormalTok(i = \DecValTok1\NormalTok; i <= sqr_lim; i++) \{ \NormalTokx2 += \DecValTok2 \NormalTok* i - \DecValTok1\NormalTok; \NormalToky2 = \DecValTok0\NormalTok; \KeywordTokfor \NormalTok(j = \DecValTok1\NormalTok; j <= sqr_lim; j++) \{ \NormalToky2 += \DecValTok2 \NormalTok* j - \DecValTok1\NormalTok; \NormalTokn = \DecValTok4 \NormalTok* x2 + y2; \KeywordTokif \NormalTok((n <= limit) && (n \NormalTokis_prime[n] = !is_prime[n]; \CommentTok// n = 3 * x2 + y2; \NormalTokn -= x2; \CommentTok// Оптимизация \KeywordTokif \NormalTok((n <= limit) && (n \NormalTokis_prime[n] = !is_prime[n]; \CommentTok// n = 3 * x2 - y2; \NormalTokn -= \DecValTok2 \NormalTok* y2; \CommentTok// Оптимизация \KeywordTokif \NormalTok((i > j) && (n <= limit) && (n \NormalTokis_prime[n] = !is_prime[n]; \NormalTok\} \NormalTok\}
 \CommentTok// Отсеиваем кратные квадратам простых чисел в интервале [5, sqrt(limit)]. \CommentTok// (основной этап не может их отсеять) \KeywordTokfor \NormalTok(i = \DecValTok5\NormalTok; i <= sqr_lim; i++) \{ \KeywordTokif \NormalTok(is_prime[i]) \{ \NormalTokn = i * i; \KeywordTokfor \NormalTok(j = n; j <= limit; j += n) \{ \NormalTokis_prime[j] = \KeywordTokfalse\NormalTok; \NormalTok\} \NormalTok\} \NormalTok\}
 \CommentTok// Вывод списка простых чисел в консоль. \NormalTokprintf(\StringTok"2, 3, 5"\NormalTok); \KeywordTokfor \NormalTok(i = \DecValTok6\NormalTok; i <= limit; i++) \{ \CommentTok// добавлена проверка делимости на 3 и 5. В оригинальной версии алгоритма потребности в ней нет. \KeywordTokif \NormalTok((is_prime[i]) && (i \NormalTokprintf(\StringTok", \NormalTok\} \NormalTok\} \endHighlighting \endShaded

Версия алгоритма на языке Pascal


 \beginShaded \beginHighlighting[] \KeywordTokprogram \NormalTokatkin; \KeywordTokvar \NormalTokis_prime:\KeywordTokarray\NormalTok[\DecValTok1\NormalTok..\DecValTok10001\NormalTok] \KeywordTokof \DataTypeTokboolean\NormalTok; jj:\DataTypeTokint64\NormalTok; \KeywordTokprocedure \NormalTokdosieve(limit:\DataTypeTokint64\NormalTok); \KeywordTokvar \NormalToki,k,x,y:\DataTypeTokint64\NormalTok; n:\DataTypeTokint64\NormalTok; \KeywordTokbegin \KeywordTokfor \NormalToki:=\DecValTok5 \KeywordTokto \NormalToklimit \KeywordTokdo \NormalTokis_prime[i]:=\KeywordTokfalse\NormalTok; \KeywordTokfor \NormalTokx:=\DecValTok1 \KeywordTokto \NormalToktrunc(sqrt(limit)) \KeywordTokdo \KeywordTokfor \NormalToky:=\DecValTok1 \KeywordTokto \NormalToktrunc(sqrt(limit)) \KeywordTokdo \KeywordTokbegin \NormalTokn:=\DecValTok4\NormalTok*sqr(x)+sqr(y); \KeywordTokif \NormalTok(n<=limit) \KeywordTokand \NormalTok((n \KeywordTokmod \DecValTok12 \NormalTok= \DecValTok1\NormalTok) \KeywordTokor \NormalTok(n \KeywordTokmod \DecValTok12 \NormalTok= \DecValTok5\NormalTok)) \KeywordTokthen \NormalTokis_prime[n]:=\KeywordToknot \NormalTokis_prime[n]; \NormalTokn:=n-sqr(x); \KeywordTokif \NormalTok(n<=limit) \KeywordTokand \NormalTok(n \KeywordTokmod \DecValTok12 \NormalTok= \DecValTok7\NormalTok) \KeywordTokthen \NormalTokis_prime[n]:=\KeywordToknot \NormalTokis_prime[n]; \NormalTokn:=n\DecValTok-2\NormalTok*sqr(y); \KeywordTokif \NormalTok(x>y) \KeywordTokand \NormalTok(n<=limit) \KeywordTokand \NormalTok(n \KeywordTokmod \DecValTok12 \NormalTok= \DecValTok11\NormalTok) \KeywordTokthen \NormalTokis_prime[n]:=\KeywordToknot \NormalTokis_prime[n]; \KeywordTokend\NormalTok; \KeywordTokfor \NormalToki:=\DecValTok5 \KeywordTokto \NormalToktrunc(sqrt(limit)) \KeywordTokdo \KeywordTokbegin \KeywordTokif \NormalTokis_prime[i] \KeywordTokthen \KeywordTokbegin \NormalTokk:=sqr(i); \NormalTokn:=k; \KeywordTokwhile \NormalTokn<=limit \KeywordTokdo \KeywordTokbegin \NormalTokis_prime[n]:=\KeywordTokfalse\NormalTok; \NormalTokn:=n+k; \KeywordTokend\NormalTok; \KeywordTokend\NormalTok; \KeywordTokend\NormalTok; \NormalTokis_prime[\DecValTok2\NormalTok]:=\KeywordToktrue\NormalTok; \NormalTokis_prime[\DecValTok3\NormalTok]:=\KeywordToktrue\NormalTok; \KeywordTokend\NormalTok; \KeywordTokbegin \NormalTokdosieve(\DecValTok10000\NormalTok); \KeywordTokfor \NormalTokjj:=\DecValTok1 \KeywordTokto \DecValTok10000 \KeywordTokdo \KeywordTokif \NormalTokis_prime[jj] \KeywordTokthen \NormalTokwriteln(jj); \KeywordTokend\NormalTok. \endHighlighting \endShaded