M.Hiroi's Home Page

Common Lisp Programming

お気楽 ISLisp プログラミング超入門

Copyright (C) 2024 Makoto Hiroi
All rights reserved.

簡単なプログラム

●有理数 (分数)

●プログラムリスト

;;;
;;; rational.lsp : 有理数 (ISLisp 用)
;;;
;;; Copyright (c) 2024 Makoto Hiroi
;;;
;;; Released under the MIT license
;;; https://opensource.org/license/mit/
;;;

;;; メソッドの宣言
(defgeneric +/ (x y))
(defgeneric -/ (x y))
(defgeneric */ (x y))
(defgeneric // (x y))
(defgeneric inv (x))
(defgeneric =/ (x y))
(defgeneric /=/ (x y))
(defgeneric </ (x y))
(defgeneric <=/ (x y))
(defgeneric >/ (x y))
(defgeneric >=/ (x y))
(defgeneric is-integer (x))
(defgeneric to-integer (x))
(defgeneric to-float (x))
(defgeneric to-string (x))
(defgeneric display/ (x))
(defgeneric repeat-decimal (x))
(defgeneric from-repeat-decimal (x y))

;;; 有理数
(defclass <rational> ()
  ((num :accessor numerator   :initform 0 :initarg num)
   (den :accessor denominator :initform 1 :initarg den)))

;;; コンストラクタ
(defun make-ratio (a b)
  (if (= b 0)
      (error "divizion by zero"))
  (let ((f (if (< b 0) -1 1))
        (z (gcd (abs a) (abs b))))
    (create (class <rational>) 'num (div (* f a) z) 'den (div (* f b) z))))

;;; 型述語
(defun rationalp (x)
  (eq (class-of x) (class <rational>)))

;;; 算術演算子
(defmethod +/ ((x <rational>) (y <rational>))
  (let ((x1 (numerator x)) (x2 (denominator x))
        (y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (+ (* x1 y2) (* x2 y1)) (* x2 y2))))

(defmethod +/ ((x <rational>) (y <integer>))
  (let ((x1 (numerator x)) (x2 (denominator x)))
    (make-ratio (+ x1 (* x2 y)) x2)))

(defmethod +/ ((x <integer>) (y <rational>))
  (let ((y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (+ (* x y2) y1) y2)))

(defmethod +/ ((x <integer>) (y <integer>))
  (make-ratio (+ x y) 1))

(defmethod -/ ((x <rational>) (y <rational>))
  (let ((x1 (numerator x)) (x2 (denominator x))
        (y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (- (* x1 y2) (* x2 y1)) (* x2 y2))))

(defmethod -/ ((x <rational>) (y <integer>))
  (let ((x1 (numerator x)) (x2 (denominator x)))
    (make-ratio (- x1 (* x2 y)) x2)))

(defmethod -/ ((x <integer>) (y <rational>))
  (let ((y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (- (* x y2) y1) y2)))

(defmethod -/ ((x <integer>) (y <integer>))
 (make-ratio (- x y) 1))

(defmethod */ ((x <rational>) (y <rational>))
  (let ((x1 (numerator x)) (x2 (denominator x))
        (y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (* x1 y1) (* x2 y2))))

(defmethod */ ((x <rational>) (y <integer>))
  (let ((x1 (numerator x)) (x2 (denominator x)))
    (make-ratio (* x1 y) x2)))

(defmethod */ ((x <integer>) (y <rational>))
  (let ((y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (* x y1) y2)))

(defmethod */ ((x <integer>) (y <integer>))
  (make-ratio (* x y) 1))

(defmethod // ((x <rational>) (y <rational>))
  (let ((x1 (numerator x)) (x2 (denominator x))
        (y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (* x1 y2) (* x2 y1))))

(defmethod // ((x <rational>) (y <integer>))
  (let ((x1 (numerator x)) (x2 (denominator x)))
    (make-ratio x1 (* x2 y))))

(defmethod // ((x <integer>) (y <rational>))
  (let ((y1 (numerator y)) (y2 (denominator y)))
    (make-ratio (* x y2) y1)))

(defmethod // ((x <integer>) (y <integer>))
  (make-ratio x y))

;;; 逆数
(defmethod inv ((x <integer>)) (make-ratio 1 x))

(defmethod inv ((x <rational>))
  (make-ratio (denominator x) (numerator x)))

;;; 有理数の比較
(defun compare-ratio (x y)
  (let ((x1 (numerator x)) (x2 (denominator x))
        (y1 (numerator y)) (y2 (denominator y)))
    (- (* x1 y2) (* x2 y1))))

;;; 比較演算子の定義
(defmethod =/ ((x <rational>) (y <rational>))
  (= (compare-ratio x y) 0))

(defmethod /=/ ((x <rational>) (y <rational>))
  (not (= (compare-ratio x y) 0)))

(defmethod </ ((x <rational>) (y <rational>))
  (< (compare-ratio x y) 0))

(defmethod >/ ((x <rational>) (y <rational>))
  (> (compare-ratio x y) 0))

(defmethod <=/ ((x <rational>) (y <rational>))
  (<= (compare-ratio x y) 0))

(defmethod >=/ ((x <rational>) (y <rational>))
  (>= (compare-ratio x y) 0))

;;; 整数の判定
(defmethod is-integer ((x <rational>))
  (= (denominator x) 1))

;;; 整数に変換
(defmethod to-integer ((x <rational>))
  (div (numerator x) (denominator x)))

;;; 浮動小数点数に変換
(defmethod to-float ((x <rational>))
  (quotient (numerator x) (denominator x)))

;;; 文字列に変換
(defmethod to-string ((x <rational>))
  (string-append (convert (numerator x) <string>) "/" (convert (denominator x) <string>)))

;;; 表示
(defmethod display/ ((x <rational>))
  (format (standard-output) "~D/~D" (numerator x) (denominator x)))

;;; 検索
(defun find-index (x xs i)
  (cond
   ((null xs) -1)
   ((eql x (car xs)) i)
   (t
    (find-index x (cdr xs) (+ i 1)))))

;;; 循環小数 m/n = ((...) (...))
(defmethod repeat-decimal ((x <rational>))
  (let ((m (numerator x))
        (n (denominator x))
        (xs nil) (ys nil) (zs nil))
    (while (not zs)
      (let ((p (div m n))
            (q (mod m n)))
        (setq ys (cons p ys))
        (cond
         ((= q 0)
          (setq zs (list (reverse ys) '(0))))
         (t
          (let ((i (find-index q xs 0)))
            (cond
             ((>= i 0)
              (setq zs (list (reverse (subseq ys (+ i 1) (length ys)))
                             (reverse (subseq ys 0 (+ i 1))))))
             (t
              (setq xs (cons q xs))
              (setq m (* q 10)))))))))
    zs))

;;; 循環小数を分数に直す
(defun digit-to-number (xs a)
  (if (null xs)
      a
    (digit-to-number (cdr xs) (+ (* a 10) (car xs)))))

(defmethod from-repeat-decimal ((xs <list>) (ys <list>))
  (let ((p0 (digit-to-number xs 0))
        (q0 (expt 10 (- (length xs) 1)))
        (p1 (digit-to-number ys 0))
        (q1 (- (expt 10 (length ys)) 1)))
    (// (+ (* q1 p0) p1) (* q0 q1))))

●使用例

ISLisp>(load "rational.lsp")
T
ISLisp>(defglobal a (// 1 2))
A
ISLisp>(display/ a)
1/2NIL
ISLisp>(rationalp a)
T
ISLisp>(to-string a)
"1/2"
ISLisp>(defglobal b (// 1 3))
B
ISLisp>(display/ b)
1/3NIL
ISLisp>(display/ (+/ a b))
5/6NIL
ISLisp>(display/ (-/ a b))
1/6NIL
ISLisp>(display/ (-/ b a))
-1/6NIL
ISLisp>(display/ (*/ a b))
1/6NIL
ISLisp>(display/ (// a b))
3/2NIL
ISLisp>(display/ (// b a))
2/3NIL

ISLisp>(=/ a b)
NIL
ISLisp>(=/ a a)
T
ISLisp>(/=/ a b)
T
ISLisp>(/=/ a a)
NIL
ISLisp>(</ a b)
NIL
ISLisp>(</ b a)
T
ISLisp>(>/ a b)
T
ISLisp>(>/ b a)
NIL

ISLisp>(is-integer a)
NIL
ISLisp>(is-integer (+/ a a))
T
ISLisp>(to-integer (+/ a a))
1
ISLisp>(to-integer (+/ a b))
0
ISLisp>(to-integer (// a b))
1
ISLisp>(to-float (+/ a b))
0.833333
ISLisp>(to-float (// a b))
1.5

ISLisp>(repeat-decimal (// 1 2))
((0 5) (0))
ISLisp>(repeat-decimal (// 1 3))
((0) (3))
ISLisp>(repeat-decimal (// 1 7))
((0) (1 4 2 8 5 7))
ISLisp>(repeat-decimal (// 355 113))
((3) (1 4 1 5 9 2 9 2 0 3 5 3 9 8 2 3 0 0 8 8 4 9 5 5 7 5 2 2 1 2 3 8 9 3 8 0 5 
3 0 9 7 3 4 5 1 3 2 7 4 3 3 6 2 8 3 1 8 5 8 4 0 7 0 7 9 6 4 6 0 1 7 6 9 9 1 1 5 
0 4 4 2 4 7 7 8 7 6 1 0 6 1 9 4 6 9 0 2 6 5 4 8 6 7 2 5 6 6 3 7 1 6 8))
ISLisp>(display/ (apply #'from-repeat-decimal (repeat-decimal (// 1 2))))
1/2NIL
ISLisp>(display/ (apply #'from-repeat-decimal (repeat-decimal (// 1 3))))
1/3NIL
ISLisp>(display/ (apply #'from-repeat-decimal (repeat-decimal (// 1 7))))
1/7NIL
ISLisp>(display/ (apply #'from-repeat-decimal (repeat-decimal (// 355 113))))
355/113NIL
-- [注意] --------
Easy-ISLisp の場合、ver5.37 では正常に動作しないことがあります。rational.lsp を使用する場合は ver5.39 にバージョンアップしてください。迅速に対応していただいた笹川さんに感謝いたします。

●分数を使ったパズル

それでは rational.lsp の簡単な例題として、分数を使ったパズルを紹介しましょう。rational.lsp を使うと、簡単に解法プログラムを作ることができます。

●パズル「小町分数」

それでは問題です。

下図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。3 つの分数を足すと 1 / N になる配置を求めてください。

\( \dfrac{A}{BC} + \dfrac{D}{EF} + \dfrac{G}{HI} = \dfrac{1}{N} \)

ex)
3/27 + 6/54 + 9/81 = 1/3
3/54 + 6/72 + 9/81 = 1/4

図 1 : 小町分数 (1)

このパズルの元ネタは N = 1 の場合で、参考文献『超々難問数理パズル 解けるものなら解いてごらん』に掲載されています。ちなみに、3 つの分数の和が整数になる場合、その値は 1 しかありません。また、値が 1 / N (N は整数) になる場合は 2, 3, 4, 6, 10 の 5 通りです。

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

リスト : パズル「小町分数」の解法

(load "combination.lsp")  ; Easy-ISLisp は (import "combination") とする
(load "rational.lsp")

(defun check (a b c d e f g h i)
  (let ((x (// a (+ (* b 10) c)))
        (y (// d (+ (* e 10) f)))
        (z (// g (+ (* h 10) i))))
    (if (and (</ x y) (</ y z))
        (let ((s (+/ (+/ x y) z)))
          (if (= (numerator s) 1)
              (format (standard-output) "~D/~D~D + ~D/~D~D + ~D/~D~D = ~A~%"
                      a b c d e f g h i (to-string s)))))))

(defun solver ()
  (permutation (lambda (xs) (apply #'check xs)) 9 '(1 2 3 4 5 6 7 8 9)))

拙作のライブラリ combination.lsp の関数 permutation を使って順列を生成し、分数の和が 1 / N になる数字を探します。単純な生成検定法ですが、分数 x, y, z を x < y < z に限定することで、重複解を生成しないように工夫しています。それから、このプログラムでは 3 つの分数の和が 1 になる場合も解を出力します。プログラムを実行するときはご注意くださいませ

●実行結果

ISLisp>(load "komachi_ratio.lsp")
T
ISLisp>(solver)
1/24 + 3/56 + 7/98 = 1/6
1/26 + 7/84 + 5/39 = 1/4
1/32 + 5/96 + 7/84 = 1/6
1/48 + 7/96 + 5/32 = 1/4
1/56 + 3/72 + 9/84 = 1/6
1/96 + 5/48 + 7/32 = 1/3
1/96 + 7/84 + 5/32 = 1/4
2/95 + 1/38 + 4/76 = 1/10
3/48 + 9/72 + 5/16 = 1/2
3/54 + 6/72 + 9/81 = 1/4
4/57 + 2/19 + 6/38 = 1/3
5/63 + 2/18 + 7/49 = 1/3
7/68 + 5/34 + 9/12 = 1/1
NIL

●パズル「小町分数 (2)」

もうひとつ「小町分数」を出題しましょう。

下図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。3 つの分数を足すと 1 / 2 になる配置を求めてください。

\( \dfrac{A}{B \times C} + \dfrac{D}{E \times F} + \dfrac{G}{H \times I} = \dfrac{1}{2} \)

図 2 : 小町分数 (2)

このパズルの元ネタも値が 1 になる場合で、参考文献『超々難問数理パズル 解けるものなら解いてごらん』に掲載されています。この問題で値が 1 / N (N は整数) になる場合は 1 と 2 の 2 通りしかないようです。

●プログラムと実行結果

リスト : 小町分数 (2)

(defun check2 (a b c d e f g h i)
  (let ((x (// a (+ (* b c))))
	(y (// d (+ (* e f))))
	(z (// g (+ (* h i)))))
    (if (and (< a d) (< d g) (< b c) (< e f) (< h i))
        (let ((s (+/ (+/ x y) z)))
          (if (= (numerator s) 1)
	      (format (standard-output) "~D/(~D*~D) + ~D/(~D*~D) + ~D/(~D*~D) = ~A~%"
                      a b c d e f g h i (to-string s)))))))

(defun solver2 ()
  (permutation (lambda (xs) (apply #'check2 xs)) 9 '(1 2 3 4 5 6 7 8 9)))

ISLisp>(solver2)
1/(2*4) + 5/(3*6) + 7/(8*9) = 1/2
1/(3*6) + 5/(8*9) + 7/(2*4) = 1/1
NIL

●連分数

連分数 (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}\)

OK!-ISLisp で検算してみましょう。 

ISLisp>(load "rational.lsp")
T
ISLisp>(display/ (+/ 1 (inv (+/ 1 (inv (+/ 1 (inv 5)))))))
17/11NIL
ISLisp>(display/ (+/ 3 (inv (+/ 7 (inv 16)))))
355/113NIL

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

●プログラムの作成

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

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

(defun cont-frac (x)
  (if (not (rationalp x))
      (error "must be rational")
    (for ((a (numerator x))
          (b (denominator x))
          (xs nil))
       ((<= b 0) (reverse xs))
       (let ((p (div a b))
             (q (mod a b)))
         (setq xs (cons p xs))
         (setq a b)
         (setq b q)))))

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

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

(defun cont-frac-reduce (xs)
  (if (null (cdr xs))
      (car xs)
    (+/ (car xs)
       (inv (cont-frac-reduce (cdr xs))))))

;;; 畳み込み
(defun foldr (f xs)
  (if (null (cdr xs))
      (car xs)
    (funcall f (car xs) (foldr f (cdr xs)))))

;;; 畳み込みを使う
(defun cont-frac-reduce1 (xs)
  (foldr (lambda (x a) (+/ x (inv a))) xs))

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

なお、参考 URL 『ペル方程式に関する基本的な性質まとめ』によると、連分数 \([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}} \)

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

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

;;; 畳み込み
(defun foldl (f a xs)
  (if (null xs)
      a
    (foldl f (funcall f a (car xs)) (cdr xs))))

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

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

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

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

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

ISLisp>(display/ (cont-frac-reduce '(1 1 1 5)))
17/11NIL
ISLisp>(display/ (cont-frac-reduce '(-2 2 5)))
-17/11NIL
ISLisp>(display/ (cont-frac-reduce '(3 7 16)))
355/113NIL

ISLisp>(display/ (cont-frac-reduce1 '(1 1 1 5)))
17/11NIL
ISLisp>(display/ (cont-frac-reduce1 '(-2 2 5)))
-17/11NIL
ISLisp>(display/ (cont-frac-reduce1 '(3 7 16)))
355/113NIL

ISLisp>(display/ (cont-frac-reduce2 '(1 1 1 5)))
17/11NIL
ISLisp>(display/ (cont-frac-reduce2 '(-2 2 5)))
-17/11NIL
ISLisp>(display/ (cont-frac-reduce2 '(3 7 16)))
355/113NIL

ISLisp>(display/ (cont-frac-reduce '(1 2 3 4 5 6 7 8)))
81201/56660NIL
ISLisp>(cont-frac (// 81201 56660))
(1 2 3 4 5 6 7 8)

ISLisp>(display/ (cont-frac-reduce1 '(1 2 3 4 5 6 7 8)))
81201/56660NIL
ISLisp>(display/ (cont-frac-reduce2 '(1 2 3 4 5 6 7 8)))
81201/56660NIL

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

●無理数の連分数展開

有理数の連分数展開は有限回で終了しますが、無理数の連分数展開は無限に続きます。たとえば、\(\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}\)

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

ISLisp>(abs (- (to-float (cont-frac-reduce '(1 2 2 2 2))) (sqrt 2)))
0.000420459
ISLisp>(abs (- (to-float (cont-frac-reduce '(1 2 2 2 2 2 2))) (sqrt 2)))
1.23789e-005
ISLisp>(abs (- (to-float (cont-frac-reduce '(1 2 2 2 2 2 2 2 2))) (sqrt 2)))
3.64404e-007

ISLisp>(abs (- (to-float (cont-frac-reduce '(1 1 2 1 2))) (sqrt 3)))
0.00477808
ISLisp>(abs (- (to-float (cont-frac-reduce '(1 1 2 1 2 1 2))) (sqrt 3)))
0.00034349
ISLisp>(abs (- (to-float (cont-frac-reduce '(1 1 2 1 2 1 2 1 2))) (sqrt 3)))
2.46638e-005

ISLisp>(abs (- (to-float (cont-frac-reduce '(2 4 4 4 4))) (sqrt 5)))
2.40373e-006
ISLisp>(abs (- (to-float (cont-frac-reduce '(2 4 4 4 4 4 4))) (sqrt 5)))
7.46507e-009
ISLisp>(abs (- (to-float (cont-frac-reduce '(2 4 4 4 4 4 4 4 4))) (sqrt 5)))
2.31837e-011

ISLisp>(abs (- (to-float (cont-frac-reduce '(2 1 1 1 4))) (sqrt 7)))
0.00289417
ISLisp>(abs (- (to-float (cont-frac-reduce '(2 1 1 1 4 1 1 1 4))) (sqrt 7)))
1.14008e-005
ISLisp>(abs (- (to-float (cont-frac-reduce '(2 1 1 1 4 1 1 1 4 1 1 1 4))) (sqrt 7)))
4.48856e-008

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

浮動小数点数 \(\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 k e)
  (block exit
    (for ((k k (- k 1))
          (xs nil))
         ((= k 0) (reverse xs))
      (let* ((p (floor x))
             (q (- x p)))
        (setq xs (cons p xs))
        (cond
         ((< (abs q) e)
	  (return-from exit (reverse xs)))
         (t
	  (setq x (quotient 1 q))))))))

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

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

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

ISLisp>(cont-frac-flt (sqrt 2) 8 1e-16)
(1 2 2 2 2 2 2 2)
ISLisp>(cont-frac-flt (sqrt 3) 8 1e-16)
(1 1 2 1 2 1 2 1)
ISLisp>(cont-frac-flt (sqrt 5) 8 1e-16)
(2 4 4 4 4 4 4 4)
ISLisp>(cont-frac-flt (sqrt 7) 8 1e-16)
(2 1 1 1 4 1 1 1)
ISLisp>(cont-frac-flt (sqrt 7) 16 1e-16)
(2 1 1 1 4 1 1 1 4 1 1 1 4 1 1 1)

ISLisp>(cont-frac-flt 0.5 8 1e-16)
(0 2)
ISLisp>(cont-frac-flt 0.25 8 1e-16)
(0 4)
ISLisp>(cont-frac-flt 0.125 8 1e-16)
(0 8)

ISLisp>(cont-frac-flt 1.1 8 1e-16)
(1 9 1 112589990684262 2 2 6 46912496118442)
ISLisp>(cont-frac-flt 1.1 8 1e-14)
(1 9 1)
ISLisp>(display/ (cont-frac-reduce '(1 9 1)))
11/10NIL

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

●プログラムリスト

;;;
;;; contfrac.lsp : 連分数 (ISLisp 用)
;;;
;;; Copyright (c) 2024 Makoto Hiroi
;;;
;;; Released under the MIT license
;;; https://opensource.org/license/mit/
;;;
(load "rational.lsp")

;;; 分数を連分数に展開する
(defun cont-frac (x)
  (if (not (rationalp x))
      (error "must be rational")
    (for ((a (numerator x))
          (b (denominator x))
          (xs nil))
       ((<= b 0) (reverse xs))
       (let ((p (div a b))
             (q (mod a b)))
         (setq xs (cons p xs))
         (setq a b)
	 (setq b q)))))

;;; 連分数を分数に変換する

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

;;; 畳み込み
(defun foldr (f xs)
  (if (null (cdr xs))
      (car xs)
    (funcall f (car xs) (foldr f (cdr xs)))))

(defun foldl (f a xs)
  (if (null xs)
      a
    (foldl f (funcall f a (car xs)) (cdr xs))))

;;; 畳み込みを使う
(defun cont-frac-reduce1 (xs)
  (foldr (lambda (x a) (+/ x (inv a))) xs))

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

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

;;; 浮動小数点数を正則連分数に展開する
(defun cont-frac-flt (x k e)
  (block exit
    (for ((k k (- k 1))
          (xs nil))
         ((= k 0) (reverse xs))
      (let* ((p (floor x))
             (q (- x p)))
        (setq xs (cons p xs))
        (cond
         ((< (abs q) e)
	  (return-from exit (reverse xs)))
         (t
	  (setq x (quotient 1 q))))))))