ペル型方程式をsympyのdiop_DNで解く

ペル型方程式

以下の形の2次不定方程式をペル方程式(Pell's equation)と呼びますが、
$\displaystyle \large \ \ \ \ x^2-Dy^2=1$
今回は一歩進んで右辺が$N$のペル型方程式(Pell-like equation)をの解を生成するプログラムをpythonのsympyを使って書きました。

$\displaystyle \large \ \ \ \ x^2-Dy^2=N$ 手順は以下の通りです。

1. diop_DNで基本解を求める
2. $N=1$の最小解を漸化式の係数とする
3. 係数Matrixを基本解に掛けて解を生成していく

1. diop_DNで基本解を求める

以下を例題を解きます。diop_DNに$D=8, N=56$を与えて呼ぶだけですね。diop_DNの説明はこちら(英文)
$\displaystyle \large \ \ \ x^2-8y^2=56$

import sympy
from sympy.solvers.diophantine.diophantine import diop_DN
D, N = 8, 56
xylist = diop_DN(D,N)
print(xylist)
# [(-8, 1), (8, 1)]

簡単に2つの基本解$(-8, 1),\ (8, 1)$が求まります。

2. N=1の最小解を漸化式の係数とする

同様にN=1にして最小解$(3,1)$が求まります。

D, N = 8, 1
(u,v) = diop_DN(D,N)[0]
print((u,v))
# (3, 1)

ここから他の整数解は英語版WikipediaのGeneralized Pell's equationに書いてあるように以下の漸化式で生成できます。

$\displaystyle \large \ \ \ \ x_{n+1}+y_{n+1}\sqrt d = (x_n+y_n\sqrt d)(u+v\sqrt d)$

行列で表すと

$\displaystyle \begin{pmatrix} x_{n+1} \\ y_{n+1}\end{pmatrix} =\begin{pmatrix} u & vD \\ v & u\end{pmatrix}\begin{pmatrix} x_{n} \\ y_{n}\end{pmatrix} $

となるのでこの行列を掛けていけば無限に解を生成することができます。上の例では行列をAとすると以下のようになります。

$A =\begin{pmatrix} 3 & 1 \times 8 \\\ 1 & 3\end{pmatrix}$

3. 係数Matrixを基本解に掛けて解を生成していく

例えば基本解の(8,1)に$A$を掛けると$(32,11)$となり、やはり$x^2-8y^2=56$の解となることが確認できます。

$\displaystyle \ \ \ \ \begin{pmatrix} 3 & 8 \\ 1 & 3\end{pmatrix} \times \begin{pmatrix} 8 \\ 1\end{pmatrix} = \begin{pmatrix} 32 \\ 11\end{pmatrix} $

これを繰り返して解を生成していくプログラムは以下のようになります。yieldを使っているのでnextを使って好きなだけ解を生成して行けます。

def pelEq(D,N):
  # (u,v): solution to x**2-D*y**2=1
  (u,v) = diop_DN(D, 1)[0]  
  # Matrix to generate other solution
  A = sympy.Matrix([[u, v*D], [v, u]]) 
  # List of fundamental solutions 
  xyl = [sympy.Matrix([[x1], [y1]]) for (x1, y1) in diop_DN(D, N)]

  while True:
    yield [(xy[0],xy[1]) for xy in xyl]
    for i in range(len(xyl)):
      xyl[i] = A * xyl[i]

D, N = 8, 56
ansPelEq = pelEq(D, N)
for i in range(10):
  xylist = next(ansPelEq)
  print(xylist)
# [(-8, 1), (8, 1)]
# [(-16, -5), (32, 11)]
# [(-88, -31), (184, 65)]
# [(-512, -181), (1072, 379)]
# [(-2984, -1055), (6248, 2209)]
# [(-17392, -6149), (36416, 12875)]
# [(-101368, -35839), (212248, 75041)]
# [(-590816, -208885), (1237072, 437371)]
# [(-3443528, -1217471), (7210184, 2549185)]
# [(-20070352, -7095941), (42024032, 14857739)]

結果は、最初の呼び出しで基本解を、後はその基本解の各々に係数Matrixを掛けたものになります。値がマイナスの物もありますが式の性質状、絶対値を取った物も解です。基本解がたくさんある場合などは生成された解は必ずしも値の小さい順になっているとは限らないので必要に応じてソートする必要があります。