2016/05/13

誤差が生じるとき

Ruiさんがおもしろい問題をtweetしてた。もすこし厳密に書き直すとこんな感じ。

IEEE単精度浮動小数点数演算で、 x + 0.25 - 0.25 が x - 0.25 + 0.25 と同じにならない数xを求めよ。

xをひとつだけ求めるなら簡単だけれど、xを全て求めるのは難問になると思う。

実は、浮動小数点数の誤差についてはこれまで何度もはまってきたので、このくらい簡単さ、 と思って舐めてかかったんだけど、 全部のケースを網羅するのはやっぱり難問だった。 浮動小数点数はどこまでも油断ならない。

(1) 0.25が分解能の半分に近い場合

すぐに考えつくのはこれだろう。 浮動小数点数の分解能d(ここでは隣の浮動小数点数までの距離とする)の半分より 小さな数値は足しても引いても元の数に影響を与えない。

0.25がd/2未満なら演算の順番を入れ替えても結局変わらないんじゃないか、 と思うかもしれないが、 浮動小数点数には、プラス側とマイナス側で分解能が違う場合があるのだ。 ちょうど±2^mの場合である。 2^mを境に仮数部が受け持つ範囲が一桁シフトするので、2^mより上では隣との距離が 2^(m-p) (pはhidden bitを含まない仮数部の桁数)になるけれど、2^mより下ではその半分の距離になる。

この絵は昔のを流用してるのでdoubleって書いてあるけどfloatでも同じ。

(厳密に言えば、これは±2^mの両側が正規化数の場合。非正規化数領域ではそれ以上分解能を 細かくできないので両側の分解能は同じになる。 以前これに起因するバグを踏んだ)

プラス側のdが0.5、マイナス側が0.25になるような2^mを選ぶと、

• 2^m + 0.25 は隣の浮動小数点数2^m + 0.5とのちょうど中間に落ちるので、 デフォルトの偶数丸め規則により結果は2^mに丸められる
• 2^m - 0.25 は浮動小数点数として正確に表現可能。

であるので、

• 2^m + 0.25 - 0.25 だと、最初の2^m + 0.25が2^mに丸められちゃうので、 結果は2^m - 0.25の値。
• 2^m - 0.25 + 0.25 だと、最初の2^m - 0.25は正しく計算され、 その次の+0.25も正しく計算されて、結果は 2^m。

と言う具合に差がでる。xが -(2^m) でもプラスとマイナスが逆になるだけで同じ議論。

(なお、丸め設定が切り上げの場合、2^m+0.25は2^m+0.5に切り上げられるが、そこから0.25を引いた結果も2^m+0.5に切り上げられるので、やっぱり同じにならない)

後は、こういうmは何かってことになる。IEEE単精度浮動小数点数は仮数部が23bitあるので、 []内を2進表記とすると:

    [1.00000_00000_00000_00000_000] * 2^m  が x
    [1.00000_00000_00000_00000_001] * 2^m  が x+0.5

となるような値ってことだから、0.5 * 2^23が x、つまり x = 2^22。 符号を逆にした場合も含め、このケースに該当するのは x = ±2^22。

(2) 0.25が分解能に近い場合

類似のケースとして、0.25を足すことが2のべき乗の境界をまたぐというのもある。 2^mのプラス側の分解能が0.25、マイナス側の分解能が0.125だったとしよう。 xを 2^m - 0.125 とすると、これに0.25を足したものは2^m + 0.125だが丸められて2^mとなる。 ここから0.25を引くと 2^m - 0.25 になってしまう。 減算から先にやった場合は丸めが生じない。桁を揃えて書いてみると：

    [1.00000_00000_00000_00000_000] * 2^m   が x + 0.125
    [0.00000_00000_00000_00000_0001] * 2^m  ← これが0.125 = 2^-3

なので、正負も考慮するとxは ±(2^21 - 2^-3)。

(3) 0.25を基準にしてxが分解能の半分に近い場合

0.25は

    [1.00000_00000_00000_00000_000] * 2^(-2)

である。(1)とは逆に、xが小さくて、この数値の分解能の 半分以下であった場合にも、同様の可能性が生じる。例えばxが2^-26だった場合。

    [1.00000_00000_00000_00000_000] * 2^(-2)   0.25
    [0.00000_00000_00000_00000_0001] * 2^(-2)  x

0.25+xは丸められて0.25になるが、0.25-xは有効桁数に収まる。 (わかりやすく0.25-xと表記したが、-(x-0.25)と同じなので議論に影響はない。)

このケースは少々複雑だ。xが2^-26より少しでも大きければ0.25+xは 0.25の次の浮動小数点数に近くなるため、そちらに丸められ、加算を先にやっても 減算を先にやっても結果は同じになる。

しかし、xが2^-26より少しだけ小さい場合、加算は0.25に丸められて無かったことになるが、 減算は0.25よりひとつ下の浮動小数点数に丸められる。たとえばx = 0.875*2^(-26)なら、

    [1.00000_00000_00000_00000_000] * 2^(-2)   0.25
    [0.00000_00000_00000_00000_0000111] * 2^(-2)  x

0.25+xは

    [1.00000_00000_00000_00000_0000111] * 2^(-2) 正確な結果
    [1.00000_00000_00000_00000_000] * 2^(-2)     丸め結果

だが、0.25-xは

    [0.11111_11111_11111_11111_1111001] * 2^(-2) 正確な結果
    [1.11111_11111_11111_11111_111001] * 2^(-3)  正確な結果、正規化
    [1.11111_11111_11111_11111_111] * 2^(-3)     丸め結果

となる。xが次のパターンで#にひとつでも1があれば差が生じる。

    [0.00000_00000_00000_00000_00001######...] * 2^(-2)

したがって、このケースでのxは 2^-27 < x ≦ 2^-26 、およびその符号を 反転したもの。

(4) 0.25の加算で桁数が増えて丸めを生じる場合

xの分解能が0.25に対して十分にあっても、丸めが起きてしまうケースがある。 次の例を見てみよう。

    [1.11111_11111_11111_00000_001] * 2^13  x
    [0.00000_00000_00001_00000_000] * 2^13  0.25

加算すると、次々と繰り上がりがおきて、最上桁がひとつ増え、仮数部の有効桁数を 越えてしまい、最後のビットが失われる。

    [10.00000_00000_00000_00000_001] * 2^13  x + 0.25、正確な結果
    [1.00000_00000_00000_00000_0001] * 2^14  正規化
    [1.00000_00000_00000_00000_000] * 2^14   丸め結果

0.25の減算から先にやった場合は丸めが起こらないので、結果が食い違う。

これが起きるケースはどういったパターンだろうか。

• 0.25による加算の繰り上がりがカスケード的に最上位まで伝搬する
• xのLSBが1であり、0.25のビットより下位桁にある

という条件を満たせば良いので、以下の23通りのパターンがあり得る。 #の部分は0でも1でも構わない。

    [1.#####_#####_#####_#####_##1] * 2^-2   x
    [1.00000_00000_00000_00000_000] * 2^-2   0.25
    
    [1.1####_#####_#####_#####_##1] * 2^-1   x
    [0.10000_00000_00000_00000_000] * 2^-1   0.25
    
    [1.11###_#####_#####_#####_##1] * 2^0    x
    [0.01000_00000_00000_00000_000] * 2^0    0.25
    
    [1.111##_#####_#####_#####_##1] * 2^1    x
    [0.00100_00000_00000_00000_000] * 2^1    0.25
    
    …
    
    [1.11111_11111_11111_11111_111] * 2^20   x
    [0.00000_00000_00000_00000_010] * 2^20   0.25

符号が逆の場合も同様。

(z) NaNはNaNにしてNaNにあらず

あと、イレギュラーな回答として「xがNaNの場合」というのもあり得る。 NaNには何を足し引きしてもNaNだが、NaNとNaNを比較すると必ず偽になるので、 プログラムコード上で数値比較していた場合は「等しくない」となる。 ただ、ビットパターンで比べれば等しいかもしれない。等しくないかもしれない (NaNを構成するビットパターンは複数あるので)。

これは問題の「同じにならない」の解釈次第である。

他にあるかな?

多分これで全部だと思うんだけど、他のパターンを見つけたら教えてください。

なお、大元のテスト課題は、答えがpow(2, n) {n > 0} のnを埋める形だったそうなので、 (1)のパターンのみ考慮すれば良かった模様。

Tag: Programming