A. サンダラムの篩

　1934 年, インドの数学科の学生サンダラム (Sundaram) によって, 素数を求めるアルゴリズム「サンダラムの篩」が発見されました.

　このアルゴリズムでは, $1$ から $n$ までの整数から

$i + j + 2ij \leqq n \quad (i,j \in N, \ 1 \leqq i \leqq j)$

の形に表わすことができる整数を取り除きます. この形で表わせる任意の整数を $k$ とし $2k + 1$ とすれば, 奇数の合成数が得られます (後述). これがサンダラムの篩の核心になります.

素数と合成数

　$n > 1$ の整数は次の二種類に分かれます.

　　素数　　　$1$ と $n$ 以外の約数がない (例えば $2, 31, 73$ など).
　　合成数　　$1$ と $n$ 以外の約数がある (例えば $8, 49, 72$ など).

また合成数は次の三種類に分かれます.

　　偶数の合成数　　偶数と偶数の積 (例えば $2 \times 4 = 8$ など).
　　偶数の合成数　　偶数と奇数の積 (例えば $8 \times 9 = 72$ など).
　　奇数の合成数　　奇数と奇数の積 (例えば $7 \times 7 = 49$ など).

したがって整数 $n > 1$ は, 素数または合成数で, 合成数だけを篩い落せば素数が残されるという訳です.

理論的な根拠

　まず正の整数を

$1, \ 2, \ 3, \ \cdots\cdots$

のように大きさの順序に並べ, それぞれの整数を $2n+1$ に代入すれば

$3, \ 5, \ 7, \ \cdots\cdots$

で, 偶数の合成数だけを取り除いた数列が得られます. 残ったのは素数と奇数の合成数からなる数列で, あとは奇数の合成数を取り除くことができれば, 素数だけが残る数列が得られることになります.

　さて奇数の合成数は

$(2i+1)(2j+1) \quad (1 \leqq i \leqq j)$

の形で表わすことができます. ここで

$(2i+1)(2j+1)=2(i+j+2ij)+1$

とすれば

$k=i+j+2ij$

で

$(2i+1)(2j+1)=2k+1$

を得ます. 言い換えれば $k$ を

$k=4, \ 7, \ 10, \ 12, \ 13, \ 16, \ 17, \ \cdots\cdots$

のように大きさの順序に並べ, それぞれの整数を $2k + 1$ に代入すれば奇数の合成数が得られます.

　いま正の整数を

$1, \ 2, \ 3, \ \cdots\cdots$

のように大きさの順序に並べ, それらの中から $k$ を取り除けば

$1, \ 2, \ 3, \ 5, \ 6, \ 8, \ 9, \ 11, \ 14, \ 15, \ \cdots\cdots$

で, $k$ を篩った数列が得られます. この数列のそれぞれの整数を $2n+1$ に代入すれば

$3, \ 5, \ 7, \ 11, \ 13, \ 17, \ 19, \ 23, \ 29, \ 31, \ \cdots\cdots$

が得られます. ここで $k$ を $2n+1$ にしたと仮定すれば, それは奇数の合成数になるので, 最後に得られた数列は偶数の合成数と奇数の合成数を篩った数列, すなわち素数の数列 (ただし $2$ を除く) となります.

C による実装

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX_PRIMES 1000000

void
sundaram(void)
{
    int *sieve = (int *) calloc(MAX_PRIMES/2, sizeof(int));
    if (sieve == NULL)
        return;

    int i, j, k;
    for (i = 1; i < sqrt(MAX_PRIMES); i++)
        for (j = i; (k = i+j+2*i*j) < MAX_PRIMES/2; j++)
            sieve[k] = 1;

    printf("2\n");
    for (i = 1; i < MAX_PRIMES/2; i++)
        if (sieve[i] == 0)
            printf("%d\n", 2*i+1);

    free(sieve);
}

int
main(void)
{
    sundaram();

    return EXIT_SUCCESS;
}

実行例

$ time ./a.out | wc -l
  664579
    0m1.11s real     0m1.06s user     0m0.02s system
$

Cython による実装

from libc.stdlib cimport calloc, free

cdef enum:
    MAX_PRIMES = 1000000

cdef sundaram():
    cdef bint *sieve = <bint *> calloc(MAX_PRIMES//2, sizeof(bint))

    cdef i, j, k
    i = 1
    while i < MAX_PRIMES ** 0.5:
        j = i
        k = i + j + 2 * i * j
        while k < MAX_PRIMES//2:
            sieve[k] = 1
            j += 1
            k = i + j + 2 * i * j
        i += 1

    print(2)
    i = 1
    while i < MAX_PRIMES//2:
        if sieve[i] == 0:
            print(2*i+1)
        i += 1

    free(sieve)

sundaram()

実行例

$ time python run.py | wc -l
  664579
    0m10.96s real     0m10.89s user     0m0.02s system
$