Project Euler 110

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

Q110.
1/x + 1/y = 1/n の解(x, y > 0)の個数が400万を超える最小のn

108とほとんど同じ。少し工夫した。
例えば、素因数分解したときに素数が2つというのはありえない。約数の個数が100として、2か3のどちらかの指数が9以上となる。より厳しい2の指数が9以上とすると、同じ約数の個数でも、29 > 2451だから、素数が3つ以上あったほうが小さくなる。このような考えに基づくと、素数は8個以上必要であることがわかる。



from itertools import count
from math import log

def 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 primes

def 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 c

def 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):
if a == 0:
return False
elif b == 0:
return True
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 False

def Min(a, b):
if less(a, b):
return a
else:
return b

def 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 ] + a

def exp_to_num(a):
m = { }
for i in range(len(a)):
m[primes[i]] = a[i] / 2
return m

def min_number_core(M, m):
print M, m
return reduce(Min, map(exp_to_num, gen_div(M, m)), 0)

def min_number(M):
min_num_primes = calc_min_num_primes(M)
return reduce(Min, map(lambda i: min_number_core(M, i),
range(min_num_primes, max_num_primes + 1)), 0)

def val(a):
def pow(x): return x[0] ** x[1]
return reduce(lambda x, y: x * y, map(pow, a.iteritems()))

def calc_min_num_primes(M):
for n in count(2):
min_e = int(M ** (1.0 / n))
if (min_e + 1) / 2 * log(2) < log(primes[n]):
return n

N = 4 * 10 ** 6
M = 2 * N + 1
max_num_primes = int(log(M) / log(3)) + 1
primes = make_primes(max_num_primes)
print max_num_primes

print val(min_number(M))