2011/01/12

拡張浮動小数点数の扱い

IEEE倍精度浮動小数点数で最大の非正規化数の10進表現である 2.2250738585072011e-308 を パーズしようとすると PHPがハングする問題 についてあれこれ。

直接の原因はx86プロセッサで浮動小数点数演算が80bit拡張精度で行われることで、 似たような話についてはWiLiKi:Gauche:拡張浮動小数点演算の謎でも触れた。

これについて、「gccのバグだろ?」と主張してるblogがあった。

• The PHP strtod() denial of service bug

そこからリンクされているgccのバグレコード323はもう10年の歴史を持つ:

• Bug 323 - optimized code gives strange floating point results

ここから参照されてる論文が、この問題の良いまとめになっていた。

• David Monniaux: The pitfalls of verifying floating-point computations

以下私見。

なんで80bitで計算するとまずいの?

具体的な問題はいくつかの異なる要因から生じている。

1. 二重丸め問題。計算の途中結果を80bitで丸め、それをさらにdouble型の変数にストアする 時に64bitで丸めると、最後に64bitで一括して丸めた場合とは異なる結果になる場合がある。 説明の簡略化のため、本来の仮数部精度を4bit、拡張精度を6bitとして例をあげてみよう。 #bは2進数を表す。
   • x = #b0.111 のとき、1/x = #b1.00010001... となる。有効数字4bitで これに近いのは #b1.000 か #b1.001 でどちらかに丸めることになるが、 1/xは中間 #b1.0001よりも後者に近いので結果は #b1.001になるべき。
   • ところが、コプロセッサは有効数字6bitで計算したので1/xの結果は #b1.00010となった。 これを4bitの変数に代入すると、この値はちょうど #b1.000と #b1.001の中間なので 偶数側に丸める規則を使えば #b1.000となる。
2. 指数部の範囲の問題。doubleの範囲では計算の途中結果がオーバーフローしたり アンダーフローしたりするはずなのに、拡張精度では表現できる場合。 doubleで計算してれば無限大とかNaNとか0になるはずの結果が、普通の数値になる。
3. ソース上で同一に見える数値が、それが置かれるレジスタによって異なる表現を持つ問題。 x/yを2箇所で計算してて、xもyも変えてなければ結果は == になると期待してしまうが、 一方が64bitに丸められてメモリにストアされ、 もう一方が80bitのレジスタに載ったままである場合、 前者を80bitに拡張して浮動小数点数ユニットで比較すると等しくならない。
   • WiLiKi:Gauche:拡張浮動小数点演算の謎の、途中結果をprintfしたら 同じなのにその後の計算結果が違ってくるというのもこのカテゴリで、 printfが受け取る「途中結果」と、実際に計算に使われてる「途中結果」が ソース上では同じなのに実は違うものであることによる。
   • 上のgccバグレポートにはもうちょっと深刻な例が。 この問題のせいで、(x < y && y < x) がtrueになる場合があり、 アルゴリズムが正しく動かない。

どうすれば防げるの?

上記のようにいくつかの異なる問題が混在しているので、一律にこれをやっとけば 解決という対応は難しい。

• 浮動小数点ユニットのモードレジスタを変更して仮数部を doubleの精度に制限することはできるが、指数部の範囲は拡張精度のままなので、 上記2の問題は残る。
• 変数にvolatileつけたりgccの-ffloat-storeオプション つけるのは、変数にアサインされる時点でdoubleになることを保証するけれど、 計算の途中結果が80bitのままになることは防げない。
• 実行時に全部64bitに載ることを保証できたとしても、コンパイラは最適化の 過程でいくつかの値をコンパイル時に計算してしまうかもしれず、その計算は 拡張精度で行われるかもしれない。
• 全部long doubleで計算することにすれば予想外の結果は出ない。 long doubleの精度がアーキテクチャによって異なることと、 アプリのコードが全部自分でコントロールできない場合に難しいことが 問題になるかも。
• C99には中間結果もちゃんと64bitでやるように指示するpragmaが指定されてるけど あんまり実装されていない (少なくともDavid Monniauxの論文の時点では)。

なんでgccは直さないの?

ここはgccチームの見解じゃなくて勝手な推測だけど。

• 「中間結果も全てIEEE doubleで計算した場合に得られるはずの理論値と一致すること」と、 「ハードが持つ最大限の精度をできるだけ使って得られる値」と、どちらが「より正しい」 かはアプリケーションによる。そして、前者を保証するオーバヘッドはかなり大きい
• もっとも、x87由来の浮動小数点数ユニットを使うのをやめて、SSEの方を使えば 64bitのまま計算できる。歴史的にみると、ハードウェア浮動小数点数ユニットが コプロセッサとして実現されだした頃 (8087とか、68881とか) は、精度について 「大は小を兼ねる」と思われていた節がある。 次第に、一貫した精度で計算することが重要ってことがわかってきて、 浮動小数点数ユニットはそういうモードを備えるようになってきた。 それならx87のフェードアウトを待つってことでもいいんじゃないか。

浮動小数点数演算てそもそも不正確なものでしょ。こまけぇことはいいんだよ

浮動小数点数演算でも、特定の条件のもとでは正確な値が計算できるし、 どうなったらどのくらい精度が失われるかも評価できる。なので、 「不正確だから」で考えを止めてしまうのはあんまり良くない。

この80bit問題の根っこは、浮動小数点数が不正確だからということよりも、 浮動小数点数演算に使われる精度をプログラマが完全に制御できない、 あるいはプログラマの直感と違う精度が使われる、というところにあると思う。

良くあるmyth:

• 「浮動小数点数演算は常に誤差を持つ」: 結果が53bitに収まる整数演算(加減乗算、 および余りの出ない除算)とか、指数部が溢れない範囲での2^n でのスケーリングとかは 誤差なくできる。
• 「浮動小数点数演算を10進数で読み書きすると不正確」: 任意の10進小数を有限桁数の2進小数に変換することはできない、という意味では そうなんだけど。2進浮動小数点数を、「実数全体を有限個の区間に分割し、 それぞれを代表値で表現することにしたもの」と考えた場合、 ある区間内の10進小数は、常にその代表値である2進浮動小数点数表現へと 変換でき、また2進浮動小数点数表現は、常にその区間内の10進小数表現へと 変換できる。実数の各区間に対して曖昧さなく2進浮動小数点数を1対1対応 させられるってこと。
  ただ、それはちゃんとアルゴリズムを実装しないとだめ。Cのprintfとか strtodはちゃんとしてない。

Tag: Programming