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のソースを見ると、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のソースはここで確認できます。
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になっております。
しかしちゃんと動作しているようなので安心しました。