読者です 読者をやめる 読者になる 読者になる

ペル方程式の最小解

Project Eulerをずっとやってるんだけど,やっとProblem 67まで解いた.
他の人を見てると自分の 遅さにうんざりする.
でもがんばる(*'ω')b


で,今日とけたProblem 66について
404 Not Foundとか ペル方程式の解の列挙方法 - まめめもとか404 Not Foundを参考に解いた
Problem 64でやった連分数展開を使う
気をつけなきゃいけないのは,連分数展開したとき繰り返し部分の周期(参考:404 Not Found
周期が奇数のときと偶数のときで若干計算が違う
そして,ソース

from math import sqrt, floor

def is_square(n):
    sqrt_n = math.sqrt(n)
    if math.floor(sqrt_n) == sqrt_n:
        return True
    return False

def expand_square_root(n):
    if is_square(n):
        return
    a = a0 = int(floor(sqrt(n)))
    s, t = 0, 1
    while True:
        s = a * t - s
        t = (n - s * s) / t
        a = int(floor((a0 + s) / t))
        yield a
        if t == 1:
            break

def min_sol_of_pells(D):
    P = [1, int(floor(sqrt(D)))]
    Q = [0, 1]
    f = False
    while True:
        for q in expand_square_root(D):
            f = not f
            P.append(q * P[-1] + P[-2])
            Q.append(q * Q[-1] + Q[-2])
        if not f:
            break
    return P[-2], Q[-2]