プロジェクトオイラー
http://projecteuler.net/index.php
Q221.
A = pqr, 1/A = 1/p + 1/q + 1/r
となる正の整数Aの15万番目を求めよ。
整理して、
(p + q)(p + r) = p2 + 1
ここで、p > 0, q, r < 0。
右辺を素因数分解してもいいが、やはりコストが大きい。Q216と同じような手法で、エラトステネスのふるい的に素因数分解していく。
def numerize(b, n):
m = 1
for p, e in b:
m *= p ** (n % (e + 1))
n /= (e + 1)
return mdef gen_divisors(b):
def mul(x, y): return x * y
n = reduce(mul, map(lambda t: t[0] ** t[1], b))
m = reduce(mul, map(lambda t: t[1] + 1, b))
for i in xrange(m):
k = numerize(b, i)
if k * k <= n:
yield kdef proc_p(c, p, n):
if n + p > N:
return
for m in xrange(n + p, N + 1, p):
q = c[m][-1][0]
e = 0
while q % p == 0:
q /= p
e += 1
c[m][-1] = (p, e)
if q != 1:
c[m].append( (q, 1) )N = 15 * 10 ** 4
ary = [ ]
set_primes = set()
c = [ [ (n * n + 1, 1) ] for n in xrange(N + 1) ]
for n in xrange(1, N + 1):
p = c[n][-1][0]
if p in set_primes:
continue
proc_p(c, p, n)
if p != 2:
proc_p(c, p, p - n)
set_primes.add(p)for p in xrange(1, N + 1):
n = p * p + 1
for d in gen_divisors(c[p]):
q = -p - d
r = -p - n / d
A = p * q * r
ary.append(A)
ary.sort()
print ary[N-1]