HDU5572 An Easy Physics Problem【计算几何】

首先考虑射线A->A+V与圆是否会有交点

如果没有：
判断点B是否在射线A->A+V上。

如果有交点：
求出与A点更近的那个交点pInter
首先判断点B是否在线段（A，pInter）上，如果是则会相遇。
如果不是，继续判断：
求出A点关于直线（O，pInter）对称的点pS，判断点B是否在射线pInter->pS上。

（训练的时候想着不求对称点直接根据角度判断就WA了，现在想想精度真的炸了呀，以后写题尽量都不要用到角度吧，毕竟范围[0,PI]。还有这个代码精度1e-13也WA了，要到1e-8就过了）

代码：

    #include<bits/stdc++.h>
    using namespace std;
    #define sq(x) ((x)*(x))
    const double eps = 1e-8;
    //#define double long double
    int sgn(double x) {
        if (fabs(x) < eps) return 0;
        return x < 0? -1: 1;
    }
    struct point {
        double x, y;
        void output() {
            printf("x%f y%f\n", x, y);
        }
        void input() {
            scanf("%lf%lf", &x, &y);
        }
        point (double x=0., double y=0.): x(x), y(y) {}
        point operator+ (const point &rhs) const {
            return point (x+rhs.x, y+rhs.y);
        }
        point operator- (const point &rhs) const {
            return point (x-rhs.x, y-rhs.y);
        }
        double operator* (const point &rhs) const {
            return x * rhs.x + y * rhs.y;
        }
        double operator^ (const point &rhs) const {
            return x * rhs.y - y * rhs.x;
        }
        point operator* (double k) const {
            return point (x*k, y*k);
        }
        point operator/ (double k) const {
            return point (x/k, y/k);
        }
    };
    point operator* (double a, point b) {
        return point(a*b.x, a*b.y);
    }
    double vlen(point p) {
        return sqrt(p*p);
    }
    point unit(point p) {
        return p*(1.0/vlen(p));
    }
    point lerp(point a, point b, double t) {
        return a * (1-t) + b*t;
    }
    double project(point a, point b) {
        return a * unit(b);
    }
    
    double angleVV(point a, point b) {
        return acos(project(a,b)/vlen(a));
    }
    bool Quadratic(double A, double B, double C, double *t0, double *t1) {
        double discrim=sq(B)-4.f*A*C;
        if (discrim < 0.) return false;
        double rootDiscrim=sqrt(discrim);
        double q;
        if (B<0) q= -.5f * (B - rootDiscrim);
        else     q= -.5f * (B + rootDiscrim);
        *t0 = q/A;
        *t1 = C/q;
        if (*t0 > *t1) swap(*t0, *t1);
        return true;
    }
    bool CL_inter(point o, double r, point a, point b, double *t0, double *t1) {
        point d = b-a;
        double A=sq(d);
        double B=d*(a-o)*2.;
        double C=sq(a-o)-sq(r);
        return Quadratic(A,B,C,t0,t1);
    }
    double slen(point a) {
        return a*a;
    }
    bool sameLine(point a, point b) {
        return !sgn(a*b - vlen(a)*vlen(b));
    }
    bool PointOnS(point p, point a, point b) {
        return !sgn((p-a)^(b-a)) && sgn((p-a)*(p-b))<=0;
    }
    bool twoside(point a, point b, point c) {
        double pl = a^b, pr = a^c;
        return sgn(pl)*sgn(pr)<0 || (!sgn(pl) && !sgn(pr));
    }
    point ABClineInter(double a1, double b1, double c1, double a2, double b2, double c2) {
        return point((c1*b2-c2*b1)/(a2*b1-a1*b2),(c1*a2-c2*a1)/(a1*b2-a2*b1));
    }
    point Symmetry(point p, point a, point b) {
        double A = a.y-b.y, B = b.x-a.x, C = -A*a.x-B*a.y;
        double a1 = A, b1 = B, c1 = A*p.x+B*p.y+2.*C;
        double a2 = -B, b2 = A, c2 = B*p.x-A*p.y;
        return ABClineInter(a1, b1, c1, a2, b2, c2);
    }
    
    int T, ks;
    point a, b, o, v;
    double r;
    int main(){
        scanf("%d",&T);
        while (T--) {
            o.input(); scanf("%lf", &r);
            a.input(); v.input(); b.input();
            printf("Case #%d: ", ++ks);
            double t0,t1;
            if (!CL_inter(o,r,a,a+v,&t0,&t1) || !sgn(t1-t0)) {
                if (sameLine(v,b-a)) {
                    puts("Yes");
                } else puts("No"); continue;
            }
            point pInter = lerp(a,a+v,t0), pInter_ = lerp(a,a+v,t1);
            if (vlen(a-pInter) > vlen(a-pInter_)) swap(pInter, pInter_);
            if (!sameLine(v,pInter-a)) {
                if (sameLine(v,b-a)) {
                    puts("Yes");
                } else puts("No"); continue;
            }
            if (PointOnS(b, a, pInter)) {
                puts("Yes"); continue;
            }
            point pS = Symmetry(a, pInter, o);
            if (sameLine(b-pInter,pS-pInter)) {
                puts("Yes");
            } else puts("No");
        }
    }