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