[SymPy] ペル方程式

[SymPy] ペル方程式

2元2次不定方程式

 2元2次不定方程式
 
\[ax^2 + bxy + cy^2 + dx + ey + f = 0\]
は必ず $x^2-Dy^2=N$ という形に帰着されることが知られています。find_DN() を使うと、変換されたときの $D,\ N$ を得ることができます。また、diop_DN(D, N) によって基本解を求めることができます。

# 2元2次不定方程式

import sympy
from sympy.solvers.diophantine import find_DN, diop_DN

# 記号x,yを定義
sympy.var("x y")

# x**2-Dy = N の形式に変換されたときの(D,N)を取得
dn = find_DN(x**2 - 3*x*y + y**2 - 2*x + 6*y - 3)

# x**2-Dy = N の解
sol = diop_DN(dn[0], dn[1])

print("(D, N) = {}".format(dn))
print("基本解:{}".format(sol))
(D, N) = (5, 220)
基本解:[(-15, 1), (25, 9), (15, 1), (45, 19), (20, 6), (60, 26)]

 diophantine() を使って一般解を得ることもできますが、ずらずらと長い式が現れてコンピュータに大きな負荷がかかります。特に LaTeX 形式で表示させようとすると、そのままフリーズすることもあるので(私のモバイルPCはフリーズしました)、高性能なコンピュータ以外では実行しないでください。
 

ペル方程式

ペル方程式の漸化式

 $x^2-Dy^2=N$ において $N=1$ とおいた
 
\[x^2-Dy^2=1\]
という形の不定方程式をペル方程式とよびます。ペルの方程式については、漸化式で解を求める方法が知られているので、簡単に解説しておきます。この方程式には $x=\pm 1,\ y=0$ という自明な解があります。また、$D$ が平方数 $a^2$ であるときには、
 
\[(x+ay)(x-ay)=1\]
となるので、
 
\[\begin{cases}x+ay=1\\[6pt]
x-ay=1\end{cases}\]
 または
 
\[\begin{cases}x+ay=-1\\[6pt]
x-ay=-1\end{cases}\]
 すなわち、$x=\pm 1,\ y=0$ が唯一の解となります。
 したがって以降は $D$ が平方数でない場合を考えます。
 また、一組の整数解を $(x,\ y)$ とすると、
 
\[(-x,\ y),\ (x,\ -y),\ (-x,\ -y)\]
も解となるのは明らかなので、$x,\ y\gt 0$ の場合を考えます。いま1組の解 $(x_1,\ y_1)$ が得られたとして
 
\[(x_1+y_1\sqrt{D})^2=x_2+y_2\sqrt{D}\tag{A}\]
とおけば、
 
\[(x_1-y_1\sqrt{D})^2=x_2-y_2\sqrt{D}\tag{B}\]
が成り立ちます。(A) と (B) を辺々掛けて、$x_1^2-Dy_1^2=1$ を用いると
 
\[x_2^2-Dy_2^2=1\]
となって、$(x_2,\ y_2)$ もまた $x^2-Dy^2=1$ の解であることがわかります。こうして得られた2組の解 $(x_2,\ y_2),\ (x_2,\ y_2)$ を用いて
 
\[(x_1+y_1\sqrt{D})(x_2+y_2\sqrt{D})=x_3+y_3\sqrt{D}\]
とおくと、先ほどと同様の手順で $(x_3,\ y_3)$ も解であることがわかります。したがって一般に、$(x_n,\ y_n)$ と $(x_{n+1},\ y_{n+1})$ の間に
 
\[x_{n+1}+y_{n+1}\sqrt{D}=(x_1+y_1\sqrt{D})(x_n+y_n\sqrt{D})\]
という関係が成り立ちます。この式を整理すると、
 
\[\begin{align*}&x_{n+1}=x_1x_n+y_1y_n D\\[6pt]
&y_{n+1}=y_1x_n+x_1x_n\end{align*}\]
という漸化式が得られます。行列で表すと
 
\[\binom{x_{n+1}}{y_{n+1}}=\begin{pmatrix}
x_1 & y_1D\\y_1 & x_1\end{pmatrix}\binom{x_n}{y_n}\]
となります。
 

SimPyによるペル方程式の解の取得

 例として $x^2-2y=1$ の解を SimPy で求めてみます。
 diop_DN(D, N) を使えば1組の解を求めることができます。

# ペルの方程式

import sympy
from sympy.solvers.diophantine import diop_DN

# 記号x,yを定義
sympy.var("x y")

# x**2-2y = 1 の1組の解を取得
sol = diop_DN(2, 1)[0]

print(sol)
(3, 2)

 diop_DN() の戻り値はタプルを要素にもつリストの形式です。上のコードではインデックス 0 番を指定してタプルを取り出しています。

 sympy.Matrix() を使って5組の解を求めてみます。

# ペルの方程式

import sympy
from sympy.solvers.diophantine import diop_DN

# 記号x,yを定義
sympy.var("x y")

# x**2-2y = 1 の1組の解を取得
sol = diop_DN(2, 1)[0]

x1 = sol[0]
y1 = sol[1]

# 行列の定義
mt_a = sympy.Matrix([[x1, y1*2], [y1, x1]])
mt_sol = sympy.Matrix([[x1], [y1]])

# 5組の解を表示
for k in range(5):
    mt_sol = mt_a * mt_sol
    print((mt_sol[0], mt_sol[1]), end = ",")
(17, 12), (99, 70), (577, 408), (3363, 2378), (19601, 13860)