Project Euler 187

プロジェクトオイラー
http://projecteuler.net/index.php

Q187.
1億より小さい2つの素数の積で表される数の個数。

まず、1万までの素数を求める。そこまでの素数の積で表される個数はすぐにわかる。一方の素数が1万以上の数の個数を求めるために、5000万未満の素数が必要となる。
エラトステネスのふるいを使えば速くなるが、メモリがそれなりに必要となる。どうやらPythonでは1要素あたり4バイト要るらしく、200MBを使う。True/Falseでも4バイト使う。試しに、1ビットで素数・素数でないを表そうと、32個の素数・非素数を32ビット整数で表そうとしたら、どうも多倍長整数になってしまうようで遅くなる。欲張らずに16個にすると、True/Falseのときの1.4倍の時間で収まった。

def is_prime(n):

    for p in primes:

        if p * p > n:

            return True

        elif n % p == 0:

            return False

    return True
def make_primes(limit):

    for n in xrange(3, limit + 1, 2):

        if is_prime(n):

            primes.append(n)
def find_prime(m, ip):

    for i in xrange(ip, -1, -1):

        if m >= primes[i]:

            return i
def sieve(n):

    a = [ 0 ] * ((n + 15) / 16)

    for p in primes:

        if p * p > n:

            break

        if p == 2:

            continue

        for m in xrange(3 * p, n, 2 * p):

            a[m>>4] |= 1 << (m & 15)

    return a
N = 10 ** 8

M = int((N + 0.5) ** 0.5)

primes = [ 2 ]

make_primes(M)

counter = len(primes) * (len(primes) + 1) / 2

primes_more = sieve(N / 2)

ip = len(primes) - 1

for p in xrange((M + 1) / 2 * 2 + 1, N / 2, 2):

    if (primes_more[p>>4] & (1 << (p & 15))) == 0:

        ip = find_prime(N / p, ip)

        counter += ip + 1

print counter