2015上海区域赛现场赛第5题。
题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=5572
题意：在平面上，已知圆(O, R)，点B、A（均在圆外），向量V。
一个小球从A出发，沿速度向量V匀速前进，若撞到圆则发生弹性碰撞，沿“反射方向”继续匀速前进。问A点能否经过B点。
题目读懂了，把所有情况都考虑全，流程图就出来了，然后直接套模版即可（注意有些功能可剪裁~）
我的第一道计算几何，WA了一个星期啊。。。原因有姿势不对，精度不对，还有输出不对的。
不过借此积累了入门的一些模版，向量运算等等，见代码。
下面这个版本可以说是参考网上群巨的代码拼凑着写出的，不过期间发现OJ上的数据没有覆盖所有情况，导致大家各种姿势过的都有，看得云里雾里。
之后自己一遍遍整理思路，考虑各种情况，把别人的代码自己不太赞同的地方强行改成了自己的思路，于是有了下面这个AC的代码。
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std; const double eps = 1e-; int cmp(double x){
return x < -eps ? - : (x>eps);
} double pow2(double x){
return x * x;
} double mySqrt(double x){
return sqrt(max((double), x));
} struct Vec
{
double x, y;
Vec(){}
Vec(double xx, double yy):x(xx), y(yy){} Vec& operator=(const Vec& v){
x = v.x;
y = v.y;
return *this;
} double norm(){
return sqrt(pow2(x) + pow2(y));
}
//返回单位向量
Vec unit(){
return Vec(x, y) / norm();
}
//判等
friend bool operator==(const Vec& v1, const Vec& v2){
return cmp(v1.x - v2.x)== && cmp(v1.y - v2.y)==;
} //+, -, 数乘
friend Vec operator+(const Vec& v1, const Vec& v2){
return Vec(v1.x + v2.x, v1.y + v2.y);
}
friend Vec operator-(const Vec& v1, const Vec& v2){
return Vec(v1.x - v2.x, v1.y - v2.y);
}
friend Vec operator*(const Vec& v, const double c){
return Vec(c*v.x, c*v.y);
}
friend Vec operator*(const double c, const Vec& v){
return Vec(c*v.x, c*v.y);
}
friend Vec operator/(const Vec& v, const double c){
return Vec(v.x/c, v.y/c);
}
}; int T;
int ans;
Vec O, A, V, B, C, B1;
int R; //点乘
double dot(const Vec v1, const Vec v2){
return v1.x*v2.x + v1.y*v2.y;
}
//叉乘的模长
double product(const Vec v1, const Vec v2){
return v1.x*v2.y - v1.y*v2.x;
} //点p到直线v1,v2的投影
Vec project(Vec& v1, Vec& v2, Vec& p){
Vec v = v2 - v1;
return v1 + v * dot(v, p-v1) / dot(v, v);
}
//点p关于直线v1,v2的对称点
Vec mirrorPoint(Vec& v1, Vec& v2, Vec& p){
Vec c = project(v1, v2, p);
//printf("project: %lf, %lf\n", c.x, c.y);
return (double)*c - p;
} //判断点p是否在线段v1,v2上
bool onSeg(Vec& v1, Vec& v2, Vec& p){
if(cmp(product(p-v1, v2-v1))== && cmp(dot(p-v1, p-v2))<=)
return true;
return false;
} bool calc_C(){
//将AC表示为关于t的参数方程
//则与圆的方程联立得到关于t的一元二次方程, a,b,c为一般式的三个系数
double a = pow2(V.x) + pow2(V.y);
double b = *V.x*(A.x - O.x) + *V.y*(A.y - O.y);
double c = pow2(A.x - O.x) + pow2(A.y - O.y) - pow2(R);
double delta = pow2(b) - *a*c;
if(cmp(delta) <= ) return false;
else{
double t1 = (-b - mySqrt(delta))/(*a);
double t2 = (-b + mySqrt(delta))/(*a);
double t;
if(cmp(t1) >= ) t = t1;
if(cmp(t2) >= && t2 < t1) t = t2;//取较小的那个，即为离A近的那个交点
C.x = A.x + V.x*t;
C.y = A.y + V.y*t;
return true; //有交点
}
} int main()
{
freopen("5572.txt", "r", stdin);
scanf("%d", &T);
for(int ca = ; ca <= T; ca++){
scanf("%lf%lf%d", &O.x, &O.y, &R);
scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y);
scanf("%lf%lf", &B.x, &B.y);
if(calc_C()){
if(onSeg(A, C, B)) ans = ;
else{
//printf("has intersection (%lf, %lf)\n", C.x, C.y);
Vec A1 = mirrorPoint(O, C, A);
// printf("A: %lf, %lf\n", A.x, A.y);
// printf("A1: %lf, %lf\n", A1.x, A1.y);
//printf("B1 (%lf, %lf)\n", B1.x, B1.y);
if(A==A1){
Vec u = B - O;
Vec v = C - O;
// printf("OB: %lf %lf\n", u.unit().x, u.unit().y);
// printf("OC: %lf %lf\n", v.unit().x, v.unit().y);
if(u.unit() == v.unit()){
ans = ;
}else ans = ;
}else {
Vec u = B - C;
Vec v = A1 - C;
if(u.unit() == v.unit()){
ans = ;
}
else ans = ;
}
}
}else{//运动方向不变，则AB与V同向才可碰到B
//printf("no intersection\n");
Vec temp = B - A;
if(temp.unit() == V.unit())
ans = ;
else ans = ;
}
printf("Case #%d: %s\n", ca, ans ? "Yes" : "No");
}
return ;
}
ver1
然后呢，意识到自己实在太弱了，之前的尝试好比盲人摸象，该找些系统的资料好好学一学，于是学习了《训练指南》白书的计算几何入门部分，作者的讲解深入浅出，代码规范清晰且经得起推敲，真是初学者的福音。（即使一些代码没给出，也都可以利用高中学过的知识推出）。于是然后就有了下面这个版本，完全按照上面流程图的思路实现。
这个版本开始一直WA，不明原因，开始比对上一版本做修改，最后发现把eps由1e-10改成1e-8就过了。几何题精度真是个问题，应该需要试才能知道。
//按照训练指南白书上模版写的新姿势，更好理解
#include <cstdio>
#include <vector>
#include <cmath>
using namespace std; const double eps = 1e-; //1e-10会WA，注意调整精度，过大过小都不行
int dcmp(double x){
if(fabs(x) < eps) return ;
return x < ? - : ;
}
double mySqrt(double x){
return sqrt(max((double), x));
} struct Point
{
double x, y;
Point(double x=, double y=):x(x), y(y){}
Point& operator = (Point p){
x = p.x;
y = p.y;
return *this;
}
}; typedef Point Vector; Vector operator + (Vector A, Vector B){ return Vector(A.x + B.x, A.y + B.y);}
Vector operator - (Point A, Point B){ return Vector(A.x - B.x, A.y - B.y);}
Vector operator * (Vector A, double p){ return Vector(A.x * p, A.y * p);}
Vector operator / (Vector A, double p){ return Vector(A.x / p, A.y / p);} double dot(Vector A, Vector B){ return A.x * B.x + A.y * B.y;}
double length(Vector A){ return mySqrt(dot(A, A));}
double cross(Vector A, Vector B){ return A.x * B.y - A.y * B.x;} struct Line
{
Point p;
Vector v;
Line(Point p, Vector v):p(p), v(v){}
Point getPoint(double t){
return Point(p.x + v.x*t, p.y + v.y*t);
}
}; struct Circle
{
Point c;
double r;
Circle(Point c, double r):c(c), r(r){}
}; int getLineCircleIntersection(Line L, Circle C, Point& P){ //返回t较小的那个点
double a = L.v.x;
double b = L.p.x - C.c.x;
double c = L.v.y;
double d = L.p.y - C.c.y; double e = a*a + c*c;
double f = *(a*b + c*d);
double g = b*b + d*d - C.r*C.r; double delta = f*f - *e*g; if(dcmp(delta) <= ) return ; double t = (-f - mySqrt(delta)) / (*e);
P = L.getPoint(t);
return ;
} bool onRay(Point A, Line L){//点A在射线L(p,v)上，不含端点
Vector w = A - L.p;
if(dcmp(cross(w, L.v))== && dcmp(dot(w, L.v))>) return true;
return false;
} bool onSeg(Point A, Point B, Point C){//点A在线段BC上
return dcmp(cross(B-A, C-A))== && dcmp(dot(B-A, C-A))<;
} Point project(Point A, Line L){
return L.p + L.v * (dot(L.v, A - L.p)/dot(L.v, L.v));
}
Point mirrorPoint(Point A, Line L){
Vector D = project(A, L);
//printf("project: %lf, %lf\n", D.x, D.y);
return D + (D - A);
} int main()
{
int T;
int ans = ;
double R;
Point O, A, B;
Vector V;
freopen("5572.txt", "r", stdin);
scanf("%d", &T);
for(int ca = ; ca <= T; ca++){
scanf("%lf%lf%lf", &O.x, &O.y, &R);
scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y);
scanf("%lf%lf", &B.x, &B.y);
Line LA = Line(A, V);
Circle yuanO = Circle(O, R);
Point C;
if(getLineCircleIntersection(LA, yuanO, C)){
if(onSeg(B, A, C)) ans = ;
else{
Line OC = Line(O, Vector(C.x - O.x, C.y - O.y));
Point A1 = mirrorPoint(A, OC);
// printf("%lf, %lf\n", C.x, C.y);
// printf("%lf, %lf\n", A1.x, A1.y);
Line CB = Line(C, Vector(B.x - C.x, B.y - C.y)); if(onRay(A1, CB)){
ans = ; }
else ans = ;
}
}else{
if(onRay(B, LA)) ans = ;
else ans = ;
}
printf("Case #%d: %s\n", ca, ans ? "Yes" : "No");
} return ;
}