@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余る素数なら解がある。その他の場合は、解があるかどうかわかりません。