Project Euler 183

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

Q183.
自然数Nに対して、自然数kを変化させて、(N/k)^kの最大値を求める。このとき、D(N)は、この最大値が無限小数ならN、有限小数なら-Nを返す。D(5) + ... + D(10000)を求めよ。

特になんのひねりもなく解答できる。微分するとkはeでPが最大になることがわかるので、eをはさむ2つのkで比較すればよい。このときkが2と5しか因子を含まなければ有限小数
しかし、Pの計算で、オーバーフローしてしまったので、klog(N/k)で代用した。本当は、有限小数になるkはごく限られているので、その付近だけ計算すればいいが。



from math import e, log
from fractions import Fraction
from itertools import imap

def P(n, k):
r = float(n) / k
return k * log(r)

def div_pow(n, p):
while n % p == 0:
n /= p
return n

def is_finite(f):
d = f.denominator
d = div_pow(d, 2)
d = div_pow(d, 5)
return d == 1

def D(n):
k1 = int(n / e)
k2 = k1 + 1
P1 = P(n, k1)
P2 = P(n, k2)
if P1 > P2:
k = k1
else:
k = k2

if is_finite(Fraction(n, k)):
return -n
else:
return n

N = 10000
print sum(imap(D, xrange(5, N + 1)))