【Project Euler】Problem 66: ディオファントス方程式

• 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。

問題：ディオファントス（ペル）方程式

原文 Problem 66:Diophantine equation

問題の要約：　以下の不定方程式で$D<10^3$で$x$の最小解の最大値を求めよ

$\ \ \ \ \ \displaystyle \large x^2-Dy^2=1$

原文ではディオファントス方程式となっていますが、その中のペル方程式を扱ってますので、以下ペル方程式と呼びます。

まず過去の投稿「ペル型方程式をsympyのdiop_DNで解く」でsympyのdiop_DNを紹介しましたが、それがそのまま使って$N=7$とすると$D=5$が求まります。ただこれも$N=1000$とすると42秒と1分ルールには収まるもののあまり速くないので。別解を考えます。

from sympy.solvers.diophantine.diophantine import diop_DN

N,mx = 7, 0
for D in range(2,N+1):
  (x,y) = diop_DN(D,1)[0]
  if x > mx: 
    mx, mxD = x, D
print(mxD, mx)
#5 9

ペル方程式の基本解を平方根の連分数展開から求める

以前の投稿【Project Euler】Problem 64: 平方根の連分数でも参照した文献「平方根の連分数とペル方程式 有澤健治著」の87ページにある以下の定理によりペル方程式の解を平方根の連分数展開を使って求めることが出来るということです。

定理 2a.
ペル方程式 $p^2-mq^2=1$
$n=\lfloor\sqrt{m}\rfloor $、$\sqrt{m}$の連分数展開を$[n,n_1,n_2,\dots,]$、$l$を連分数の循環部分の周期とする。
そこで$p_k/q_k=[n,n_1,n_2,\dots ,n_{kl-1}]$ $(k=1,2,\dots)$とすると
$p_k^2-mq_k^2=(-1)^{kl}$である

従ってまず$D$の連分数展開を求めて循環部分の周期が偶数ならば$(-1)^{l}=1$となるので、$k=1$として

• $p_1/q_1=[n,n_1,n_2,\dots ,n_{l-1}]$

となり$p_1, q_1$がペル方程式の解となります。

循環部分の周期が奇数ならば、$(-1)^{kl}=1$となる最小の$k$は$2$なので

• $p_2/q_2=[n,n_1,n_2,\dots ,n_{2l-1}]$

の$p_2, q_2$が解となります。

$D$の連分数展開は【Project Euler】Problem 64: 平方根の連分数で作った関数cfsqrtを使い、連分数から分数に変換する関数をcont2fracとすると、ペル方程式の最小解を求める関数pellEqは以下のようになります。これをdiop_DNの代わりに使うと$N=1000$でも1秒以内に求めることが出来ました。

# Solve Pell Equation
def pellEq(D):
  cf = cfsqrt(D) # continued fraction of square root of D (Problem 64)
  if len(cf) == 1: return (1,0)  ## square numbers
  cf1 = [cf[0]]+cf[1]+cf[1][:-1] if len(cf[1])%2 == 1 else [cf[0]]+cf[1][:-1]
  return cont2frac(cf1)

# Convert continued fraction to fraction
def cont2frac(cont):
    num, den = 1, 0
    for u in reversed(cont):
      num, den = den + num*u, num
    return num, den