Sorensen's Incremental Sieve は素数を順番に列挙するアルゴリズムです．このアルゴリズムは incremental, compact という 2 つの性質をもっています．incremental というのは $n$ までの素数を列挙したとき $n+1$ に対する処理が高速（$ t(n) $ 回の算術演算）であること，compact というのは $n$ まで素数を列挙するのに $O^*(\sqrt{n})$ のメモリしか使わないことを指します．Sorensen's Incremental Sieve は $O(\log n / \log \log n)$-incremental time，および $O(\sqrt{n} \log n)$-bit spaceのアルゴリズムです．

方針

素数をバケットのリストで管理します．次に処理する数 $n$ に対して $\sqrt{n}$ 未満の素数 $p$ を列挙し，$(-n) \mod p$ 番目のバケットに挿入します．このとき 0 番目のバケットが空であることと $n$ が素数または素数の平方数であることが同値となるので素数判定が $O(1)$ 回の算術演算でできます．また，数を $n$ から $n+1$ に変化させたときのバケットたちの変化は簡単に管理できます．

バケットたちの変化を具体例で見ましょう．$n = 26$ とします．このとき考える素数は $p = 2, 3, 5$ であり，それらによる $n$ の剰余は $0$, $2$, $1$ です．したがってそれらを 0 番目, 1 番目, 4 番目のバケットに挿入することで (2) (3) () () (5) が得られます．ここから順にカウントアップしていくとバケットたちの様子は次のようになります．

26: (2) (3) () () (5)
27: (3) (2) () (5) ()
28: (2) () (3,5) () ()
29: () (2,3,5) () () ()
30: (2,3,5) () () () ()
31: () (2) (3) () (5)
32: (2) (3) () (5) ()
33: (3) (2) (5) () ()

これを見ればわかる通り，「各ステップで全ての数はバケットを右から左に 1 つずつ移動する．ただし 0 番目のバケットの数は $p - 1$ に移動する」という法則で変化しています．あとはこれを実装するだけです．すべてのバケットを同時に移動させるのは添字にオフセットを加えればよいので，普通のバケットの配列を用いて実装できます．

実装例

#[derive(Debug)]
pub struct IncrementalSieveSorensen {
    n: usize,
    s: usize,
    r: usize,
    stacks: Vec<Vec<usize>>,
    pos: usize,
}
impl IncrementalSieveSorensen {
    pub fn new(start: u64) -> Self {
        let start = start.max(2) as usize;
        let r = (start as f64).sqrt() as usize + 1;
        let mut is_prime = vec![true; r];
        let mut stacks = vec![vec![]; r + 2];
        for p in 2..r {
            if is_prime[p] {
                let j = (p - start % p) % p;
                stacks[j].push(p);
                for q in (p * p..r).step_by(p) {
                    is_prime[q] = false;
                }
            }
        }
        Self {
            n: start,
            s: r * r,
            r,
            stacks,
            pos: 0,
        }
    }

    fn is_next_prime(&mut self) -> bool {
        let m = self.stacks.len();
        let mut is_prime = true;
        while let Some(p) = self.stacks[self.pos].pop() {
            self.stacks[(p + self.pos) % m].push(p);
            is_prime = false;
        }
        if self.n == self.s {
            if is_prime {
                self.stacks[(self.r + self.pos) % m].push(self.r);
                is_prime = false;
            }
            self.r += 1;
            self.s = self.r * self.r;
        }
        self.n += 1;
        self.pos += 1;
        if self.pos == m {
            self.pos = 0;
            self.stacks.push(vec![]);
            self.stacks.push(vec![]);
        }
        is_prime
    }
}
impl Iterator for IncrementalSieveSorensen {
    type Item = u64;

    fn next(&mut self) -> Option<u64> {
        while !self.is_next_prime() {}
        Some(self.n as u64 - 1)
    }
}

TM's memo

Sorensen's Incremental Prime Sieve

方針

実装例