丸め誤差の恐怖〜浮動小数点数の扱いに注意する〜
丸め誤差は予想以上に厄介であることを思い知らされた。
こいつのせいで予想外のバグとの格闘に時間を割かれてしまったのだ。
自分への戒めとしてこの文章を残す。
丸め誤差とは?
簡単な例として10/3を考えよう。10を3で割り切ると
3.33333....
と永遠に続き割り切れない。これらの情報すべてを正確に残すには無限の容量が必要なので適当なところで打ち切る必要がある。有効数字10ケタまでを情報として残そう!みたいな。
こうやって有効数字の桁数で切り捨てた時に発生する誤差が丸め誤差である。
ただ、これは10進数の場合であり、コンピュータは2進数なのでちょっと異なる。
例えば0.1は10進数だと誤差なしで書けるけど、2進数で表現すると0.000110011...みたいな循環少数になってしまう。
僕達はふだんついつい10進数脳で物事を考えてしまう上に、丸め誤差もあまり気にせずプログラムを書きがちだけど、これは要注意である。
例
自分がハマった例を書く。
x、yと2つの量(3つの要素をもつベクトル)があったとして
x=(3,3,4)
y=(1,2,7)
とする。ここでx-yの期待値を考えてみよう。xもyも3つの要素をすべて足すと10になるのでx-yの期待値は0になるはずである。実際に手計算してみると
{(3-1)+(3-2)+(4-7)}/3 = 0
まぁ当たり前の結果だ。
次にこれをc言語で書いてみる。
double a = 3/3. + 3/3. + 4/3.;
double b = 1/3. + 2/3. + 7/3.;
cout<<a-b<<endl;
んで出力は0になるかというと思えば
-4.440892e-16
とか出てくるのである。0じゃない。
3で割った時に丸め誤差がab間で発生してしまい、差をとっても0にならなくなってしまうのだ。
「10の-16乗でしょ?ほとんど0なんだから問題ないじゃん」
とか思うかもしれないがそうでもない。
例えばa-bの逆数をとって
c=1/(a-b)
を計算したいとする。
(a-b)が0だとこの量は発散してしまうので、
if((a-b)!=0) c = 1./(a-b);
とか書けばa-bが0じゃないときだけ計算してくれることになるが、先ほどの値を使うと0じゃないので
普通にこのif文をtrueとみなしてしまい、結果とんでもない値がcに代入されてしまう。
解決策
足してから割ればうまくいった。
double a = (3 + 3 + 4)/3.;
double b = (1 + 2 + 7)/3.;
cout<<a-b<<endl;
一度足してしまえば割り算は一回なので、誤差が発生するリスクも少ない。
今回の場合はどちらも10/3.を計算することになるのでぴったり一致する。
誤差の桁の意味
せっかくなのでもう少し踏み込んでみよう。
最初の例だと
a-b=-4.440892e-16
なので、15桁は一致していて、16桁目がずれていることになる。
今回は浮動小数点数のdoubleを使っているが、doubleの場合
(符号1ビット)*(仮数部52ビット)*(指数部11ビット)
の62bit = 8byteからなる。
今回の例だと指数部は問題ないので、大事なのは仮数部の部分で、ここが有効数字の桁数を決めてると考えて良い。
2^52=4.5035996e+15
なので、仮数部は10進数で15桁くらいの精度で計算できることになる。
そう考えると16桁のオーダーで誤差が発生するのもうなずける気がする。
んーと、合ってるかな??
*余談
丸め誤差と似たような量として打ち切り誤差ってのもある。
これサイトによって書いてあることが違って、あるサイ卜だと
「無限級数などの無限に続く計算を途中で打ち切った時に生じる誤差」
っとか書いてある一方で
「循環少数などの計算を途中で打ち切った時に生じる誤差」
って書いてあるページもある。
前者だと足し算の項とかを足すのを途中でやめてるわけだから今回の例とはちょっと異なるし、後者だと今回の無限小数の例も打ち切り誤差ってことになる。
どっちなんだ。手元の本だと前者に近くて、丸め誤差と打ち切り誤差は区別されている。
まぁ無限に続いているものを途中で打ち切ることで生じる誤差、という意味では丸め誤差と似たような量なんだろう。
他にも情報落ち、ケタ落ち、オーバーフローとかいろいろあるけど、今回はこのへんで。
まとめ
・割ってから足すな、足してから割れ
おわりー