Project Euler P.243

$R(d) < \dfrac{\phi(d)}{d}$ 임을 이용해서 $ \dfrac{\phi(d_1)}{d_1} < R(d) < \dfrac{15499}{94744} < \dfrac {\phi(d_2)}{d_2} $ 인 $d_1$을 찾고 $d_1 < d < d_2$ 를 이용해서 $d$를 찾는다.

#!/usr/bin/env python
#-*- coding:utf-8 -*-
import operator

def combinations_with_replacement(iterable, r):
  pool = tuple(iterable)
  n = len(pool)
  if not n and r:
    return
  indices = [0] * r
  yield tuple(pool[i] for i in indices)
  while True:
    for i in reversed(range(r)):
      if indices[i] != n - 1:
        break
    else:
      return
    indices[i:] = [indices[i] + 1] * (r - i)
    yield tuple(pool[i] for i in indices)

def factors(n):
  l = []
  if n % 2 == 0:
    l.append(2)
    while n % 2 == 0:
      n /= 2
  i = 3
  while n > 1:
    if n % i == 0:
      l.append(i)
    while n % i == 0:
      n /= i
    i += 2
  return l

def phi(n):
  f = factors(n)
  p = n
  for factor in f:
    p = p * (factor - 1) / float(factor)
  return p

def gen_primes():
  D = {}
  q = 2

  while True:
    if q not in D:
      yield q
      D[q * q] = [q]
    else:
      for p in D[q]:
        D.setdefault(p + q, []).append(p)
      del D[q]
    q += 1

def r(d):
  d = float(d)
  return (1.0 / (d-1)) * phi(int(d))

def solve(limit = (15499.0 / 94744.0)):
  primes = gen_primes()
  rr = 1.0
  prime = primes.next()
  pl = [ prime ]
  while rr >= limit:
    prime = primes.next()
    pl.append(prime)
    tmp = r( reduce(operator.mul, pl) )
    rr = tmp
  pl.pop()

  additionals = []
  for i in range(1, len(pl)):
    for a in combinations_with_replacement(pl, i):
      m = reduce(operator.mul, a)
      if m not in additionals:
        additionals.append(m)
  additionals.sort()
  for m in additionals:
    d = reduce(operator.mul, pl) * m
    if r(d) < limit:
      print 'Answer: %s' % d
      return

if __name__ == '__main__':
  solve()