いやな誤差

ぼくは幾何の問題を解くときいつもccw周りは

#include <complex>
using namespace std;
typedef complex<double> P;
const double EPS=1e-8, PI=acos(-1.0);

double dot(P a,P b) { return real(conj(a)*b); }
double cross(P a,P b) { return imag(conj(a)*b); }
int ccw(P a, P b, P c) {
    b -= a; c -= a;
    if (cross(b, c) > +EPS)   return +1;    // counter clockwise
    if (cross(b, c) < -EPS)   return -1;    // clockwise
    if (dot(b, c) < 0)     return +2;       // c--a--b on line
    if (norm(b) < norm(c)) return -2;       // a--b--c on line
    return 0;                               // a--c--b
}

みたいにしていた.cross(b, c)をEPSで評価してるのは誤差を考慮するため.


なんとなくいや〜な感じはしていたのだけど,大体の問題で上手くいってるし別にいいかなぁと思って放置していた.なんとなくいやな感じがするというのは,cross(b,c)は長さに関して2次式で誤差の積み重ねが強いのにそれを右辺が吸収できていない感じがしていたからだ.
具体的には,もし真の値 a=(x1,y1), b=(x2,y2) に適当な誤差が加わって a'=(x1+ε1, y1+δ1), b'=(x2+ε2, y2+δ2) となっていたとすると,これで外績をとったときの誤差は a'×b' - a×b ≒ (ε1*y2+δ2*x1) - (δ1*x2+ε2*y1) のようになって,もとの誤差に長さ(大きい値かもしれない!)が掛かった風になってしまう.これはあまりよくない.


というわけで今日は多分これが原因でWAを喰らってしまったので,修正して次のようにすることにした.absを取っているので若干遅くはなるが,次のものなら誤差にはそこそこ強そう,な気がする.本当に直ったのかよくわからないけどWAがACになったのでとりあえず満足することにしている.小数誤差の類のエラーは本当に直ったのかそれとも駄目なのかが目で確認しずらいのでなかなか扱いにくいもんです.

#include <complex>
using namespace std;
typedef complex<double> P;
const double EPS=1e-8, PI=acos(-1.0);

double dot(P a,P b) { return real(conj(a)*b); }
double cross(P a,P b) { return imag(conj(a)*b); }
int ccw(P a, P b, P c) {
    b -= a; c -= a;
    double len = abs(b)*abs(c)
    if (cross(b, c) > +EPS*len)   return +1;    // counter clockwise
    if (cross(b, c) < -EPS*len)   return -1;    // clockwise
    if (dot(b, c) < 0)     return +2;       // c--a--b on line
    if (norm(b) < norm(c)) return -2;       // a--b--c on line
    return 0;                               // a--c--b
}