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)))