rationalizeのアルゴリズムが全く分からなかった

Common Lispの関数rationalizeの話題です。

関数rationalizeとは、浮動小数を分数に変換する関数です。
例えば浮動小数の0.1を分数に直すと次のようになります。

* (rationalize 0.1)
1/10

実は浮動小数の0.1を1/10に変換するって難しい事なのです。
0.25を1/4に変換するなら簡単です。
0.25は基数2の浮動小数では自然に表せる数ですから。

0.1は違います。
浮動小数の0.1というのは、基数が2の浮動小数では 正確に表すことができずに仮数が循環してしまうものとして有名です。 例えば十進数で10を3で割ると0.3333...になるのと同じです。

0.1を単純に分数へ変換するにはrational関数を使います。

* (rational 0.1)
13421773/134217728

すごい値ですね。
でもちゃんと0.1なんです。

* (coerce 13421773/134217728 'float)
0.1

関数rationalの実装は簡単です。 仮数×2の指数乗をそのまま分数形式にして通分するだけですから。

しかし関数rationalizeは違います。
どうも浮動小数0.1と誤差の範囲から適切なものを探し出しているようなのです。
ではどうやって実装されているのでしょうか。

考えてもわかりませんでした。
仕方がないのでcmuclとsbclのソースを見てみました。

• cmucl
  • https://gitlab.common-lisp.net/cmucl/cmucl/blob/21d/src/code/float.lisp#L1618
• sbcl
  • https://sourceforge.net/p/sbcl/sbcl/ci/sbcl-1.5.0/tree/src/code/float.lisp#l863

cmuclのソースを見ると、1618行目から下記のように始まっています。

;;; RATIONALIZE  --  Public
;;;
;;; The algorithm here is the method described in CLISP.  Bruno Haible has
;;; graciously given permission to use this algorithm.  He says, "You can use
;;; it, if you present the following explanation of the algorithm."

ヒュー！！　カッコいい！！

たぶんcmuclさんがrationalizeのアルゴリズムを使用したいってお願いして、 clispさんがいいよって言ったのかなって思います。

Bruno Haibleという名前には見覚えがあります。
なにせclispを起動するときに毎回出てくるのですから。
いつもお世話になっております。

$ clisp
  i i i i i i i       ooooo    o        ooooooo   ooooo   ooooo
  I I I I I I I      8     8   8           8     8     o  8    8
  I  \ `+' /  I      8         8           8     8        8    8
   \  `-+-'  /       8         8           8      ooooo   8oooo
    `-__|__-'        8         8           8           8  8
        |            8     o   8           8     o     8  8
  ------+------       ooooo    8oooooo  ooo8ooo   ooooo   8

Welcome to GNU CLISP 2.49 (2010-07-07) <http://clisp.cons.org/>

Copyright (c) Bruno Haible, Michael Stoll 1992, 1993
Copyright (c) Bruno Haible, Marcus Daniels 1994-1997
Copyright (c) Bruno Haible, Pierpaolo Bernardi, Sam Steingold 1998
Copyright (c) Bruno Haible, Sam Steingold 1999-2000
Copyright (c) Sam Steingold, Bruno Haible 2001-2010

Type :h and hit Enter for context help.

[1]>

こういうところで、Lisp仲間がつながっているのですね。

clispのソースはここで確認できます。

• clisp
  • https://sourceforge.net/p/clisp/clisp/ci/clisp_2_49_92-2018-02-18/tree/src/realelem.d#l1544

cmucl, sbcl, clispのソースには、コメントにアルゴリズムが載っています。

そのアルゴリズムを見ていたのですが、全く分からない。
sbclのソースには「An introduction to number theory」を見ろと書かれており、 さらにharvard.eduのリンク切れのサイトが紹介されている。 何となく思うんですが、 分数の割り算でつまずいている自分が 理解するのは非常に困難なんじゃないか。

理解できないなら既存のソースをたどって作るしかありません。 cmuclのソースはPublic Domainとの事なので、ありがたく参考にさせていただきます。

nptではcmuclのlispコードをもとにC言語で実装しました。 - https://github.com/nptcl/npt/blob/v0.1.6/src/real_decode.c#L746

他の関数はC言語に極力逆らわないように実装しているのですが、 関数rationalizeに関しては、C言語の皮をかぶったlispになっております。 しかしちゃんと動作しているようなので安心しました。