HDU 5572 An Easy Physics Problem【计算几何】

计算几何的题做的真是少之又少。
之前wa以为是精度问题，后来发现是情况没有考虑全。。。

---

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=5572

题意：

给定起点A和方向V，路径中遇到圆柱体会发生折射，问能否到达终点B。

分析：

将路径表示为a+t∗v得到关于t的二元方程组，求出Δ。
Δ小于等于0时，表示不会发生折射。直接判断ab是否共线。
Δ大于0时，求出根。

1. 根小于0说明路上不会发生折射，判断ab是否共线。
2. 根大于等于0，仍然要判断ab是否共线，并注意此时b应该a与交点的线段上。如果b未在线段上，则判断折射后能否经过b。

代码：

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>#include<cmath>using namespace std;typedef long long ll;const double eps= 1e-8;const double INF = 1e20;double add(double a, double b){    if(fabs(a + b) < eps * (fabs(a) + fabs(b))) return 0;    else return  a + b;}struct Point{    double x, y;    Point(){}    Point(double x, double y):x(x), y(y){}    Point operator + (Point p){        return Point(add(x, p.x), add(y, p.y));    }    Point operator - (Point p){        return Point(add(x, -p.x), add(y, -p.y));    }    double dot(Point p){        return add(x * p.y,  - y * p.x);    }    Point operator * (double d){        return Point(x * d, y * d);    }};struct Circle{    Point o;    double r;    Circle(){}    Circle(double x, double y, double r):o(x, y), r(r){}};Point a, b, v;Circle c;inline int dcmp(double a){if(fabs(a) < eps) return 0; return a < 0 ? -1:1;}inline bool online(Point a, Point v, Point b){return dcmp(v.dot(b-a)) == 0;}inline double getpos(Point a, Point b){    if(dcmp(a.x) == 0) return b.y / a.y;    else return b.x / a.x;}bool judge(Point a, Point v, Circle c, Point b){    double aa = v.x * v.x + v.y * v.y;    double bb = 2 * v.y * (a.y - c.o.y) + 2 * v.x * (a.x - c.o.x);    double cc = (a.x  - c.o.x) * (a.x - c.o.x) + (a.y - c.o.y) * (a.y - c.o.y) - c.r * c.r;    double delta = bb * bb - 4 * aa * cc;    double t1, t2, t;    if(dcmp(delta) <= 0){        if(online(a, v, b)){            t = getpos(v, b - a);            if(dcmp(t) >= 0) return true;        }    }    double anst = INF;    t1 = (- bb + sqrt(delta))/(2 * aa);    t2 = (- bb - sqrt(delta))/(2 * aa);    if(dcmp(t1) >= 0){        if(dcmp(t2) >= 0) anst = t2;        else anst = t1;        if(online(a, v, b)){            double t = getpos(v, b - a);            if(dcmp(t) >= 0 && t <= anst) return true;        }        Point tmp = a + v * anst;        Point temp = c.o - tmp;        Point revers = Point(-temp.y, temp.x);        double k = temp.dot(tmp - b) / revers.dot(temp);        Point tt = tmp + revers * k;        b = tt * 2 - b;        if(online(a, v, b)){            double tmp = getpos(v, b - a);            if(dcmp(tmp) >= 0) return true;        }    }     if(online(a, v, b)){        double t = getpos(v, b - a);        if(dcmp(t) >= 0) return true;     }     return false;}int main (void){    int T;cin>>T;    int cnt = 1;    while(T--){       double QX, QY, R;        cin>>QX>>QY>>R;        c = Circle(QX, QY, R);        double AX, AY, VX, VY, BX, BY;        cin>>AX>>AY>>VX>>VY>>BX>>BY;        a = Point(AX, AY);        v = Point(VX, VY);        b = Point(BX, BY);        cout<<"Case #"<<cnt<<": ";        if(judge(a, v, c, b)) cout<<"Yes"<<endl;        else cout<<"No"<<endl;        cnt++;    }    return 0;}