プロジェクトオイラー
http://projecteuler.net/index.php
Q108.
1/x + 1/y = 1/n の解(x, y > 0)の個数が1000を超える最小のn
(x - n)(y - n) = n2となるから、n2の約数を数えればよい。あとから見るように対象となるnは偶数だから、重複を考えて、約数は2001以上あればよい。
nを1から順に見ていっても求まるが、時間がかかる。逆に素因数分解した形から見ていく。
m = Πpieiなら、mの約数の個数はΠ(ei + 1)となる。約数の個数に素数の値自体は関係ないから、素数は小さいものから使い、eも降順とすればよい。n2を素因数分解したとき、e+1は3以上だから、37 = 2187 > 2001だから、素数はせいぜい7つあればよい。すなわち、{ 2, 3, 5, 7, 11, 13, 17 }を使う。分割数のときと同じようにギリギリ2001以上になるeのリストを発生させ、最も小さいnを求める。ただし、使う素数の個数が少ないときに、元の数が非常に大きくなる可能性があるので、そこには配慮して、なるべく大きな数の計算が行われないように工夫した。
from math import logdef make_primes(n):
primes = [ 2 ]
m = 3
while len(primes) < n:
for p in primes:
if p * p > m:
primes.append(m)
break
elif m % p == 0:
break
m += 2
return primesdef div(a, b):
c = { }
for p, e in a.iteritems():
c[p] = e
if b.has_key(p):
c[p] -= b[p]
for p, e in b.iteritems():
if not c.has_key(p):
c[p] = -e
return cdef key_with_max_value(a):
return reduce(lambda x, y: max(x, y, key=lambda x: abs(x[1])),
a.iteritems())[0]def less(a, b):
c = div(a, b)
pmax = key_with_max_value(c)
num = 1
den = 1
for p, e in c.iteritems():
if p != pmax:
if e > 0:
num *= p ** e
else:
den *= p ** -e
emax = c[pmax]
if emax > 0:
while emax > 0:
num *= pmax
if num > den:
return False
emax -= 1
return True
else:
while emax < 0:
den *= pmax
if num < den:
return True
emax += 1
return Falsedef Min(a, b):
if less(a, b):
return a
else:
return bdef gen_div(n, m, r = 0):
if r == 0:
r = n
elif r > n:
r = n
if m == 1:
if n <= r and n & 1 and n >= 3:
yield [ n ]
elif n >= 3:
for k in xrange(3, r + 1, 2):
for a in gen_div((n + k - 1) / k, m - 1, k):
yield [ k ] + adef exp_to_num(a):
m = { }
for i in range(len(a)):
m[primes[i]] = a[i] / 2
return mdef min_number_core(M, m):
return reduce(Min, map(exp_to_num, gen_div(M, m)))def min_number(M):
a = min_number_core(M, 1)
for i in range(2, 8):
b = min_number_core(M, i)
a = Min(a, b)
return reduce(Min, map(lambda i: min_number_core(M, i),
range(1, max_num_primes + 1)))def val(a):
def pow(x): return x[0] ** x[1]
return reduce(lambda x, y: x * y, map(pow, a.iteritems()))N = 1000
M = 2 * N + 1
max_num_primes = int(log(M) / log(3)) + 1
primes = make_primes(max_num_primes)print val(min_number(M))