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;
}