ペル方程式の最小解を求めます。
x2 - Dy2 = 1
これを変形すると、
x / y = √(D + 1 / y2)
となり、x / y は√Dの近似になっています。x, yが大きくなればなるほど√Dに近づいていきます。
√Dを連分数で表したとき、途中で打ち切ったものの内でペル方程式を満たすものがあります。実はその最小がペル方程式の最小解にもなります。例えば√2を連分数で表すと[1;2,2,...]なので、
1 ×
1 + 1 / 2 = 3 / 2 ○
(x, y) = (3, 2)が最小解となります。
√Dを連分数で表す計算方法を考えましょう。通常の連分数計算を行います。連分数計算を途中で打ち切ると次の形になります。
(bn + cn√D) / dn
これを超えない整数は、
an+1 = [ (bn + [cn√D]) / dn ]
これを引いて逆数を取ると、
(-dn(bn - an+1dn) + dncn√D) / (an+12D - (bn - an+1dn)2)
このanを使って連分数を作って、ペル方程式を満たすかどうかを調べます。
本当はもっとうまい方法があるようですが、これでも十分でしょう。
Project Eulerに必要な用語集
http://d.hatena.ne.jp/inamori/20091225/p1
from math import sqrt
from fractions import gcddef 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 / divdef is_satisfied_equation(x, y, D):
return x * x - D * y * y == 1def 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, q2print init_Pell_eq_solution(1141)