Project Euler 進歩状況

今日は問66を解きました。
Problem 66 - PukiWiki
問題文に出てくる  X^2 - D Y^2 = 1 は、ペル方程式と言うんだそおです。

最初は X = \sqrt{1 + D Y^2}とおいて、X が整数になるような、最小のYを探すという方針でプログラムを書きましたが時間がかかりすぎるので途中で諦めました。

ググって見ると,Dの大きさに比べて(X,Y)は非常に大きくなり得るため、上のような方法では時間がかかりすぎるとのこと。

解法としては連分数を用いる方法が一般的のようです。
 \sqrt{D} の連分数による近似を、分子=X, 分母=Y と置いて上の方定式を満たすまで続けます。
正直なんでこの方法で解が見つかるのかは理解できていません。(-.-;)
機会があれば勉強して見ようと思います。

今日までで、85問解きました。

以下 python によるコードです。

from math import floor, sqrt
def pellequation(num):
    if num**0.5 == int(num**0.5):
        return 0, 0
    flag, b, c = False,0, 1
    P1, P2 = 1, int(sqrt(num))
    Q1, Q2 = 0, 1
    while True:
        i = int( (sqrt(num)+b) / c)
        if P2**2 - Q2**2*num == 1:
            return P2, Q2
        if flag:
            P1, P2 = P2, i*P2 + P1
            Q1, Q2 = Q2, i*Q2 + Q1
        flag = True
        b  = i*c - b
        c = (num - b**2) / c
 
ans, max_x = 0, 0
for i in range(1, 1001):
    tmp = pellequation(i)
    if max_x < tmp[0]:
        max_x = tmp[0]
        ans = i
 
print ans