拡張Pell方程式の解

@Mi_Sawaさんの指摘を受けて書き直しました。

拡張Pell方程式を

x2 - Dy2 = -1 (Dは平方数でない自然数

とします。Pell方程式

x2 - Dy2 = 1 (Dは平方数でない自然数

はどのDに対しても解を持つことが知られています。では、拡張Pell方程式はどうでしょう。

まず、Pell方程式の初期解を(x0, y0)とします(初期解の求め方はペル方程式(2) 参照)。このとき(a, b)を拡張Pell方程式の初期解とすると、

a2 - Db2 = -1

ここから、

(a + √Db)(a - √Db) = -1
(a + √Db)2(a - √Db)2 = 1
(a2 + Db2 + 2√Dab)(a2 + Db2 - 2√Dab) = 1

だから、Pell方程式の初期解は(a2 + Db2, 2a b)となります。よって、

a2 + Db2 = x0

これと上の、

a2 - Db2 = -1

より、

a2 = (x0 - 1) / 2
b2 = (x0 + 1) / 2D

が共に平方数であることが条件になります。
さて、具体的にどのようなDのときに拡張Pell方程式が解を持つでしょうか。100まででそのようなDを抽出してみました。

2 5 10 13 17 26 29 37 41 50 53 58 61 65 73 74 82 85 89 97

どうやら2つのガウス整数に分解できる整数のようです。しかし、よく見ると、52のような2を2つ含むものはありません。また、58がありません。4で割って3でない素数は全て含まれているようです。ところが1000まで調べてみると、

421 541 661

の3つだけが4で割って1余る素数であるにもかかわらず拡張Pell方程式の解が無いようです。
と思ったのですが、そうではないらしいです。実はD = 421のとき初期解が超巨大になります。

(3879474045914926879468217167061449, 189073995951839020880499780706260)

だから、平方数の判定がおかしくなるわけです。ここで、数値計算による誤差の手法が使えるわけです。これを使うと正しく平方数の判定ができます。

from itertools import *
from math import sqrt
from fractions import gcd
import time

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 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

# x^2 - Dy^2 = -1
def init_extended_Pell_eq_solutions(D):
    x0, y0 = init_Pell_eq_solution(D)
    # a^2 + Db^2 = x0
    # 2ab = y0
    # a^2 - Db^2 = -1
    # a^2 = (x0 - 1) / 2
    if x0 % 2 == 0 or (x0 + 1) % (2 * D) != 0:
        return None
    a2 = (x0 - 1) / 2
    a = int(sqrt(a2))
    if a * a != a2:
        return None
    b = int(sqrt((a2 + 1) / D))
    return a, b

20000まで4で割って1余る素数なら解があることを確認しました。30000だと解が大きすぎて浮動小数点数に変換できず落ちてしまいます。対策のしようはありますが。
それから2pは、pが8で割って5余る素数なら解があるそうです。これは500000まで確認しました。

Dに4で割って3余る素数が入っていると解はありません。例えば3が入っていると、

a2 ≡ -1 (mod 3)

となります。7や11でも同様に解がありません。
また、4の倍数の場合も解がありません。これも簡単ですね。

まとめると、Dが2を2つ以上または4で割って3余る素数を持っていると解はない。Dが4で割って1余る素数なら解がある。D = 2pとしてpが8で割って5余る素数なら解がある。その他の場合は、解があるかどうかわかりません。