Pythonで素数判定（2）

昨日ひどかったので考え直してみました。要するに試し割りの回数がわかればいいのです。もちろんわからないので、n前後の試し割りの回数の平均を求めます。ガウスの故事に倣って1000個ずつ取って平均します。なるべく長いスパンで見たいので、平均を取る範囲の始まりを1.2倍ずつしていきます。ただし、最初は1000から始まって2倍ずつします。32ビット整数の範囲ギリギリまで計算しましたが、計算時間は一瞬でした。

from itertools import *
from math import sqrt

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

def pows(first, last, ratio):
    n = first
    while n < last:
        yield n
        n = int(n * ratio)

def mean_num_divs(first, last):
    a = [ len(primes) ] * (last - first)
    for k, (p, p_next) in enumerate(izip(primes, islice(primes, 1, None))):
        for l in xrange(max(first, p * p), min(last, p_next * p_next)):
            a[l-first] = k + 1
    
    b = [ last ] * (last - first)
    for k, p in enumerate(takewhile(lambda p: p * p < last, primes)):
        for l in xrange((first + p - 1) / p * p, last, p):
            b[l-first] = min(b[l-first], k + 1)
    
    return sum(min(m, n) for m, n in izip(a, b)) / float(last - first)

N = 21 * 10 ** 8
L = 1000
primes = sieve(int(sqrt(N)))
gen_starts = chain(pows(1000, 10000, 2.0), pows(10000, N, 1.2))
for start in gen_starts:
    print start + L / 2, mean_num_divs(start, start + L)

対数軸表示してみましたが、直線に乗りませんね。

しかし、なぜか(logn)³を掛けると、かなりきれいに直線に乗ります。