ペル方程式（2）

ペル方程式の最小解を求めます。

x² - Dy² = 1

これを変形すると、

x / y = √(D + 1 / y²)

となり、x / y は√Dの近似になっています。x, yが大きくなればなるほど√Dに近づいていきます。
√Dを連分数で表したとき、途中で打ち切ったものの内でペル方程式を満たすものがあります。実はその最小がペル方程式の最小解にもなります。例えば√2を連分数で表すと[1;2,2,...]なので、

1 ×
1 + 1 / 2 = 3 / 2 ○

(x, y) = (3, 2)が最小解となります。
√Dを連分数で表す計算方法を考えましょう。通常の連分数計算を行います。連分数計算を途中で打ち切ると次の形になります。

(b_n + c_n√D) / d_n

これを超えない整数は、

a_n+1 = [ (b_n + [c_n√D]) / d_n ]

これを引いて逆数を取ると、

(-d_n(b_n - a_n+1d_n) + d_nc_n√D) / (a_n+1²D - (b_n - a_n+1d_n)²)

このa_nを使って連分数を作って、ペル方程式を満たすかどうかを調べます。
本当はもっとうまい方法があるようですが、これでも十分でしょう。

Project Eulerに必要な用語集
http://d.hatena.ne.jp/inamori/20091225/p1

from math import sqrt

from fractions import gcd
def int_sqrt(n):

    return int(sqrt(n))
def gen_cfraction_coeff(D):

    b, c, d = 0, 1, 1

    while True:

        a = (b + int_sqrt(c * c * D)) / d

        yield a

        b -= a * d

        b, c, d = -b * d, c * d, c * c * D - b * b

        div = gcd(d, gcd(c, b))

        if div != 1:

            b, c, d = b / div, c / div, d / div
def is_satisfied_equation(x, y, D):

    return x * x - D * y * y == 1
def init_Pell_eq_solution(D):

    p1, q1, p2, q2 = 0, 1, 1, 0

    for a in gen_cfraction_coeff(D):

        p2, p1 = a * p2 + p1, p2

        q2, q1 = a * q2 + q1, q2

        if is_satisfied_equation(p2, q2, D):

            return p2, q2
print init_Pell_eq_solution(1141)