浮動小数点数の表現に関する、特徴的な部分や罠にはまりそうな部分の非包括的すぎるメモ。浮動小数点数がおおざっぱに x * (2 ** y) みたいに表現されていることは知っているけど、詳細はよく知らんという向け。
浮動小数点数の符号化方式として標準的な IEEE754 では、
± (1.xxxx) * 2 ** (yyyy - bias) // xxxx, yyyy は二進数
の形で符号、仮数部 xxxx 、指数部 yyyy を符号化する。仮数部の 1 は符号化しないのがポイント。 1-bit 節約できる以上に、仮数部が自然に [1, 2) の範囲に制限され、任意のビット列 xxxx yyyy と浮動小数点数が (だいたい) 1:1 対応するのが気持ちいい。この形で表される数を正規化数と呼ぶ。
ただ、このままでは表現できる値の絶対値に下限ができてしまう。0 も表現できないし、仮に 0 に対応する表現を例外的に設けたとしても、 0 と最小の正規数の間がもの凄ーく開いてしまう。
そこで非正規化数 (denormal number / denormalized number) というものを考える。これは 0 と最小正規化数の間を固定小数点でつなぐもので、指数部が最小のときは例外的に仮数部の解釈を
± (0.xxxx) * 2 ** (最小の指数 + 1)
と変えることで導入される。
刻み幅は以下のようなイメージ。
0 最小正規化数
: :
非正規化数なし -+---+++++-+-+-+-+---+---+---+---+-------+-------+-------+-----
: :
非正規化数あり -+-+-+-+-+-+-+-+-+---+---+---+---+-------+-------+-------+-----
: :
非正規化数の分かりやすいメリットとしては、例えば x == y と x - y == 0.0 が (後述の NaN, Inf を除いて) 等価になる。
もちろん非正規化数が 0 近傍の問題を何でも解決してくれるわけではない。非正規化数は固定小数点表現なので、乗除算については気をつける必要がある。例として、信号処理でよく使われる sinc 関数の実装を考える。
double sinc( double x ) {
if( x == 0.0 )
return 1.0;
else
return sin( x ) / x;
}
この関数は x が非正規化数に突入すると精度的にまずい。つまり、
double sinc( double x ) {
// std::numeric_limits<double>::min() は double の最小正規化数。
if( x < numeric_limits<double>::min() )
return 1.0;
else
return sin( x ) / x;
}
こう実装するのが正しい（ただし sinc(x) = 1 + O(x^2) なので、この場合は 1e-8 くらいで分けても十分だったりする。例が良くない…）。
IEEE754 には正規化数と非正規化数の他にいくつか特殊な数表現がある。
-0, +0, -Inf, +Inf, qNaN, sNaN
Inf は無限大に相当する表現で、 1.0 / 0.0 や結果が正規数の最大値を超えた場合に返される。 -0, +0, -Inf, +Inf に関しては、そこそこ予想される範囲内の挙動をする (と思う) 。 1.0 / (-0) == -Inf とか。なので詳細略。
問題は NaN 。 Not a Number の略。こいつは 0.0 / 0.0 や Inf - Inf など結果が不定な場合に返される他、 sqrt(-1.0) なんかの未定義の演算でも生成される。 NaN の何が怖いかというと
NaN < NaN, NaN == NaN, NaN != NaN, NaN > NaN は全て偽
であること。つまり浮動小数点数 x に関して
(x == x) や (x <= 0.0 || 0.0 < x) のような式は偽に成り得る。
具体的にはまりそうな例として、ある点がマンデルブロ集合に属さないことを判定する次のプログラムを考える。
bool not_in_mandelbrot_set( double x, double y ) {
double re = 0.0, im = 0.0;
double re2, im2;
for( int i = 0; i < N; ++i ) {
re2 = re * re;
im2 = im * im;
im = 2.0 * re * im + y;
re = re2 - im2 + x;
}
return re2 + im2 > 4.0;
}
この関数には NaN に関するバグがあって、発散する点でも多くの場合 false を返してくる。ループの途中で浮動小数点で表現できる最大値を超えてしまっても、 Inf になるだけで特に問題ないんじゃないの、と一瞬思ってしまうのだけど、それは間違い。
re2, im2 が Inf → re2 - im2 が NaN → re2, im2 が NaN → re2 + im2 > 4.0 が偽
なので、最後の return は、
return !(re2 + im2 <= 4.0);
と書かないといけない。
浮動小数点数に関するコンパイラの最適化は、一見しょぼく見えることが多い。これはコンパイラが悪いのではなくて、浮動小数点数の性質によるもの。例えば、精度も含めて考えると、一般に浮動小数点数は結合則を満たさないので、
y = x + 3.0 + 4.0;
これは最適化されない。次のように書いておけば最適化される。
y = x + (3.0 + 4.0);
上述の NaN などの特殊な数も最適化を阻害する。その他
y = x + 0.0;
こんなのも (-0) + (+0) = (+0) なので最適化されない。
実際問題として、結果が NaN や Inf になった時点で例外が飛んでくれた方が嬉しいときもある。どうも標準化された方法はないっぽいのだけれど、 UNIX で比較的ポータブルに使えそうな feenableexcept() という関数がある。
#define FE_INVALID #define FE_DENORMAL #define FE_DIVBYZERO #define FE_OVERFLOW #define FE_UNDERFLOW #define FE_INEXACT int feenableexcept( int );y.fujii <y-fujii at mimosa-pudica.net>