浮動小数点について復習した

SICPの練習問題1.7を解くに辺り、浮動小数点の復習を行った。
というか、割と適当に流してたけど、この問題は奥が深いんだな。

SICPの練習問題1.7自体はこちらを見てもらうとして、まずSICPを解くにあたり自分が利用しているGaucheでは、IEEEの64bit浮動小数点数を用いているらしい。

浮動小数点の不正確な実数
実装に使われるC言語のdouble型で表現されます。通常IEEE 64bit浮動小数点数です。
Gauche ユーザリファレンス: 6.2 等価

で、浮動小数点の復習。
日経PC21 / 演算誤差の正体 - IEEE 754 浮動小数点数の仕組みが、分かりやすかった。
今回の問題を解くためのポイントは符号部1bit、指数部11bit、仮数部52bitで構成されているという点。
数値の絶対値部分を表現するのは仮数部の52bitを用いる。

練習問題1.7では、更新の打ち切りのための閾値に0.001が用いられていて、これを2進数で表すためには、2^{-9} = 0.0019531252^{-10} = 0.0009765625から仮数部に10ビットが必要となる。
すなわち残りの42ビットのみを利用して表現できる数が、0.001と同時に表現できる数値となる。すなわち2^{41} = 2199023255552より、10^{12}桁の数字までなら更新式が成り立つ。
それ以上の値では、0.001が表現できなくなり、更新式が終了しなくなる。


ここまで考えてきて、色々とすっきりしたけど、唯一分からない現象が残っている。

  • 偶数倍の値の時(sqrt 1e14), (sqrt 1e16)...なら、(sqrt 1e48)まで処理が行える(奇数の場合はダメ)

色々考えたけどよくわからない。
この辺りの話って皆つまずかないのかな?すごく気持悪い。

追記

@haruyamaさんに教えてもらったリンクより、

x86系の浮動小数点演算ユニットのレジスタは80bitの拡張浮動小数点数 (仮数部64bit)で、floatでもdoubleでも演算は一度80bitに直してから行われ、 結果が単精度または倍精度に丸められる。 ところが、ある種の最適化が行われた場合、 途中結果が浮動小数レジスタに置いたままにされるため、80bitのまま 演算が進行する。

その値をスタック経由で関数に渡す時にはdoubleに丸められるので、 printfに渡して印字させるだけでは違いがわからなかったのだ。
http://practical-scheme.net/wiliki/wiliki.cgi?Gauche:%E6%8B%A1%E5%BC%B5%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E6%BC%94%E7%AE%97%E3%81%AE%E8%AC%8E

これで、想定していた数値よりも大きなものが処理できている事を部分的には理解できた。
また続きを書くかもしれない。