【HDU 5572 An Easy Physics Problem】计算几何基础
2015上海区域赛现场赛第5题。
题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=5572
题意：在平面上，已知圆(O, R)，点B、A（均在圆外），向量V。
一个小球从A出发，沿速度向量V匀速前进，若撞到圆则发生弹性碰撞，沿“反射方向”继续匀速前进。问A点能否经过B点。
题目读懂了，把所有情况都考虑全，流程图就出来了，然后直接套模版即可（注意有些功能可剪裁~）
我的第一道计算几何，WA了一个星期啊。。。原因有姿势不对，精度不对，还有输出不对的。
不过借此积累了入门的一些模版，向量运算等等，见代码。
下面这个版本可以说是参考网上群巨的代码拼凑着写出的，不过期间发现OJ上的数据没有覆盖所有情况，导致大家各种姿势过的都有，看得云里雾里。
之后自己一遍遍整理思路，考虑各种情况，把别人的代码自己不太赞同的地方强行改成了自己的思路，于是有了下面这个AC的代码。
1 #include <cstdio> 2 #include <cmath> 3 #include <algorithm> 4 using namespace std; 5 6 const double eps = 1e-8; 7 8 int cmp(double x){ 9 return x < -eps ? -1 : (x>eps); 10 } 11 12 double pow2(double x){ 13 return x * x; 14 } 15 16 double mySqrt(double x){ 17 return sqrt(max((double)0, x)); 18 } 19 20 struct Vec 21 { 22 double x, y; 23 Vec(){} 24 Vec(double xx, double yy):x(xx), y(yy){} 25 26 Vec& operator=(const Vec& v){ 27 x = v.x; 28 y = v.y; 29 return *this; 30 } 31 32 double norm(){ 33 return sqrt(pow2(x) + pow2(y)); 34 } 35 //返回单位向量 36 Vec unit(){ 37 return Vec(x, y) / norm(); 38 } 39 //判等 40 friend bool operator==(const Vec& v1, const Vec& v2){ 41 return cmp(v1.x - v2.x)==0 && cmp(v1.y - v2.y)==0; 42 } 43 44 //+, -, 数乘 45 friend Vec operator+(const Vec& v1, const Vec& v2){ 46 return Vec(v1.x + v2.x, v1.y + v2.y); 47 } 48 friend Vec operator-(const Vec& v1, const Vec& v2){ 49 return Vec(v1.x - v2.x, v1.y - v2.y); 50 } 51 friend Vec operator*(const Vec& v, const double c){ 52 return Vec(c*v.x, c*v.y); 53 } 54 friend Vec operator*(const double c, const Vec& v){ 55 return Vec(c*v.x, c*v.y); 56 } 57 friend Vec operator/(const Vec& v, const double c){ 58 return Vec(v.x/c, v.y/c); 59 } 60 }; 61 62 int T; 63 int ans; 64 Vec O, A, V, B, C, B1; 65 int R; 66 67 //点乘 68 double dot(const Vec v1, const Vec v2){ 69 return v1.x*v2.x + v1.y*v2.y; 70 } 71 //叉乘的模长 72 double product(const Vec v1, const Vec v2){ 73 return v1.x*v2.y - v1.y*v2.x; 74 } 75 76 //点p到直线v1,v2的投影 77 Vec project(Vec& v1, Vec& v2, Vec& p){ 78 Vec v = v2 - v1; 79 return v1 + v * dot(v, p-v1) / dot(v, v); 80 } 81 //点p关于直线v1,v2的对称点 82 Vec mirrorPoint(Vec& v1, Vec& v2, Vec& p){ 83 Vec c = project(v1, v2, p); 84 //printf("project: %lf, %lf\n", c.x, c.y); 85 return (double)2*c - p; 86 } 87 88 //判断点p是否在线段v1,v2上 89 bool onSeg(Vec& v1, Vec& v2, Vec& p){ 90 if(cmp(product(p-v1, v2-v1))==0 && cmp(dot(p-v1, p-v2))<=0) 91 return true; 92 return false; 93 } 94 95 bool calc_C(){ 96 //将AC表示为关于t的参数方程 97 //则与圆的方程联立得到关于t的一元二次方程, a,b,c为一般式的三个系数 98 double a = pow2(V.x) + pow2(V.y); 99 double b = 2*V.x*(A.x - O.x) + 2*V.y*(A.y - O.y); 100 double c = pow2(A.x - O.x) + pow2(A.y - O.y) - pow2(R); 101 double delta = pow2(b) - 4*a*c; 102 if(cmp(delta) <= 0) return false; 103 else{ 104 double t1 = (-b - mySqrt(delta))/(2*a); 105 double t2 = (-b + mySqrt(delta))/(2*a); 106 double t; 107 if(cmp(t1) >= 0) t = t1; 108 if(cmp(t2) >= 0 && t2 < t1) t = t2;//取较小的那个，即为离A近的那个交点 109 C.x = A.x + V.x*t; 110 C.y = A.y + V.y*t; 111 return true; //有交点 112 } 113 } 114 115 int main() 116 { 117 freopen("5572.txt", "r", stdin); 118 scanf("%d", &T); 119 for(int ca = 1; ca <= T; ca++){ 120 scanf("%lf%lf%d", &O.x, &O.y, &R); 121 scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y); 122 scanf("%lf%lf", &B.x, &B.y); 123 if(calc_C()){ 124 if(onSeg(A, C, B)) ans = 1; 125 else{ 126 //printf("has intersection (%lf, %lf)\n", C.x, C.y); 127 Vec A1 = mirrorPoint(O, C, A); 128 // printf("A: %lf, %lf\n", A.x, A.y); 129 // printf("A1: %lf, %lf\n", A1.x, A1.y); 130 //printf("B1 (%lf, %lf)\n", B1.x, B1.y); 131 if(A==A1){ 132 Vec u = B - O; 133 Vec v = C - O; 134 // printf("OB: %lf %lf\n", u.unit().x, u.unit().y); 135 // printf("OC: %lf %lf\n", v.unit().x, v.unit().y); 136 if(u.unit() == v.unit()){ 137 ans = 1; 138 }else ans = 0; 139 }else { 140 Vec u = B - C; 141 Vec v = A1 - C; 142 if(u.unit() == v.unit()){ 143 ans = 1; 144 } 145 else ans = 0; 146 } 147 } 148 }else{//运动方向不变，则AB与V同向才可碰到B 149 //printf("no intersection\n"); 150 Vec temp = B - A; 151 if(temp.unit() == V.unit()) 152 ans = 1; 153 else ans = 0; 154 } 155 printf("Case #%d: %s\n", ca, ans ? "Yes" : "No"); 156 } 157 return 0; 158 }
然后呢，意识到自己实在太弱了，之前的尝试好比盲人摸象，该找些系统的资料好好学一学，于是学习了《训练指南》白书的计算几何入门部分，作者的讲解深入浅出，代码规范清晰且经得起推敲，真是初学者的福音。（即使一些代码没给出，也都可以利用高中学过的知识推出）。于是然后就有了下面这个版本，完全按照上面流程图的思路实现。
这个版本开始一直WA，不明原因，开始比对上一版本做修改，最后发现把eps由1e-10改成1e-8就过了。几何题精度真是个问题，应该需要试才能知道。
1 //按照训练指南白书上模版写的新姿势，更好理解 2 #include <cstdio> 3 #include <vector> 4 #include <cmath> 5 using namespace std; 6 7 const double eps = 1e-8; //1e-10会WA，注意调整精度，过大过小都不行 8 int dcmp(double x){ 9 if(fabs(x) < eps) return 0; 10 return x < 0 ? -1 : 1; 11 } 12 double mySqrt(double x){ 13 return sqrt(max((double)0, x)); 14 } 15 16 struct Point 17 { 18 double x, y; 19 Point(double x=0, double y=0):x(x), y(y){} 20 Point& operator = (Point p){ 21 x = p.x; 22 y = p.y; 23 return *this; 24 } 25 }; 26 27 typedef Point Vector; 28 29 Vector operator + (Vector A, Vector B){ return Vector(A.x + B.x, A.y + B.y);} 30 Vector operator - (Point A, Point B){ return Vector(A.x - B.x, A.y - B.y);} 31 Vector operator * (Vector A, double p){ return Vector(A.x * p, A.y * p);} 32 Vector operator / (Vector A, double p){ return Vector(A.x / p, A.y / p);} 33 34 double dot(Vector A, Vector B){ return A.x * B.x + A.y * B.y;} 35 double length(Vector A){ return mySqrt(dot(A, A));} 36 double cross(Vector A, Vector B){ return A.x * B.y - A.y * B.x;} 37 38 struct Line 39 { 40 Point p; 41 Vector v; 42 Line(Point p, Vector v):p(p), v(v){} 43 Point getPoint(double t){ 44 return Point(p.x + v.x*t, p.y + v.y*t); 45 } 46 }; 47 48 struct Circle 49 { 50 Point c; 51 double r; 52 Circle(Point c, double r):c(c), r(r){} 53 }; 54 55 int getLineCircleIntersection(Line L, Circle C, Point& P){ //返回t较小的那个点 56 double a = L.v.x; 57 double b = L.p.x - C.c.x; 58 double c = L.v.y; 59 double d = L.p.y - C.c.y; 60 61 double e = a*a + c*c; 62 double f = 2*(a*b + c*d); 63 double g = b*b + d*d - C.r*C.r; 64 65 double delta = f*f - 4*e*g; 66 67 if(dcmp(delta) <= 0) return 0; 68 69 double t = (-f - mySqrt(delta)) / (2*e); 70 P = L.getPoint(t); 71 return 1; 72 } 73 74 bool onRay(Point A, Line L){//点A在射线L(p,v)上，不含端点 75 Vector w = A - L.p; 76 if(dcmp(cross(w, L.v))==0 && dcmp(dot(w, L.v))>0) return true; 77 return false; 78 } 79 80 bool onSeg(Point A, Point B, Point C){//点A在线段BC上 81 return dcmp(cross(B-A, C-A))==0 && dcmp(dot(B-A, C-A))<0; 82 } 83 84 Point project(Point A, Line L){ 85 return L.p + L.v * (dot(L.v, A - L.p)/dot(L.v, L.v)); 86 } 87 Point mirrorPoint(Point A, Line L){ 88 Vector D = project(A, L); 89 //printf("project: %lf, %lf\n", D.x, D.y); 90 return D + (D - A); 91 } 92 93 int main() 94 { 95 int T; 96 int ans = 0; 97 double R; 98 Point O, A, B; 99 Vector V; 100 freopen("5572.txt", "r", stdin); 101 scanf("%d", &T); 102 for(int ca = 1; ca <= T; ca++){ 103 scanf("%lf%lf%lf", &O.x, &O.y, &R); 104 scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y); 105 scanf("%lf%lf", &B.x, &B.y); 106 Line LA = Line(A, V); 107 Circle yuanO = Circle(O, R); 108 Point C; 109 if(getLineCircleIntersection(LA, yuanO, C)){ 110 if(onSeg(B, A, C)) ans = 1; 111 else{ 112 Line OC = Line(O, Vector(C.x - O.x, C.y - O.y)); 113 Point A1 = mirrorPoint(A, OC); 114 // printf("%lf, %lf\n", C.x, C.y); 115 // printf("%lf, %lf\n", A1.x, A1.y); 116 Line CB = Line(C, Vector(B.x - C.x, B.y - C.y)); 117 118 if(onRay(A1, CB)){ 119 ans = 1; 120 121 } 122 else ans = 0; 123 } 124 }else{ 125 if(onRay(B, LA)) ans = 1; 126 else ans = 0; 127 } 128 printf("Case #%d: %s\n", ca, ans ? "Yes" : "No"); 129 } 130 131 return 0; 132 }
