前回の投稿でペル型方程式の解を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を解くのに役に立ちます。