【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
", 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)
", C.x, C.y);
127                 Vec A1 = mirrorPoint(O, C, A);
128                   // printf("A: %lf, %lf
", A.x, A.y);
129                   // printf("A1: %lf, %lf
", A1.x, A1.y);
130                 //printf("B1 (%lf, %lf)
", B1.x, B1.y);
131                 if(A==A1){
132                     Vec u = B - O;
133                     Vec v = C - O;
134                     // printf("OB: %lf %lf
", u.unit().x, u.unit().y);
135                     // printf("OC: %lf %lf
", 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
");
150             Vec temp = B - A;
151             if(temp.unit() == V.unit())
152                 ans = 1;
153             else ans = 0;
154         }
155         printf("Case #%d: %s
", ca, ans ? "Yes" : "No");
156     }
157     return 0;
158 }

ver1

然后呢，意识到自己实在太弱了，之前的尝试好比盲人摸象，该找些系统的资料好好学一学，于是学习了《训练指南》白书的计算几何入门部分，作者的讲解深入浅出，代码规范清晰且经得起推敲，真是初学者的福音。（即使一些代码没给出，也都可以利用高中学过的知识推出）。于是然后就有了下面这个版本，完全按照上面流程图的思路实现。

这个版本开始一直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
", 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
", C.x, C.y);
115                 // printf("%lf, %lf
", 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
", ca, ans ? "Yes" : "No");
129     }
130 
131     return 0;
132 }