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: 数値

Last modified : 2012/02/01 20:08:31 UTC