あまり速くならない。なぜこのコードでこんなに時間がかかるのか。未だにPythonはどこがネックになるのかわからない。
この問題は、
xy2 = (x + 1)z2 + x(x + 1)
という形に帰着する。
s = y / (x + 1)
t = z / x
とおくと、
(x + 1)s2 = xt2 + 1
xを固定すると、ペル方程式に似た形になっている。
x = 1とすると、
2s2 = t2 + 1
これは、（拡張）ペル方程式になっている。
x = 2とすると、
3s2 = 2t2 + 1
これは、ペル方程式の解法と同じように、√3 - √2のべき乗を計算する。
(√3 - √2)2 = 5 - 2√6
(√3 - √2)3 = 9√3 - 11√2
3 • 92 = 2 • 112 + 1が確認できる。このように、奇数乗を計算していけばよい。
x = 3とすると、
4s2 = 3t2 + 1
(2 - √3)2 = 7 - 4√3
で、べき乗を計算していけばよい（この場合、(s, t) = (7/2, 4)が解になっているが、yに戻すときに整数になる）。
x = 4とすると、
5s2 = 4t2 + 1
(√5 - 2)2 = -4√5 + 9
(√5 - 2)3 = 17√5 - 38
だから、これも奇数乗を計算していく。
これで行けるかと思ったが、例えば、x = 48とすると、
49s2 = 48t2 + 1
(7 - 4√3)2 = 97 - 56√3
ではダメで、
(2 - √3)n
が解になる。一般に、
√x = k12s1
√(x+1) = k22s2
となる最大のs1,s2に対し、
s1a2 = s2b2 + 1
となる最小のa,bを見つけて、
(a√s1 - b√s2)n
を計算する。s1が1でなければnは奇数で、1ならその制限はない。
