\(N\) 以下の素数を \(O(N/\log\log(N))\) で列挙するアルゴリズムである、wheel sieve （車輪の篩？）を紹介します。

エラトステネスの篩は、素数か判断する数を \(2\) と奇数に限定すれば、定数倍高速化できました。更に、素数か判断する数を \(2,3\) と \(6k+1,6k+5\) の形に限定することもできます。各 \(N\) に対して、計算量を最小化する最適な前列挙が存在するはずです。この考えを元にしたのが wheel sieve です。

アルゴリズム

素数を小さい方から順に \(p_1,p_2,\ldots\) と置きます。

まず、\(w:=p_1 p_2 \cdots p_m \leq {\sqrt{N}}\) となる素数 \(p_1,p_2,\ldots,p_m\) を最大個求めます。\(N\) 以下の素数を列挙するのに掛かる計算量はエラトステネスの篩で \(O(N\log\log(N))\) ですから、\(p_1,\ldots,p_m\) は \(O(\sqrt{N}\log\log(N))\) で列挙できます。また、同時に \(w\) 以下の自然数のうち、\(w\) と互いに素なものの集合 \(R=\{r_1,r_2,\ldots,r_s\}\) が得られます（昇順とする）。

\(w\) と互いに素な自然数の集合は \(S:=\{r+kw : k\in\mathbb{Z}_{\geq 0},r \in R\}\) の形で書けます。\(p_1,p_2,\ldots,p_m\) 以外に素数判定したい数は全てこの形です。\(|R|/w=O(1/\log\log(N))\) なので（フェルミ推定）、素数判定する数は \(O(N/\log\log(N))\) 個あります。

\(p_1,p_2,\ldots,p_m\) の合成数は \(S\) に含まれないので、素数 \(p\) で篩うとき、\(\{sp : s \in S\}\) の合成数だけ見ればいいです。素数 \(p\) で篩う数は \(O\left(\frac{N}{p\log\log(N)}\right)\) 個あるので、全体では \(O(N)\) です（フェルミ推定）。

更に高速化しましょう。同じ数が二度篩われないために、\(\{ps : s \in S, \mathrm{lpf}(ps)=p \}\) の合成数（\(\mathrm{lpf}(x):=x\)の最小素因数）だけ見たいです。これは逆に \(s \in S\) を固定し、\(\{ps : p \leq \mathrm{lpf}(s)\}\) の合成数を篩えばよいです。計算量は \(O(N/\log\log(N))\) になります。

実装

\(O(N / \log\log (N))\) ：Python での実装は次の通りです。

import itertools
import math
import operator
N = 10 ** 5  # [1:N]に含まれる素数を列挙

tiny_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]  # small primes
ws = list(itertools.accumulate(tiny_primes, func=operator.mul))
m = 0
while ws[m + 1] <= math.sqrt(N):
    m += 1
w = ws[m]
tiny_primes = tiny_primes[:m + 1]
is_coprime = [True] * w
is_coprime[0] = False
for p in tiny_primes:
    for i in range(p, w, p):
        is_coprime[i] = False

coprimes = [i for i in range(1, w) if is_coprime[i]]
m = len(coprimes)
id = [-1] * w
for i in range(m):
    id[coprimes[i]] = i

dt = [(w + coprimes[(i + 1) % m] - coprimes[i]) % w for i in range(m)]


def get_id(v):
    return v // w * m + id[v % w]


def get_val(v):
    return v // m * w + coprimes[v % m]


n = (N + w - 1) // w * m
pos = coprimes[1]
primes = []

lpf = [get_val(i) for i in range(n)]  # 最小素因数
for i in range(1, n):
    if pos > N:
        break
    if lpf[i] == pos:
        primes.append(pos)
    for p in primes:
        if p > lpf[i] or p * pos > N:
            break
        lpf[get_id(pos * p)] = p
    pos += dt[i % len(dt)]
primes = tiny_primes + primes  # [1:N] に含まれる素数

\(O(N)\) の実装は最後のループを次で置き換えればよいです。区間篩のようにすると、メモリ \(O(\sqrt{N})\) にもできます（見つけた素数はハードディスクに書き込む）。\(O(N/\log\log(N))\) にすると、メモリ \(\Theta(N/\log\log(N))\) から落ちないので、実は一長一短です。

is_prime = [True] * n
is_prime[0] = False
for i in range(1, n):
    if pos > N:
        break
    if is_prime[i]:
        primes.append(pos)
        j, multiplier = 1, coprimes[1]
        while pos * multiplier <= N:
            is_prime[get_id(pos * multiplier)] = False
            multiplier += dt[j % m]
            j += 1
    pos += dt[i % m]
primes = tiny_primes + primes

参考文献

Explaining Wheel Sieve.