誤差の話（１）

雑記 小ネタ

プログラミングコンテストで、特に幾何の問題とか解いてると「EPS変えたら通った」みたいなことがよくある。僕のような自分に優しい人間は「まあ本質的な部分は合ってたんだしいいか」と考えてしまいがち。この態度を「甘い」ととる人もいればそうでない人もいるだろう（妥協すべきラインというのはどこかに存在するので）。
ロバストな計算幾何というのは研究分野が存在するくらいなので、本質的に難しいものだと思われる。でも、最低限これくらいは知っといたほうがいいんじゃないの、というのはあると思うので思考をメモしておく。タイトルに（１）とついてるのは後から追加するかもしれない、という意味。
浮動小数点数（IEEE 754）の仕様とかは色々調べるのが面倒なので後で書く。TopCoderの問題では非正規化数のせいでTLEとか出たし、知っておいて損はないと思う。最低限持っておく知識としては、仮数部が52ビットなので、一回の計算毎に相対誤差が2^(-52)位積もっていくということとかだろうか。二分探索の終了条件間違えて無限ループ、とかもこの辺に関係する話。あとはintの範囲なら誤差なしで扱えるとか（long doubleなら64ビット整数）、その中でなら余りのでない割り算が正確にできるとか。
浮動小数点数の比較に==をつかっては駄目、というのはよく知られた事実。で、値が同じかどうか調べるのに

const double  EPS = 1e-9;

bool equal(double a, double b){
    return fabs(a-b) < EPS;
}

的なものを書いてる人とか結構見る。だけど、こういう絶対誤差基準だと値が大きいものに対しては判定が厳しく、値の小さいものに対しては判定が緩くなってしまう。なので、相対誤差ではかった方がいいんじゃないかと思ったり。

const double EPS = 1e-9;

bool relative_equal(double a, double b){
    if(a == 0) return fabs(b) < EPS;
    return fabs((a-b)/a) < EPS;
}

みたいな。実際のところとしては、大抵の幾何問題では与えられる座標の値が絶対値1000以下の整数だったりする場合が多いのでとくに気にする必要はないのだけど。
でも上のような判定をしても情報落ちは避けられない。たとえば

int main(){
	double a = 1.3, b = 1e8, c = a + b;
	c = c - b;
	printf("%d\n", equal(a, c));
	printf("%d\n", relative_equal(a, c));
	return 0;
}

上のコードの実行結果は

0
0

である。これを回避するには相当複雑な手順が必要なので、余計に大きい数を足したり引いたりしたりするのは可能な限り避けるのが一番な気がする。
誤差死を確実に回避する手段は、全て整数のまま処理することだと思う。もちろん、整数のままだと無理な問題も多いけど、できるだけ整数のまま処理する努力とかは必要。
例えば、２円（中心の座標、半径ともに整数）が外接するかどうかとかは、

int sq(int x){
    return x*x;
}

fabs(sqrt(sq(x1-x2) + sq(y1-y2)) - (r1 + r2)) < EPS

とするよりも、

sq(x1-x2) + sq(y1-y2) == sq(r1 + r2)

とする方がよい。
他には、格子点上の２点の中点を出す必要のある問題で、座標をすべて２倍するとかもそう（本当は、２のベキの割り算なので大丈夫な気もするけど）。大抵の場合、浮動小数点を出す必要が出てくるのは割り算とsqrtが原因だったりする。なので、逆の操作であるn倍とn乗を上手く使ってやると、整数のまま話ができて嬉しいということもある。
これ以外にも二分探索の中でEPS使ったら誤差死する、とかの話題もあるけどとりあえず今回はここまでにしておく。