Gauche:数値の入出力
浮動小数点数の読み込み
参考:Will Clinger, "How to Read Floating Point Numbers Accurately", ACM SIGPLAN '90, pp.92--101, 1990.
(2002/04/03 03:36:25 PST) Gauche-0.5.2までは、浮動小数点数は文字列として 読み込んだあとsscanfしていた。 とりあえずmantissaとexponentをintegerで読んで自分で計算するように してみる (number.c, v1.66)。 意外や意外、ベンチマークでは34%も速くなった (on Linux/i686)。 以前のルーチンではDStringを使っていたせいか?
(2002/04/09 13:28:20 PDT) double precisionの有効桁数を越えた場合にAlgorithmRを使うようにしてみる。 正確に読むようにはなったが、多くの場合でbignum演算をするようになったので 性能が心配。
(2002/04/16 14:30:28 PDT) 0.5.3リリース時で、ランダムな浮動小数点数の読み込み速度は0.5.2のだいたい半分。 もっともこれはほとんど全ての数値が最大桁数まである場合で、 入力桁数が少なければalgorithmRは通らないからむしろ速くなる。
Bignumの読み込み
(2002/04/03 03:36:25 PST) 改善すべき点
- intの範囲を越えたらbignumにスイッチするが、その後はひとけたづつ
bignumの乗算と加算が入るのでやな感じ。intで表せるchunk毎に
bignum演算を行うようにできないか。
- 2002/04/04 00:18:51 PST 実装してみた (number.c, v1.69)。bignumの読み込み速度2.6倍に。 (max96bit長のランダムな10進表記のbignumをstring->numberで読ませて計測 on Linux/i686)。
浮動小数点数の書き出し
参考: Guy L. Steele Jr. and Jon L. White, "How to Print Floating-Point Numbers Accurately", in the Proceedings of ACM SIGPLAN '90, pp.112-126, 1990. Robert G. Burger and R. Kent Dybvig, "Printing Floating-Point Numbers Quickly and Accurately", PLDI '96, pp.108--116, 1996.
(2002/04/06 02:29:05 PST) 素直にBurger&Dybvigを実装してみた(number.c, v1.75)が…
gosh> 1e23 0.9999999999999999e23 gosh> (* 5960464477539062 (expt 2.0 24)) 0.9999999999999999e23 gosh> (* 5960464477539063 (expt 2.0 24)) 1.0000000000000001e23
うーむ…
(2002/04/06 02:40:26 PST) あ、嘘。rounding modeをちゃんと考慮してなかった。 直した(number.c, v1.76)らちゃんと動いた。
gosh> 1e23 1.0e23 gosh> (* 5960464477539062 (expt 2.0 24)) 1.0e23 gosh> (* 5960464477539063 (expt 2.0 24)) 1.0000000000000001e23
次はオプティマイズかな。
(2002/04/06 03:19:02 PST) あら。SourceForgeが落ちている。v1.76がチェックインできない。 書き出しのベンチマークを取ってみたら、Burger&Dybvigはprintfの1/24しか 速度が出てない。全然最適化してないから無理もないが。
(2002/04/16 14:30:28 PDT) 一度はSCMにならってBurger&Dybvigをdoubleで実装してみたが、 プラットフォームによって動作が変わったりするのでやめた。 printfより速かったのだが。
bignum演算をかりかりチューニングして、0.5.3の時点でprintfの1/10くらい まできた。多分もっといける。
要注意な数値
1.0e23
1.0e23は、倍精度浮動小数点数(mantissa 53bit)で表せる二つの数値のちょうど 真中にある。
99999999999999991611392 = 5960464477539062 * 2^24 (a) 100000000000000008388608 = 5960464477539063 * 2^24 (b)
even方向へ丸めるならば、1.0e23は(a)のように読まれるべきで、 それを書き出すと1.0e23になるべき。
3.141592653589765e-20
Burger&Dybvigをfloating point numberで近似した時に、このへんの値が 正しくprintされない。
Burger&Dybvig 近似 5219866133129672 * 2^-117 : 3.141592653589764e-20 3.141592653589764e-20 OK 5219866133129673 * 2^-117 : 3.1415926535897647e-20 3.141592653589765e-20 NG 5219866133129674 * 2^-117 : 3.141592653589765e-20 3.1415926535897654e-20 NG 5219866133129675 * 2^-117 : 3.141592653589766e-20 3.141592653589766e-20 OK
理由: scalingの際のroundingによる誤差。
v = 3.141592653589765e-20 -> #x141b2f769cf0ae * 2^-117 scale = 10^19 = #x8ac7230489e80000 r = v * scale = #xa0d97bb4e7856ee35a1d0f9100000 * 2^-117 (round to 53bit mantissa) -> #xa0d97bb4e78570000000000000000 * 2^-117 = 0.31415926535897653604...
しかし、bignumを使わずfloating point numberのみでprintルーチンを構成 しているSCMではちゃんとfloating point numberのwrite-read invarianceを 実現している。その謎を解析してみた→ Gauche:拡張浮動小数点演算の謎
Tag: 数値