【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的代码。

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const double eps = 1e-8;

int cmp(double x){
    return x < -eps ? -1 : (x>eps);
}

double pow2(double x){
    return x * x;
}

double mySqrt(double x){
    return sqrt(max((double)0, 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)==0 && cmp(v1.y - v2.y)==0;
    }

    //+, -, 数乘
    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
", c.x, c.y);
    return (double)2*c - p;
}

//判断点p是否在线段v1,v2上
bool onSeg(Vec& v1, Vec& v2, Vec& p){
    if(cmp(product(p-v1, v2-v1))==0 && cmp(dot(p-v1, p-v2))<=0)
        return true;
    return false;
}

bool calc_C(){
    //将AC表示为关于t的参数方程
    //则与圆的方程联立得到关于t的一元二次方程, a,b,c为一般式的三个系数
    double a = pow2(V.x) + pow2(V.y);
    double b = 2*V.x*(A.x - O.x) + 2*V.y*(A.y - O.y);
    double c = pow2(A.x - O.x) + pow2(A.y - O.y) - pow2(R);
    double delta = pow2(b) - 4*a*c;
    if(cmp(delta) <= 0) return false;
    else{
        double t1 = (-b - mySqrt(delta))/(2*a);
        double t2 = (-b + mySqrt(delta))/(2*a);
        double t;
        if(cmp(t1) >= 0) t = t1;
        if(cmp(t2) >= 0 && 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 = 1; 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 = 1;
            else{
                //printf("has intersection (%lf, %lf)
", C.x, C.y);
                Vec A1 = mirrorPoint(O, C, A);
                  // printf("A: %lf, %lf
", A.x, A.y);
                  // printf("A1: %lf, %lf
", A1.x, A1.y);
                //printf("B1 (%lf, %lf)
", B1.x, B1.y);
                if(A==A1){
                    Vec u = B - O;
                    Vec v = C - O;
                    // printf("OB: %lf %lf
", u.unit().x, u.unit().y);
                    // printf("OC: %lf %lf
", v.unit().x, v.unit().y);
                    if(u.unit() == v.unit()){
                        ans = 1;
                    }else ans = 0;
                }else {
                    Vec u = B - C;
                    Vec v = A1 - C;
                    if(u.unit() == v.unit()){
                        ans = 1;
                    }
                    else ans = 0;
                }
            }
        }else{//运动方向不变，则AB与V同向才可碰到B
            //printf("no intersection
");
            Vec temp = B - A;
            if(temp.unit() == V.unit())
                ans = 1;
            else ans = 0;
        }
        printf("Case #%d: %s
", ca, ans ? "Yes" : "No");
    }
    return 0;
}

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

这个版本开始一直WA，不明原因，开始比对上一版本做修改，最后发现把eps由1e-10改成1e-8就过了。几何题精度真是个问题，应该需要试才能知道。

//按照训练指南白书上模版写的新姿势，更好理解
#include <cstdio>
#include <vector>
#include <cmath>
using namespace std;

const double eps = 1e-8; //1e-10会WA，注意调整精度，过大过小都不行
int dcmp(double x){
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
double mySqrt(double x){
    return sqrt(max((double)0, x));
}

struct Point
{
    double x, y;
    Point(double x=0, double y=0):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 = 2*(a*b + c*d);
    double g = b*b + d*d - C.r*C.r;

    double delta = f*f - 4*e*g;

    if(dcmp(delta) <= 0) return 0;

    double t = (-f - mySqrt(delta)) / (2*e);
    P = L.getPoint(t);
    return 1;
}

bool onRay(Point A, Line L){//点A在射线L(p,v)上，不含端点
    Vector w = A - L.p;
    if(dcmp(cross(w, L.v))==0 && dcmp(dot(w, L.v))>0) return true;
    return false;
}

bool onSeg(Point A, Point B, Point C){//点A在线段BC上
    return dcmp(cross(B-A, C-A))==0 && dcmp(dot(B-A, C-A))<0;
}

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
", D.x, D.y);
    return D + (D - A);
}

int main()
{
    int T;
    int ans = 0;
    double R;
    Point O, A, B;
    Vector V;
    freopen("5572.txt", "r", stdin);
    scanf("%d", &T);
    for(int ca = 1; 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 = 1;
            else{
                Line OC = Line(O, Vector(C.x - O.x, C.y - O.y));
                Point A1 = mirrorPoint(A, OC);
                // printf("%lf, %lf
", C.x, C.y);
                // printf("%lf, %lf
", A1.x, A1.y);
                Line CB = Line(C, Vector(B.x - C.x, B.y - C.y));
                
                if(onRay(A1, CB)){
                     ans = 1;

                }
                else ans = 0;
            }
        }else{
            if(onRay(B, LA)) ans = 1;
            else ans = 0;
        }
        printf("Case #%d: %s
", ca, ans ? "Yes" : "No");
    }

    return 0;
}