【SymPy】ペル方程式の解を求める
二元二次不定方程式
\[ax^2 + bxy + cy^2 + dx + ey + f = 0\]
は必ず $x^2-Dy^2=N$ という形に帰着されることが知られています。find_DN() を使うと、変換されたときの $D,\ N$ を得ることができます。また、diop_DN(D, N) によって基本解を求めることができます。
# PELL'S EQUATION
# In[1]
import sympy
from sympy.solvers.diophantine.diophantine import find_DN, diop_DN
# 記号x,yを定義
sympy.var("x y")
# x**2-Dy**2=Nの形式に変換されたときの(D,N)を取得
dn = find_DN(x**2 - 3*x*y + y**2 - 2*x + 6*y - 3)
# x**2-Dy**2=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\]
という形の不定方程式を ペル方程式 (Pell’s equation) とよびます。ペル方程式については、漸化式で解を求める方法が知られているので、簡単に解説しておきます。この方程式には $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$ の場合を考えます。いま、ペル方程式の一組の解 $(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$ の解であることがわかります。こうして得られた二組の解 $(x_1,\ y_1),\ (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}\]
となります。
例として、ペル方程式 $x^2-2y^2=1$ の解を SimPy で求めてみます。
diop_DN(D, N) を使えば基本解を求めることができます。
# In[2]
# 記号x,yを定義
sympy.var("x y")
# ペル方程式x**2-2y**2=1の基本解を取得
sol = diop_DN(2, 1)[0]
print(sol)
# (3, 2)
diop_DN() の戻り値はタプルを要素にもつリストの形式です。上のコードではインデックス 0 番を指定してタプルを取り出しています。sympy.Matrix() を使ってペル方程式の 5 組の解を求めてみます。
# In[3]
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),
コメント
ChatGPTに「ペル方程式」について尋ねました
[Q] ペル方程式 とは何ですか？
[ChatGPT] ペル方程式とは、二次不定方程式の一種で、以下の形式を持つ方程式を指します。
x^2 – dy^2 = 1
ここで、dは正の整数で、平方数でないものとします。ペル方程式は、数学者オイラーが命名し、数学者ペルがその解法を発見したことからその名が付けられました。ペル方程式は、古代ギリシャの数学者ディオファントスが解法を見つけられなかったことで知られており、その解法は、ペル自身や数学者ラグランジュによって見出されました。ペル方程式は、現代の数学においても、その解法が研究される重要な分野の一つです。ペル方程式の解は、一般的には無限に存在しますが、その中でも最小正整数解を求めることが問題となります。また、ペル方程式には、連分数に関する興味深い性質が存在することが知られています。ペル方程式の解法には、連分数法やペルのアルゴリズムなどがあります。
下記は誤植と思われますので、ご確認ください。
PELL’S EQUATION In[1] プログラムで、
from sympy.solvers.diophantine import find_DN, diop_DN →
from sympy.solvers.diophantine.diophantine import find_DN, diop_DN
修正しました。ありがとうございます。
下記は誤植と思われますので、ご確認ください。
In[1] プログラムのコメントで、# x**2-Dy=1 → # x**2-Dy**2=1 ( 2 ヶ所 )
In[1] プログラムの下の文で、
(x_2, y_2), (x_2, y_2) → (x_1, y_1), (x_2, y_2)
x^2-2y = 1 → x^2-2y^2=1
In[2] プログラムのコメントで、# x**2-2y=1 → # x**2-2y**2=1
In[3] プログラムの実行結果で、(19601, 13860) → (19601, 13860),
修正しました。
ありがとうございます。m(_ _)m
x**2 – 3*x*y + y**2 – 2*x + 6*y – 3 = 0 の整数解を求めようとして、「Diophantine」のチュートリアルを参考にして、以下のような手順でやってみましたがうまくいきません。間違っているところをご指摘いただけるとありがたいです。
find_DN(x**2 – 3*x*y + y**2 – 2*x + 6*y – 3)
# (D, N) = (5, 220)
diop_DN(5, 220)
# [(-15, 1), (25, 9), (15, 1), (45, 19), (20, 6), (60, 26)]
A, B = transformation_to_DN(x**2 – 3*x*y + y**2 – 2*x + 6*y – 3)
print(A)
# Matrix([[1/10, -3/10], [0, -1/5]])
print(B)
# Matrix([[14/5], [6/5]])
(X, Y) として基本解の一つの (15, 1) を採用して、
[X Y]^T=A[x, y]^T + B
の連立方程式を解いて得られた (x, y) = (125, 1) を元の式に代入しても
x**2 – 3*x*y + y**2 – 2*x + 6*y – 3 = 15004
となり 0 にならない。Xの値を小さくする方向で調整すると、(x, y)=(4, 1)が一つの解であることは見つけました。
下から 2 行目の、X の値を小さくする方向で → x の値を125から小さくする方向で
に訂正します。
チュートリアルの記述の一部に誤植があります。(X, Y) と (x, y) の変換公式は
[X Y]^T=A[x, y]^T + B
ではなくて、
[x y]^T=A[X, Y]^T + B
です。以下のコードを実行して確認してみてください。
ご回答ありがとうございました。基本解のうち(X, Y)=(-15, 1), (15, 1), (20, 6), (60, 26) に対して、(x, y) が整数解をもつことがわかりました。
Diophantine の説明文で、 transformation_to_DN 関数の解説のほうでは、Transpose([x y]) = A * Transpose([X Y]) + B となっているのを見つけました。