拡張Pell方程式の解（2）

Pell方程式の初期解は、連分数から求められます（ペル方程式（2）参照）。拡張Pell方程式の初期解があるときは、この連分数がPell方程式の初期解が求まる前に求められます。D = 5でやってみましょう。

√5 = 2 + (√5 - 2) = 2 + 1/(√5 + 2) = 2 + 1/(4 + (√5 - 2)) = 2 + 1/(4 + 1/(√5 + 2))

だから、√5 = [2; 4, 4, ...]と表されます。このとき途中で分数を打ち切っていくと、

2/1, 9/4, 38/17, ...

これを

(x, y) = (2, 1), (9, 4), (38, 17), ...

と解釈して、x² - Dy²を計算すると、

-1, 1, -1

で、(9, 4)がPell方程式の初期解になりますが、拡張Pell方程式の初期解は(2, 1)になります。
次にD = 3を見ましょう。

√3 = 1 + (√3 - 1) = 1 + 1/((√3 + 1)/2) = 1 + 1/(1 + (√3 - 1)/2)

だから、分数を途中で打ち切ると、

1, 2, ...

これを、

(x, y) = (1, 1), (2, 1), ...

と解釈して、x² - Dy²を計算すると、

-2, 1, ...

で(2, 1)がPell方程式の初期解となり、拡張Pell方程式の解はありません。
この考え方に基づいてPythonでコードを組むと、平方根の計算が配慮無くできるようになり、3倍以上速くなりました。

def int_sqrt(n):
    return int(sqrt(n))

def is_square(n):
    m = int_sqrt(n)
    return m * m == 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 init_extended_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
        s = p2 * p2 - D * q2 * q2
        if s == -1:
            return p2, q2
        elif s == 1:
            return None