计算几何的题做的真是少之又少。
之前wa以为是精度问题，后来发现是情况没有考虑全。。。
题目链接：
http://acm.hdu.edu.cn/showproblem.php?pid=5572
题意：
给定起点A和方向V，路径中遇到圆柱体会发生折射，问能否到达终点B。
分析：
将路径表示为a+t∗v得到关于t的二元方程组，求出Δ。
Δ小于等于0时，表示不会发生折射。直接判断ab是否共线。
Δ大于0时，求出根。
- 根小于0说明路上不会发生折射，判断ab是否共线。
- 根大于等于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;
}