前回の投稿でペル型方程式の解をSympyのdiop_DNで生成する方法について書きました。今回はこれを利用して以下のような2元2次不定方程式を解くプログラムを書いて行きます。
$\displaystyle \large \ \ \ ax^2+bxy+cy^2+dx+ey+f=0$
手順は以下のようになります。
- 2元2次不定方程式をペル型方程式に変換、および元の解への変換Matrixを取得
- 前回の投稿で作った関数pelEqを使ってペル型方程式の解を生成
- 1.で取得したA,Bで元の解に変換
例として以下の2元2次不定方程式を使います
$\displaystyle \large \ \ \ 2x^2+4x-y^2-y=0$
1. 2元2次不定方程式をペル型方程式に変換、および元の解への変換Matrixを取得
sympyのfind_DNを使うとあっけなく$X^2-DY^2=N$の$D,N$が求まります。
(D, N) = find_DN(2*x**2 + 4*x - y**2 -y)
print(D, N)
# (8, 56)
同時にtransformation_to_DNを使ってて$X^2-DY^2=N$の解から元の$x,y$に変換する$Matrix A, B$を取得します。
(A, B) = transformation_to_DN(2*x**2 + 4*x - y**2 -y)
print(A,B)
# Matrix([[1/8, 0], [0, 1/2]]) Matrix([[-1], [-1/2]])
2. 前回の投稿で作った関数pelEqを使ってペル型方程式の解を生成
ここは前回の投稿を参考にしてください。
3. 1.で取得したA,Bで元の解$x_n,y_n$に変換
1.で取得した$A,B$を使って$X_n,Yn$を元の2元2次不定方程式の解$x_n,y_n$に変換します。式は以下ようになります。
$\displaystyle \large \ \ \ \begin{pmatrix} x_n \\ y_n\end{pmatrix} = A \times \begin{pmatrix} X_n \\ Y_n\end{pmatrix}+ B$
(x, y) = A*sympy.Matrix([[abs(X)], [abs(Y)]])+B
これらを組み合わせて2元2次不定方程式を解くプログラムdiophQuadが以下となります。結果の最後の値が$0$になっているので正しい答えが生成されているのが確認できます。
import sympy
from sympy.solvers.diophantine.diophantine import find_DN, diop_DN, transformation_to_DN
from sympy.parsing.sympy_parser import parse_expr
def findDNAB(eqstr):
sympy.var("x y")
eq = parse_expr(eqstr)
(D, N) = find_DN(eq)
(A, B) = transformation_to_DN(eq)
return (D, N, A, B)
# solve quadratic equation with interger solutions
def diophQuad(eq):
(D, N, A, B) = findDNAB(eq)
ansPelEq = pelEq(D, N)
while True:
xyl = next(ansPelEq)
for (X,Y) in xyl:
(x, y) = A*sympy.Matrix([[abs(X)], [abs(Y)]])+B
yield (x,y)
ansDQ = diophQuad(" 2*x**2 + 4*x - y**2 -y")
for i in range(10):
(x,y) = next(ansDQ)
print(x, y, 2*x**2 + 4*x - y**2 -y)
#0 0 0
#0 0 0
#1 2 0
#3 5 0
#10 15 0
#22 32 0
#63 90 0
#133 189 0
#372 527 0
#780 1104 0
（開発環境：Google Colab）
このアルゴリズムはProject Euler Problem 321: Swapping Countersを解くのに役に立ちます。