Pythonで素数判定（3）

数学 アルゴリズム Python

エラトステネスのふるいの速度を確認します。まずは最も単純なもの。

def sieve1(N):
    a = [ True ] * (N + 1)
    for p in takewhile(lambda p: p * p <= N, (n for n in count(2) if a[n])):
        for k in xrange(p * 2, N + 1, p):
            a[k] = False
    return [ n for n in xrange(2, N + 1) if a[n] ]

t0 = time.clock()
N = 10 ** 6
primes = sieve1(N)
print len(primes)
print time.clock() - t0

78498
0.242247248827

外側のループなのでtakewhileを使ってもだいじょうぶ。
次に、2の倍数は最初から無視する方法。配列は、1, 3, 5, …が素数性を表しています。

def sieve2(N):
    M = (N + 1) / 2
    a = [ True ] * M
    for p in takewhile(lambda p: p * p <= N, (n * 2 + 1 for n in count(1) if a[n])):
        for k in xrange(p * 3 / 2, M, p):
            a[k] = False
    return [ 2 ] + [ n * 2 + 1 for n in xrange(1, M) if a[n] ]

78498
0.104205804153

最後に、3の倍数も無視する方法です。すなわち、用意する配列は1, 5, 7, 11, …の素数性を表しています。このように等差数列になっていないので少しめんどうです。等差数列が2つあるので別々に考えなければなりません。例えばp = 7の倍数をFalseにして行くことを考えます。この配列の3番目（インデックスだとk = 2）が7なので、これが乗っている等差数列（1, 7, 13, …）ではk + 2p, k + 4p, …がp = 7の倍数となるのでFalseにします。もう一つの等差数列の7の倍数の一つは-7です。pと-pは必ず別々の等差数列に乗っています。そのインデックスは-k - 1となるので、-k - 1 + 2p, -k - 1 + 4p, …がpの倍数になります。

def sieve3(N):
    M = N / 6 * 2 + 1 if (N - 1) % 6 < 4 else (N + 1) / 3
    a = [ True ] * M
    for p in takewhile(lambda p: p * p <= N,
                    (n * 3 + (n & 1) + 1 for n in count(1) if a[n])):
        n = p / 3
        for k in xrange(n + p * 2, M, p * 2):
            a[k] = False
        for k in xrange(-n - 1 + p * 2, M, p * 2):
            a[k] = False
    return [ 2, 3 ] + [ n * 3 + (n & 1) + 1 for n in xrange(1, M) if a[n] ]

78498
0.0756187572501

そこそこ速いですね。
エラトステネスのふるいはPythonでは100万までなら素数による試し割りより30倍以上速いです。これは大きくなると差が広がります。計算量を計算すればわかります。エラトステネスのふるいの計算量は簡単です。最初のコードで考えて、p の倍数はn/p 個あるので計算量は、

\sum_{p \le \sqrt{n}}{\frac{n}{p}}

これを積分で近似します。素数定理から素数の密度は1 / log pだから、

\int_2^{\sqrt{n}}{\frac{n}{p\log{p}}dp} = [n\log{\log{p}}]_2^{\sqrt{n}} = n(\log{\log{n}} - \log{2} - \log{\log{2}})

計算量はO(n log log n)となります。