お気楽 Common Lisp プログラミング入門

分数 [2]

●連分数

連分数 (continued fraction) は、分数の分母にさらに分数が含まれる形をした分数のことです。

\( b_0 + \dfrac{c_1}{b_1 + \dfrac{c_2}{b_2 + \dfrac{c_3}{b_3 + \dfrac{c_4}{b_4 + \cdots}}}} \)

連分数の中で、分子 \(c_1, c_2, \ldots\) がすべて 1 で、\(b_0\) が整数、\(b_1, b_2, \ldots\) がすべて正の整数であるような連分数を「正則連分数」といいます。正則連分数は \(b_i\) の情報を保持するだけで、連分数に復元することができます。そこで、以下の右辺のように記述することが多いようです。

\( b_0 + \dfrac{1}{b_1 + \dfrac{1}{b_2 + \dfrac{1}{b_3 + \dfrac{1}{b_4 + \cdots}}}} = [b_0; \ b_1, \ b_2, \ \ldots ] \)

●有理数の連分数展開

有理数 (分数) は「ユークリッドの互除法」を使って正則連分数に展開することができます。

\(a, b\) の割り算を \(a = q \times b + r\) で表す

\(\begin{eqnarray} \dfrac{a}{b} &=& q + \dfrac{r}{b} \\ &=& q + \dfrac{1}{\frac{b}{r}} \end{eqnarray}\)
\(\dfrac{a}{b}\) の正則連分数展開するには \(\dfrac{b}{r}\) を正則連分数展開すればよい

\(a, b\) の問題を \(b, r\) の問題に帰着させるのは「ユークリッドの互除法」と同じです。この計算を繰り返し行うと、\(a, b\) はどんどん小さくなっていき、\(r = 0\) になったとき連分数展開が終了します。

簡単な例を示します。

\(\begin{eqnarray} \dfrac{17}{11} &=& 1 + \dfrac{6}{11} = 1 + \dfrac{1}{ \frac{11}{6} } \\ &=& 1 + \dfrac{1}{1 + \dfrac{5}{6} } = 1 + \dfrac{1}{ 1 + \dfrac{1}{\frac{6}{5}} } \\ &=& 1 + \dfrac{1}{1 + \dfrac{1}{1 + \dfrac{1}{5}} } = 1 + \dfrac{1}{ 1 + \dfrac{1}{1 + \dfrac{1}{\frac{5}{1}}} } \\ &=& 1 + \dfrac{1}{ 1 + \dfrac{1}{ 1 + \dfrac{1}{ 5 + \dfrac{0}{1} } } } = 1 + \dfrac{1}{ 1 + \dfrac{1}{ 1 + \dfrac{1}{5} } } \\ \dfrac{355}{113} &=& 3 + \dfrac{16}{113} = 3 + \dfrac{1}{ \frac{113}{16} } \\ &=& 3 + \dfrac{1}{7 + \frac{1}{16}} = 3 + \dfrac{1}{ 7 + \dfrac{1}{\frac{16}{1}} } \\ &=& 3 + \dfrac{1}{ 7 + \dfrac{1}{16 + \frac{0}{16}} } = 3 + \dfrac{1}{ 7 + \dfrac{1}{16} } \\ \end{eqnarray}\)

SBCL で検算してみましょう。

* (+ 1 (/ (+ 1 (/ (+ 1 (/ 5))))))
17/11

* (+ 3 (/ (+ 7 (/ 16))))
355/113

連分数に正しく展開されていますね。

●プログラムの作成

プログラムは簡単です。次のリストを見てください。

リスト : 有理数を正則連分数に展開する

(defun cont-frac (x)
  (unless (rationalp x)
    (error "must be rational"))
  (do ((a (numerator x))
       (b (denominator x))
       (xs nil))
      ((<= b 0) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor a b)
       (push p xs)
       (setf a b b q))))

関数 cont-frac は引数 x を正則連分数展開します。最初に引数 x が有理数か rationalp でチェックします。そうであれば、分子を変数 a に、分母を変数 b にセットします。展開した結果はリスト xs にセットします。あとは、ユーグリッドの互除法と同じです。b が 0 より大きければ、floor で a と b の商と剰余を求めて変数 p, q にセットします。そして、p を xs に追加し、a を b に、b を q に更新します。最後に、xs を nreverse で反転して返します。

リスト : 正則連分数を有理数に戻す

;;; 単純な再帰呼び出し
(defun cont-frac-reduce (xs)
  (if (null (cdr xs))
      (car xs)
    (+ (car xs)
       (/ (cont-frac-reduce (cdr xs))))))

;;; 畳み込み (reduce) を使う
(defun cont-frac-reduce1 (xs)
  (reduce (lambda (x a) (+ x (/ a))) xs :from-end t))

関数 cont-frac-reduce は引数の正則連分数を通常の分数に直します。アルゴリズムは簡単です。リスト xs の要素がひとつであれば、その要素を返します。そうでなければ、cont-frac-reduce を再帰呼び出しして、その返り値の逆数と xs の先頭要素の足し算を返します。xs の先頭要素を \(x\) とし、cont-frac-reduce の返り値を \(a\) とすると、これで連分数 \(x + \frac{1}{a}\) を求めることができます。あとは、xs の後ろから順番に連分数を求めていくだけです。cont-frac-reduce1 は畳み込み reduce を使ったバージョンです。

もう一つ、簡単な方法を紹介しましょう。参考 URL 4 によると、連分数 \([a_0; a_1, a_2, \ldots, a_{n-2}, a_{n-1}]\) を約分した有理数 \(\frac{p_n}{q_n}\) の分子と分母 \(p_n, q_n\) は、以下に示す漸化式で求めることができるそうです。

整数列 \([a_0, a_1, a_2, \ldots, a_{n-2}, a_{n-1}]\) において、以下の漸化式を定義する

\( \begin{cases} p_0 = 1 & \\ p_1 = a_0 & \\ p_n = a_{n-1} p_{n-1} + p_{n-2} \quad & (\mathrm{if} \ n \geqq 2) \end{cases} \quad \begin{cases} q_0 = 0 & \\ q_1 = 1 & \\ q_n = a_{n-1} q_{n-1} + q_{n-2} \quad & (\mathrm{if} \ n \geqq 2) \end{cases} \)

整数列が正則連分数展開だとすると、その有理数は次式で求めることができる

\( [a_0; a_1, a_2, \ldots, a_{n-2}, a_{n-1}] = \dfrac{p_n}{q_n} = \dfrac{a_{n-1} p_{n-1} + p_{n-2}}{a_{n-1} q_{n-1} + q_{n-2}} \)

Common Lisp でプログラムすると、次のようになります。

リスト : 連分数を有理数に約分する (2)

;;; 漸化式を使う
(defun cont-frac-sub (xs a0 a1)
  (cdr (reduce (lambda (a x) (cons (cdr a) (+ (* x (cdr a)) (car a))))
               (cdr xs)
               :initial-value (cons a0 a1))))

(defun cont-frac-reduce2 (xs)
  (/ (cont-frac-sub xs 1 (car xs)) (cont-frac-sub xs 0 1)))

関数 cont-frac-sub は漸化式を計算します。これも畳み込み reduce を使うと簡単です。関数 cont-frac-reduce2 は cont-frac-sub を 2 回呼び出して分子 p と分母 q を求めて、p / q を計算して返すだけです。

簡単な実行例を示します。

* (cont-frac 17/11)
(1 1 1 5)
* (cont-frac -17/11)
(-2 2 5)
* (cont-frac 355/113)
(3 7 16)

* (cont-frac-reduce '(1 1 1 5))
17/11
* (cont-frac-reduce '(-2 2 5))
-17/11
* (cont-frac-reduce '(3 7 16))
355/113

* (cont-frac-reduce1 '(1 1 1 5))
17/11
* (cont-frac-reduce2 '(1 1 1 5))
17/11
* (cont-frac-reduce1 '(-2 2 5))
-17/11
* (cont-frac-reduce2 '(-2 2 5))
-17/11
* (cont-frac-reduce1 '(3 7 16))
355/113
* (cont-frac-reduce2 '(3 7 16))
355/113

* (cont-frac-reduce '(1 2 3 4 5 6 7 8))
81201/56660
* (cont-frac 81201/56660)
(1 2 3 4 5 6 7 8)

* (cont-frac-reduce1 '(1 2 3 4 5 6 7 8))
81201/56660
* (cont-frac-reduce2 '(1 2 3 4 5 6 7 8))
81201/56660

正常に動作していますね。

●無理数の連分数展開

有理数の連分数展開は有限回で終了しますが、無理数の連分数展開は無限に続きます。たとえば、\(\sqrt 2\) の連分数展開を求めてみましょう。

\(x = \sqrt 2\) とすると
\(\begin{array}{l} x^2 = 2 \\ x^2 - 1 = 1 \\ (x + 1)(x - 1) = 1 \\ x - 1 = \dfrac{1}{x + 1} \\ x = 1 + \dfrac{1}{1 + x} \end{array}\)
となる。右辺の \(x\) に \(1 + \dfrac{1}{1 + x}\) を代入していくと、次式のように連分数展開される

\(\begin{eqnarray} x &=& 1 + \dfrac{1}{1 + x} \\ &=& 1 + \dfrac{1}{2 + \frac{1}{1 + x}} \\ &=& 1 + \dfrac{1}{2 + \dfrac{1}{2 + \frac{1}{1 + x}}} \\ &=& 1 + \dfrac{1}{2 + \dfrac{1}{2 + \dfrac{1}{2 + \frac{1}{\cdots}}}} = [1; 2, 2, 2, \ldots] \end{eqnarray}\)

係数が有理数の二次方程式の根を「二次無理数」といいます。二次無理数の正則連分数展開は循環することが知られています。

\(\begin{array}{l} \sqrt 3 = [1; 1, 2, 1, 2, 1, 2, 1, 2, \ldots] \\ \sqrt 5 = [2; 4, 4, 4, 4, 4, 4, 4, 4, \ldots] \\ \sqrt 7 = [2; 1, 1, 1, 4, 1, 1, 1, 4, \ldots] \\ \end{array}\)

実際に試してみましょう。

* (sqrt 2d0)
1.4142135623730951d0
* (cont-frac-reduce '(1 2 2 2 2 2 2 2 2))
1393/985
* (float 1393/985 1d0)
1.4142131979695431d0

* (sqrt 3d0)
1.7320508075688772d0
* (cont-frac-reduce '(1 1 2 1 2 1 2 1 2))
265/153
* (float 265/153 1d0)
1.7320261437908497d0

* (sqrt 5d0)
2.23606797749979d0
* (cont-frac-reduce '(2 4 4 4 4 4 4 4 4))
219602/98209
* (float 219602/98209 1d0)
2.236067977476606d0

* (sqrt 7d0)
2.6457513110645907d0
* (cont-frac-reduce '(2 1 1 1 4 1 1 1 4))
590/223
* (float 590/223 1d0)
2.6457399103139014d0

●浮動小数点数の連分数展開

浮動小数点数 \(\alpha\) の正則連分数展開は次のアルゴリズムで求めることができます。

\(\alpha\) を超えない最大の整数を \(b_0\) とおくと
\( \alpha = b_0 + (\alpha - b_0) = b_0 + \dfrac{1}{\alpha_1} \)
ここで \(\alpha_1 = \dfrac{1}{\alpha - b_0}\)

\(\alpha_1\) に対して同様なことを行うと

\( \alpha = b_0 + \dfrac{1}{b_1 + \frac{1}{\alpha_2}}, \quad \alpha_2 = \dfrac{1}{\alpha_1 - b1} \)

この操作を n 回繰り返すと、\(\alpha\) を n 段の正則連分数に展開できる

\( \alpha = b_0 + \dfrac{1}{b_1 + \dfrac{1}{b_2 + \dfrac{1}{b_3 + \dfrac{1}{ \ddots b_{n-1} + \dfrac{1}{\alpha_n} } }}} \)

プログラムは次のようになります。

リスト : 浮動小数点数を正則連分数に展開する

(defun cont-frac-flt (x &key (k 8) (e 1e-16))
  (do ((k k (1- k))
       (xs nil))
      ((zerop k) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor x)
       (push p xs)
       (when (< (abs q) e)
         (return (nreverse xs)))
       (setf x (/ q)))))

引数 x が連分数展開を行う浮動小数点数、k が段数、e が 0 と判定する基準値です。浮動小数点数は不正確な数なので、演算誤差が発生します。k の値を増やしても、規則性は無限に続かずに途中で失われることがあります。

また、0.5 や 1.25 など浮動小数点数でも正確に表現できる場合、途中で q の値は 0 になるので、そこで展開を終了します。ただし、演算誤差が生じる場合があるので、ぴったり 0 になるとは限りません。そこで、\(|q| \lt e\) よりも小さくなったならば、0 と判定して展開を終了します。

それでは実際に試してみましょう。

* (cont-frac-flt (sqrt 2d0))
(1 2 2 2 2 2 2 2)
* (cont-frac-flt (sqrt 3d0))
(1 1 2 1 2 1 2 1)
* (cont-frac-flt (sqrt 5d0))
(2 4 4 4 4 4 4 4)
* (cont-frac-flt (sqrt 7d0))
(2 1 1 1 4 1 1 1)
* (cont-frac-flt (sqrt 7d0) :k 16)
(2 1 1 1 4 1 1 1 4 1 1 1 4 1 1 1)

* (cont-frac-flt 0.5d0)
(0 2)
* (cont-frac-flt 0.25d0)
(0 4)
* (cont-frac-flt 0.125d0)
(0 8)
* (cont-frac-flt 1.1d0)
(1 9 1 112589990684262 2 2 6 46912496118442)
* (cont-frac-flt 1.1d0 :e 1e-14)
(1 9 1)
* (cont-frac-reduce '(1 9 1))
11/10

正常に動作しているようです。興味のある方はいろいろ試してみてください。

●平方根の連分数展開

正の整数 m の平方根は二次無理数なので連分数は循環します。ただし、循環節が長くなると浮動小数点数の演算誤差により、上記プログラムでは正確に連分数展開できない場合があります。参考 URL 2 には、平方根の連分数展開 (循環節) を正確に求めるプログラム (Python) が掲載されています。ここで、簡単にアルゴリズムを紹介しましょう。

二次無理数 \(\alpha\) を次式のように連分数展開するとします。

\( \alpha_0 = c_0 + \dfrac{1}{\alpha_1}, \alpha_1 = c_1 + \dfrac{1}{\alpha_2}, \ldots, \alpha_n = c_n + \dfrac{1}{\alpha_{i+1}}, \ldots\ \)

このとき、\(\alpha_n\) を整数 \(a_n, b_n\) を使って以下のように表すと、\(\alpha\) の連分数展開を \(a_n, b_n\) の漸化式で表すことができます。

\(\begin{array}{l} a_0 = 1 \\ b_{-1} = 0 \\ \alpha_0 = \dfrac{\sqrt m + b_{-1}}{a_0} = \sqrt m, \, c_0 = \lfloor \alpha_0 \rfloor = \lfloor \sqrt m \rfloor\\ \alpha_n = \dfrac{\sqrt m + b_{n-1}}{a_n}, \, c_n = \lfloor \alpha_n \rfloor \\ \end{array}\)

\(a_n, b_n\) の漸化式は以下のように求めることができます。

\(\begin{eqnarray} \alpha_n &=& \dfrac{\sqrt m + b_{n-1}}{a_n} = c_n + \dfrac{1}{\alpha_{n+1}} \\ \dfrac{1}{\alpha_{n+1}} &=& \dfrac{\sqrt m + b_{n-1}}{a_n} - c_n = \dfrac{\sqrt m - (c_n a_n - b_{n-1})}{a_n} \\ \alpha_{n+1} &=& \dfrac{a_n}{\sqrt m - (c_n a_n - b_{n-1})} = \dfrac{\sqrt m + (c_n a_n - b_{n-1})}{\frac{m - (c_n a_n - b_{n-1})^2}{a_n}} = \dfrac{\sqrt m + b_n}{a_{n+1}} \\ \end{eqnarray}\)

\( \sqrt m + (c_n a_n - b_{n-1}) = \sqrt m + b_n \quad \therefore \ b_n = c_n a_n - b_{n-1} \)

\( \dfrac{m - (c_n a_n - b_{n-1})^2}{a_n} = a_{n+1} \quad \therefore\ a_{n+1} = \dfrac{m - {b_n}^2}{a_n}, \ c_{n+1} = \left\lfloor \dfrac{\sqrt m + b_n}{a_{n+1}} \right\rfloor \)

参考 URL 2 によると、\(c_{n+1}\) を求める計算は、\( \left\lfloor \frac{c_0 + b_n}{a_{n+1}} \right\rfloor \) に置き換えることができるそうです。証明もなされているので、興味のある方はお読みくださいませ。

このアルゴリズムを Common Lisp で実装すると、次のようになります。

リスト : 平方根の連分数展開

(defun cont-frac-root (m)
  (let* ((c0 (isqrt m))
         (cf (list c0)))
    (if (= (* c0 c0) m)
        cf
      (let ((a 1) (b 0) (c c0))
        (loop
         (setf b (- (* c a) b)
               a (floor (- m (* b b)) a)
               c (floor (+ c0 b) a))
         (push c cf)
         (when (= a 1)
           (return (nreverse cf))))))))

連分数展開の結果はリスト cf に格納します。漸化式の変数 \(a_n, \ b_n, \ c_n\) を変数 a, b, c で表しています。最初に (isqrt m) で m の平方根の整数部を求めて変数 c0 にセットし、cf を (list c0) に初期化します。c0 の二乗が m と等しい場合、m は平方数なので cf をそのまま返します。

m が平方数でなければ、\(\sqrt m\) を漸化式で連分数に展開します。変数 a, b, c を初期化して、loop で漸化式を計算します。これはアルゴリズムをそのままプログラムするだけです。a の値が 1 に戻ったら循環したので、nreverse で cf を反転して return で loop を脱出します。

簡単な実行例を示します。

* (dotimes (m 16) (format t "~D: ~S~%" (1+ m) (cont-frac-root (1+ m))))
1: (1)
2: (1 2)
3: (1 1 2)
4: (2)
5: (2 4)
6: (2 2 4)
7: (2 1 1 1 4)
8: (2 1 4)
9: (3)
10: (3 6)
11: (3 3 6)
12: (3 2 6)
13: (3 1 1 1 1 6)
14: (3 1 2 1 6)
15: (3 1 6)
16: (4)
NIL

* (cont-frac-root 97)
(9 1 5 1 1 1 1 1 1 5 1 18)
* (cont-frac-flt (sqrt 97d0) :k 24)
(9 1 5 1 1 1 1 1 1 5 1 18 1 5 1 1 1 1 1 1 5 1 11 2)

* (cont-frac-root 1234)
(35 7 1 3 1 4 4 2 9 1 1 2 3 1 1 34 1 1 3 2 1 1 9 2 4 4 1 3 1 7 70)
* (cont-frac-flt (sqrt 1234d0) :k 32)
(35 7 1 3 1 4 4 2 9 1 1 2 3 1 1 33 1 1 4 1 1 1 2 1 11 10 4 1 4 18 1 1)

cont-frac-flt では循環する前に演算誤差により規則性が失われることがありますが、cont-frac-root は循環節を正確に求めることができます。

●無理数の有理数近似

無理数を有理数で近似することを考えます。たとえば、\(\sqrt 2 = 1.4142135623730951\) の近似値として、以下のような有理数を考えることができます。

\( \sqrt 2 \fallingdotseq \dfrac{14}{10}, \ \dfrac{141}{100}, \ \dfrac{1414}{1000}, \ \dfrac{14142}{10000}, \ \ldots \)

これらの近似値の誤差を Common Lisp で求めると、次のようになります。

* (dolist (q '(10 100 1000 10000))
  (let ((r (/ (floor (* (sqrt 2d0) q)) q)))
  (format t "~G, ~G, ~G~%" (abs (- (sqrt 2d0) r)) (float (/ q)) (float (/ (* q q))))))
1.421356237309523400d-2, 0.1    , 1.00e-2
4.213562373095225400d-3, 1.00e-2, 1.0000e-4
2.13562373095221860000d-4, 1.000e-3, 1.000000e-6
1.356237309524388500000d-5, 1.0000e-4, 1.00000000e-8
NIL

近似値を \(r = \frac{p}{q}\) とすると、誤差は \(\frac{1}{q}\) の範囲に収まりますが、\(\frac{1}{q^2}\) には収まっていません。実数 \(\alpha\) の近似値 \(\frac{p}{q}\) が次の不等式を満たすとき、有理数 \(\frac{p}{q}\) を「\(\alpha\) の良い近似」といいます。

\( \left|\alpha - \dfrac{p}{q}\right| \lt \dfrac{1}{q^2} \)

誤差を \(\frac{1}{q^2}\) の範囲内に収めるのは難しいように思いますが、「\(\alpha\) が無理数なら \(\alpha\) の良い近似が無限個ある」ことが証明されています。これを「ディリクレのディオファントス近似定理」といいます。無理数 \(\alpha\) の正則連分数展開を途中で打ち切ったもの \(\frac{p_n}{q_n} = [b_0, b_1, \ldots, b_n]\) は、無理数 \(\alpha\) の良い近似になることが知られていて、その誤差は以下の式で表されます。

\( \left|\alpha - \dfrac{p_n}{q_n}\right| \lt \dfrac{1}{{q_n}^2} \)

\(\sqrt 5\) で試してみると、次のようになります。

* (cont-frac-reduce '(2 4))
9/4
* (abs (- (sqrt 5d0) 9/4))
0.013932022500210195d0
* (float 1/16 1d0)
0.0625d0

* (cont-frac-reduce '(2 4 4))
38/17
* (abs (- (sqrt 5d0) 38/17))
7.738598527309293d-4
* (float (/ (* 17 17)) 1d0)
0.0034602076124567475d0

* (cont-frac-reduce '(2 4 4 4))
161/72
* (abs (- (sqrt 5d0) 161/72))
4.31336113213554d-5
* (float (/ (* 72 72)) 1d0)
1.9290123456790122d-4

誤差は \(\frac{1}{{q_n}^2}\) よりも小さくなることがわかります。

●ペル方程式

変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 \(ax + by = c \ (c \ne 0)\) は、\(a, b, c\) が整数で、\(c\) が \(a\) と \(b\) の最大公約数の倍数であるとき、整数解を持つことが知られています。これを「ベズーの等式」といいます。二次の不定方程式では、以下に示す「ペル方程式 (Pell's Equation)」が有名です。

\( x^2 - m y^2 = 1 \)

\(m\) が平方数でない正の整数とき、ペル方程式は自明の解 (\(x = 1, y = 0\)) 以外の整数解を無限個持つことが知られています。その中で \(x + y\sqrt m\) を最小にする解を \(x_1, y_1\) とすると、ペル方程式の解は次式で表すことができます。

\( x_n + y_n \sqrt m = \left(x_1 + y_1 \sqrt m\right)^n \)

ただし、\(n\) は正の整数

ペル方程式の解はこれで全部です。最小解 \(x_1, y_1\) は \(\sqrt m\) の連分数展開で求めることができます。

\(\sqrt m\) の 1 周期分の連分数展開を \([a_0; a_1, a_2, \ldots, a_{k-1}, a_k]\) とする
その有理数近似 \([a_0; a_1, a_2, \ldots, a_{k-1}] =\frac{p}{q}\) は次式を満たす最小解になる

\( p^2 -mq^2 = \pm1 \)

右辺が -1 になる場合は \(p + q\sqrt m\) の二乗を計算する
\(\begin{array}{l} x + y\sqrt m = \left(p + q\sqrt m\right)^2 = p^2 + mq^2 + 2pq\sqrt m \\ x = p^2 + mq^2 \\ y = 2pq \end{array}\)

簡単な例を示しましょう。

\(\begin{array}{l} \sqrt 2 = [1; 2], [1] = \dfrac{1}{1} \\ 1^2 - 2 \times 1^2 = -1 \\ x = 1^2 + 2 \times 1^2 = 3 \\ y = 2 \times 1 \times 1 = 2 \\ 3^2 - 2 \times 2^2 = 1 \end{array}\)

\(\begin{array}{l} \sqrt 3 = [1; 1, 2], [1; 1] = \dfrac{2}{1} \\ 2^2 - 3 \times 1^2 = 1 \\ x = 2 \\ y = 1 \end{array}\)

\(\begin{array}{l} \sqrt 7 = [2; 1, 1, 1, 4], [2; 1, 1, 1] = \dfrac{8}{3} \\ 8^2 - 7 \times 3^2 = 64 - 63 = 1 \\ x = 8 \\ y = 3 \end{array}\)

\(\begin{array}{l} \sqrt 13 = [3; 1, 1, 1, 1, 6], [3; 1, 1, 1, 1] = \dfrac{18}{5} \\ 18^2 - 13 \times 5^2 = 324 - 325 = -1 \\ x = 18^2 + 13 \times 5^2 = 649 \\ y = 13 \times 18 \times 5 = 180 \\ 649^2 - 13 \times 180^2 = 421201 - 421200 = 1 \end{array}\)

循環節の長さが偶数の場合、\(x^2 - m y^2\) は \(+1\) になり、奇数の場合は \(-1\) になります。

これを Common Lisp でプログラムすると、次のようになります。

リスト : ペル方程式

(defun pells-equation (m)
  (let ((cf (butlast (cont-frac-root m))))
    (unless cf
      (error "pells-equation: domain error"))
    (let ((p (cont-frac-sub cf 1 (car cf)))
          (q (cont-frac-sub cf 0 1)))
      (if (= (- (* p p) (* m q q)) 1)
          (values p q)
        (values (+ (* p p) (* m q q)) (* 2 p q))))))

最初に cont-frac-root で \(\sqrt m\) の連分数展開を求め、最後の要素を butlast で削除します。もし、cf が空リストであれば、m は平方数なのでエラーを送出します。あとは、cont-frac-sub で近似値の分子と分母を求め、変数 p, q にセットします。式の値が 1 ならば、values で p と q を返します。そうでなければ、\(p + q\sqrt m\) の二乗を計算して values で返します。

それでは実行してみましょう。

* (pells-equation 2)
3
2
* (pells-equation 3)
2
1
* (pells-equation 7)
8
3
* (pells-equation 61)
1766319049
226153980
* (pells-equation 109)
158070671986249
15140424455100

正常に動作しているようです。興味のある方はいろいろ試してみてください、

●参考文献・URL

奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
平方根の連分数とペル方程式 (PDF), (有澤健治さん)
連分数展開とその計算方法, (マスオさん)
ペル方程式に関する基本的な性質まとめ, (マスオさん)
連分数 - Wikipedia
ペル方程式 - Wikipedia

●プログラムリスト

;;;
;;; frac.py : 分数 (連分数と循環小数)
;;;
;;;           Copyright (C) 2024 Makoto Hiroi
;;;

;;;
;;; 循環小数
;;;

;;; 分数を循環小数に変換する
(defun repeat-decimal (x)
  (unless (rationalp x)
    (error "must be rational number"))
  (do ((m (numerator x))
       (n (denominator x))
       (xs nil) (ys nil))
      (nil)
      (multiple-value-bind
       (p q)
       (floor m n)
       (cond
        ((zerop q)
         (push p ys)
         (return (list (nreverse ys) '(0))))
        (t
         (let ((i (position q xs)))
           (push p ys)
           (when (and i (>= i 0))
             (return (list (nreverse (subseq ys (1+ i)))
                           (nreverse (subseq ys 0 (1+ i))))))
           (push q xs)
           (setq m (* q 10))))))))

;;; 循環小数を分数に直す
(defun from-repeat-decimal (xs ys)
  (let ((p0 (reduce (lambda (a x) (+ (* a 10) x)) xs :initial-value 0))
        (q0 (expt 10 (1- (length xs))))
        (p1 (reduce (lambda (a x) (+ (* a 10) x)) ys :initial-value 0))
        (q1 (1- (expt 10 (length ys)))))
    (/ (+ (* q1 p0) p1) (* q0 q1))))

;;; 単位分数の循環節の長さを求める
(defun solver (d)
  (do ((k 0)
       (m 0)
       (d d (1- d)))
      ((>= k d) (list m k))
      (let* ((xs (repeat-decimal (/ d)))
             (l (length (second xs))))
        (when (> l k)
          (setf k l m d)))))

;;; 循環節の長さを求める
(defun repeat-length (d)
  (do ((n 1 (1+ n))
       (m (mod 10 d) (mod (* m 10) d)))
      ((= m 1) n)))

(defun solver1 (d)
  (when (evenp d) (decf d))
  (do ((k 0)
       (m 0)
       (d d (- d 2)))
      ((>= k d) (list m k))
      (unless (zerop (mod d 5))
        (let ((n (repeat-length d)))
          (when (> n k)
            (setf k n m d))))))

;;;
;;; 連分数
;;;

;;; 分数を正則連分数に変換
(defun cont-frac (x)
  (unless (rationalp x)
    (error "must be rational number"))
  (do ((a (numerator x))
       (b (denominator x))
       (xs nil))
      ((<= b 0) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor a b)
       (push p xs)
       (setf a b b q))))

;;; 正則連分数を有理数に戻す
(defun cont-frac-reduce (xs)
  (if (null (cdr xs))
      (car xs)
    (+ (car xs)
       (/ (cont-frac-reduce (cdr xs))))))

(defun cont-frac-reduce1 (xs)
  (reduce (lambda (x a) (+ x (/ a))) xs :from-end t))

(defun cont-frac-sub (xs a0 a1)
  (cdr (reduce (lambda (a x) (cons (cdr a) (+ (* x (cdr a)) (car a))))
               (cdr xs)
               :initial-value (cons a0 a1))))

(defun cont-frac-reduce2 (xs)
  (/ (cont-frac-sub xs 1 (car xs)) (cont-frac-sub xs 0 1)))

;;; 実数を正則連分数に展開する
(defun cont-frac-flt (x &key (k 8) (e 1e-16))
  (do ((k k (1- k))
       (xs nil))
      ((zerop k) (nreverse xs))
      (multiple-value-bind
       (p q)
       (floor x)
       (push p xs)
       (when (< (abs q) e)
         (return (nreverse xs)))
       (setf x (/ q)))))

;;; √ m の連分数展開
(defun cont-frac-root (m)
  (let* ((c0 (isqrt m))
         (cf (list c0)))
    (if (= (* c0 c0) m)
        cf
      (let ((a 1) (b 0) (c c0))
        (loop
         (setf b (- (* c a) b)
               a (floor (- m (* b b)) a)
               c (floor (+ c0 b) a))
         (push c cf)
         (when (= a 1)
           (return (nreverse cf))))))))

;;; ペル方程式 p^2 - m*q^2 = 1
(defun pells-equation (m)
  (let ((cf (butlast (cont-frac-root m))))
    (unless cf
      (error "pells-equation: domain error"))
    (let ((p (cont-frac-sub cf 1 (car cf)))
          (q (cont-frac-sub cf 0 1)))
      (if (= (- (* p p) (* m q q)) 1)
          (values p q)
        (values (+ (* p p) (* m q q)) (* 2 p q))))))

Common Lisp Programming