誤差にまつわるエトセトラ

今日は、浮動小数点小数の数値計算誤差について書こうと思う。


まず、以下のコードを実行してみましょう。

#define EPS 1e-7

int main() {
    // case 1
    double x = 1e16;
    printf("%lf\n", x + 1);

    // case 2
    double y1 = 123456123456.1234588623046875;
    double y2= 123456123456.1234741210937500;

    printf("%d\n", ABS(y1 - y2) < EPS);
    printf("%d\n", ABS(y1 - y2) < EPS * y1);

    // case 3
    double z1 =  1e-1072;
    double z2 =  -1e-1072;
    printf("%d\n", ABS(z1 - z2) < EPS * z1);

    return 0;
}

結果は以下のようになります。


>10000000000000000.000000
0
>1
>0
>1.000000


まず、case1から見ていきましょう。
10000000000000000.000000に1を足しても
10000000000000001.000000にはなりません。
これは、doubleの有効桁数が2^52 ~10^15だからです。最初の15桁以下は丸められてしまいます。逆に言うと、10^15以下の整数ならdoubleで正確に表すことができます。
この辺の話は昔投稿しているので、こちらをご覧ください。


次に、case2です。
これは、おもしろいです。y1とy2をビット列で表すとそれぞれ以下のようになります。


01000010 00111100 10111110 10001110 11110010 01000000 00011111 10011011
01000010 00111100 10111110 10001110 11110010 01000000 00011111 10011100


実は、この2つの数は、doubleの世界では連続する数値です。それにも関らず、(y1-y2)は1e-5程度です。
これは、大変です。連続する値なので、ほんの少しのエラーで、ある値がもう片方の値になりえます。
数値計算を実施する際は、この２つは、同じ値とみなしてよいでしょう。
しかし、ABS(y1-y2) ＜ eps (eps = 1e-7)

なんてやっちゃうと、この2つの値は、同値とは判定できません。

何がまずかったのでしょう・・。

doubleは浮動小数点なので、小数点の位置は固定ではありません。つまり、絶対的な値を誤差の許容範囲とするには無理があります。

ここでは、

ABS(y1-y2) ＜ eps * y1

とするとよいでしょう。

最後にcase3。これは、やっかいです。case2で

ABS(y1-y2) ＜ eps * y1

を判定条件に用いればよいことを示しましたが、case3では、この判定条件は適切ではありません。

同値とみなすべきz1とz2が同値とみなされません。

符号が異なる2つの微量な数値を比較する場合は、case2の場合は適さないようです。

では、どうすればいいか・・・。

case2とcase3を2つ使って判定しましょう。結構面倒ですが・・。。

競技系プログラミングでは、小数が現れないように分母を他辺に掛けるテクニックが有効のようです。また、比較する値の絶対値が比較的小さいものであれば、絶対評価（case2）だけ十分でしょう。

出典）

http://www.topcoder.com/tc?module=Static&d1=tutorials&d2=integersReals

読んで理解したことを書こうと思ったら、ほぼ和訳みたいになりました（笑）

英語が出来る人は、上の原文を読んだ方が分かりやすいはずです。