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
$
1 コメント

A. 試し割り法

キーワード: , ,
　試し割り法 (Trial division) は, 素因数分解アルゴリズムの中では最も効率の悪い方法ですが, 簡単に理解することができます. 試し割り法の基本的な考えは, 素因数に分解される整数 $n$ を $n$ より小さい整数で, ただ割っていくことです.

　$n$ の素因数を出力するだけなら, 

def junk(n):
    """"""
    for i in range(2, n):
        while n % i == 0:
            n //= i
            print(i)
    if n > 1:
        print(n)

最小の素因数 $2$ から繰り返し割っていけば, 割り切れる整数に素数が残っていくことは明らかなので、あとは残った $n$ が $n > 1$ であるならば, それが最後の素因数となります.

　$n$ を素因数分解したとき, $\sqrt{n}$ 以上の素因数は多くても一つしか含まれないので, 最大の因数は $\sqrt{n}$ としていいことがわかります。例えば $n=360$ のとき, 最大の因数は $\sqrt{n}\ \fallingdotseq\ 18.97367$ となります. ここで $\lfloor \sqrt{n} \rfloor = 18$ としてしまうと誤差がでるため, $\lfloor \sqrt{n} \rfloor + 1$ もしくは $\lceil \sqrt{n} \rceil$ とします.

Python による実装 (1) 

def trial_division(n):
    """"""
    if n == 1:
        return [1]
    factors = []
    for i in range(2, int(n**0.5) + 1 + 1):
        while n % i == 0:
            n //= i
            factors.append(i)
    if n > 1:
        factors.append(n)
    return factors

実行例 

trial_division(360)
[2, 2, 2, 3, 3, 5]

計算量

$\pi(2^{n/2})$

0 コメント