ビット配列による素数表の作成をpythonで書いたら激しく遅かった件

以前のエントリで書いたエラトステネスの篩をpythonで実装しようとして引っかかった。

まず、C言語と違ってpythonには符号なし整数が無い上に、整数と多倍長整数がシームレスなので1を32桁左にビットシフトすると0ではなく2^32になってしまう。これは、ひとつあたりの整数で表す素数表の大きさは24として解決した。

問題は、ビット配列による素数表の作成に非常に時間がかかってしまうということ。以前C言語で実装したときは、1e7までの素数表を0.728秒で作成できたのだが、同じ環境で今回のプログラムを走らせると97.601秒かかる。

おそらく、C言語と違って配列の値を上書きする形でビット演算をしても新しく整数オブジェクトが生成されるために激しく遅くなっているのだと思う。hotshotモジュールでプロファイルをとってもpython内部の挙動のプロファイルはできないと思われるので、正確な原因を突き止める手段が見つからない。ちなみにhotshotを使ってtest関数をプロファイルしたら全体の95％の時間がsieve関数の実行に使われていた。

# coding: utf8

# ビット配列によるエラトステネスの篩。
# 一つの整数に24個記録し、さらに奇数のみを記録する。
# 32にしなかったのは、Pythonでは長精度整数と通常の整数がシームレスなので
# 扱いが面倒だから

ENTRY_PER_INT = 24

def sieve(u):
    from math import ceil, sqrt
    table = [2**ENTRY_PER_INT-1] * (u / ENTRY_PER_INT / 2 + 1)

    b = int( ceil(sqrt(u)) ) + 1
    # 奇数だけ除去
    for i in xrange(3, b, 2):
        for j in xrange(i*3, u+1, i*2):
            n = j/2 - 1
            table[n / ENTRY_PER_INT] &= ~(1 << n % ENTRY_PER_INT)
    return table

def call(upper):
    # 素数表とそれに基づいた素数判定関数と走査に使うジェレネータ関数を返す
    table = sieve(upper)
    def isprime(n):
        if n==2: return True
        if n < 2: return False
        if (n & 1) == 0: return False
        if n > upper: raise ValueError, "prime table does not have such large val only below %d; passed %d" % (upper, n)
        idx = n/2 - 1
        p = table[idx / ENTRY_PER_INT] & (1 << idx % ENTRY_PER_INT)
        return p != 0
    def generator_factory():
        yield 2
        for i in xrange(3, upper, 2):
            idx = i/2 - 1
            p = table[idx / ENTRY_PER_INT] & (1 << idx % ENTRY_PER_INT)
            if p != 0:
                yield i
    return table, isprime, generator_factory

def test():
    from time import time
    for i in [4,5,6,7]:
        s = time()
        t,isp,g = call(10**i)
        e = time()
        print i, e-s

if __name__ == '__main__':
    test()

実行結果：

4 0.0811238288879
5 1.14158892632
6 11.223911047
7 97.6018071175