Project Euler 2（2）

この程度の大きさの数なら問題ありませんが、実数の計算にはどうしても誤差がつきまといます。しかし、今回の計算は誤差無しで計算することができます。

(a + b√5) / c -> (a, b, c)

と表すことにしましょう。そうすると、

(a, b, c) + (d, e, f) = (af + dc, bf + ec, cf)

となります。ただし、このあと約分が必要です。乗算は、

(a, b, c) * (d, e, f) = (ad + 5be, ae + bd, cf)

逆数は、

(a, b, c)^-1 = (ac, -bc, a² - 5b²)

これで、四則演算はできることになります。
また、整数との大きさの比較は、

(a, b, c) - (d, 0, 1) = (a - cd, b, c)

だから、(a - cd)²と5b²の比較をすればよいです。
logの計算だけはさすがにできません。

from fractions import gcd

from math import sqrt, log
def lcm(x, y):

    return x * y / gcd(x, y)
def f_reduce(x, y):

    c = lcm(x[2], y[2])

    if c != x[2]:

        d = c / x[2]

        return (x[0] * d, x[1] * d, x[2] * d)

    else:

        return x
def cancel_down(x):

    c = reduce(gcd, x)

    if c != 1:

        return (x[0] / c, x[1] / c, x[2] / c)

    else:

        return x
def f_add(x, y):

    x = f_reduce(x, y)

    y = f_reduce(y, x)

    return cancel_down( (x[0] + y[0], x[1] + y[1], x[2]) )
def f_sub(x, y):

    x = f_reduce(x, y)

    y = f_reduce(y, x)

    return cancel_down( (x[0] - y[0], x[1] - y[1], x[2]) )
def f_mul(x, y):

    z = (x[0] * y[0] + x[1] * y[1] * D, x[0] * y[1] + x[1] * y[0], x[2] * y[2])

    return cancel_down(z)
def f_inv(x):

    z = (x[0] * x[2], -x[1] * x[2], x[0] * x[0] - x[1] * x[1] * D)

    return cancel_down(z)
def f_div(x, y):

    return f_mul(x, f_inv(y))
def f_value(x):

    return (x[0] + x[1] * sqrt(D)) / x[2]
def f_pow(x, n):

    if n == 0:

        return (1, 0, 1)

    

    y = f_pow(x, n / 2)

    sq = f_mul(y, y)

    if n % 2 == 1:

        return f_mul(x, sq)

    else:

        return sq
def round_sqrt(n):

    x, y = 0, 1

    while x != y:

        x, y = y, (y + (n + y - 1) / y) / 2

    return x
def f_round(x):

    n = x[1] ** 2 * D

    return (x[0] + round_sqrt(x[1] ** 2 * D) + x[2] / 2) / x[2]
N = 4 * 10 ** 6

D = 5

a = (1, -1, 2)

b = (1, 1, 2)

f0, f1 = (1, 0, 1), (2, 0, 1)

c = f_div(f_mul(b, f_sub(f1, f_mul(f0, a))), f_sub(b, a))

m_max = int(log( (N + 0.5) / f_value(c)) / (3 * log(f_value(b))))

x = f_div(f_mul(c, f_sub(f_pow(b, 3 * m_max + 3), (1, 0, 1))),

        f_sub(f_pow(b, 3), (1, 0, 1)))

print f_round(x)