Shiro(2012/10/01 06:04:16 UTC): inexact->exactを改良したら、ひとつバグが表面化。事例としてメモ。
1.0の「次」の浮動小数点数 (Gaucheでは 1.0000000000000002 と表示される) を
表現する最も単純な有理数は
3002399751580332/3002399751580331
となる。ところが、これを非正確数へと変換すると、x86_64では 1.0000000000000002
に戻るのに、i686では 1.0000000000000004 (LSBがもう一つ大きい数) になってしまう、
という問題。
gosh> (define x 3002399751580332) x gosh> (define y 3002399751580331) y
として、64bit浮動小数点数(仮数部52+1bit)での割り算をするとどうなるか。
gosh> (quotient&remainder (* x (expt 2 52)) y) 4503599627370497 1501199875790165
なので、結果の仮数部は次のとおり。
10000000000000000000000000000000000000000000000000001
余り 1501199875790165 は除数 3002399751580331 の半分よりわずかに 小さいので、余りの桁は丸めで切り捨てられる。従って、得られる結果は 希望通り、1.0から1LSBだけ大きな数となる。
ところで、これを80bit浮動小数点数(仮数部64+1bit)で割り算するとどうなるか。
gosh> (quotient&remainder (* x (expt 2 64)) y) 18446744073709557759 3002399751578283
結果の仮数部は
10000000000000000000000000000000000000000000000000001011111111111
^
53bit目
で、余り 3002399751578283 は除数 3002399751580331 の半分より大きいため、 上の80bitの結果は切り上げられて、80bit浮動小数点数レジスタに残る値は こうなる。
10000000000000000000000000000000000000000000000000001100000000000
^
53bit目
これをさらに64bit浮動小数点数に変換すると、端数部分がちょうど真ん中なので、 偶数丸め規則により繰り上げが発生して仮数部はこうなる。
10000000000000000000000000000000000000000000000000010
てなわけで、i386では二重丸めにより1LSBだけ結果が大きくなってしまったのだった。
まあもともと浮動小数点数演算には誤差があるので、この計算自体はバグでもなんでもない。 Gaucheの「正確数→非正確数」変換ルーチンがこのコードパスを使っていた、というのがバグ。 正しく非正確数→正確数→非正確数のラウンドトリップを実現するには、 浮動小数点数演算ではない方法で除算して結果を得ないとならない。
Tag: 数値