Project Euler 234

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

Q234.
ある4以上の整数nを考える。√n以下の最大の素数をlps(n)、√n以上の最小の素数ups(n)とする。lps(n)かups(n)のどちらか一方がnで割り切れるとき、nをsemidivisibleと呼ぶ。
999966663333を超えないsemidivisibleな整数の和を求めよ。

これはやさしい。
ある素数とその次の素数の間のsemidivisibleな整数の和は簡単に求まる。



from math import sqrt

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(n):
m = int(sqrt(n + 0.5))
for p in xrange(3, m + 1, 2):
if is_prime(p):
primes.append(p)

a = [ 0 ] * ((n + 15) / 16)
for p in primes:
if p * p > n:
break
if p == 2:
continue
for k in xrange(3 * p, n, 2 * p):
a[k>>4] |= 1 << (k & 15)

for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):
if (a[k>>4] & (1 << (k & 15))) == 0:
primes.append(k)

N = 999966663333
primes = [ 2 ]
make_primes(int(sqrt(N) * 1.3))

s = 0
prime_limit = int(sqrt(N + 0.5))
for n in xrange(len(primes) - 1):
p1 = primes[n]
p2 = primes[n+1]
if p2 > prime_limit:
diff = N - p1 * p1
n1 = diff / p1
x0 = N / p2 * p2
n2 = (x0 - p1 * p1) / p2
s += (p1 * n1 * (2 * p1 + n1 + 1)
+ 2 * x0 * (n2 + 1) - n2 * (n2 + 1) * p2) / 2
if x0 >= p1 * p2:
s -= 2 * p1 * p2
break
diff = p2 * p2 - p1 * p1
n1 = diff / p1
n2 = diff / p2
s += (p1 * n1 * (2 * p1 + n1 + 1)
+ p2 * n2 * (2 * p2 - n2 - 1)) / 2 - 2 * p1 * p2
print s