HDU 5572-An Easy Physics Problem （计算几何）

题目链接：
http://acm.hdu.edu.cn/showproblem.php?pid=5572

题目大意：
有一个无体积的球，近似为质点，一开始给出球的初始位置和运动方向（为从原点到给出点的向量方向），同时平面中有一个圆柱体，给出圆心位置和半径，若球撞击圆柱体，则产生以撞击点切线为反射面的无损反弹，求问球是否会经过给定的P点

分析：
大体上可以分成两类情况，球会撞击圆柱或者不会撞击圆柱
球什么情况下会撞击圆柱呢

• 球的运动轨迹与圆柱体有交点
• 交点为两个，一个为与圆相切，不发生撞击
• 圆柱体在球的初始位置的前进方向上，而不是后方

如果球未撞击圆柱，只有一种可能，点处于球的运动轨迹前方，即一条射线上，而球撞击圆柱则有两种可能

• 球在初始位置和撞击点（与初始位置较近的一个交点）经过P点
• 球在反弹后，在撞击点的前方射线经过P点

理清楚这几种可能后，直接上计算几何板子，调一下方法即可

代码：

#include<bits/stdc++.h>

using namespace std;

const double PI = acos(-1);

bool ok;
double ox, oy, r;
double px, py;
double vx, vy;
double tx, ty;
double xta1, xta2, xta3, xta4;

const double eps = 1e-7;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 1010;

int sgn(double x){
    if (fabs(x)<eps) return 0;
    if (x<0) return -1;
    else return 1;
}

inline double sqr(double x)
{
    return x*x;
}

struct Point{
    double x,y;
    Point(){}
    Point(double _x,double _y){
        x = _x;
        y = _y;
    }
    void input(){
        scanf("%lf%lf",&x,&y);
    }
    void output(){
        printf("%2.f,%.2f\n",x,y);
    }
    bool operator == (Point b)const{
        return sgn(x-b.x)==0&& sgn(y-b.y)==0;
    }
    bool operator < (Point b)const{
        return sgn(x-b.x)==0?sgn(y-b.y)<0:x<b.x;
    }
    Point operator - (const Point & b)const{
        return Point(x-b.x,y-b.y);
    }
    //叉积
    double operator ^ (const Point &b)const{
        return x*b.y-y*b.x;
    }
    //点积
    double operator * (const Point &b)const{
        return x*b.x + y*b.y;
    }
    double len()
    {
        return hypot(x,y);
    }
    double len2()
    {
        return x*x+y*y;
    }

    double distance(Point p)
    {
        return hypot(x-p.x,y-p.y);
    }

    Point operator + (const Point & b) const
    {
        return Point(x+b.x,y+b.y);
    }

    Point operator *(const double& k)const{
        return Point(x*k,y*k);
    }

    Point operator / (const double &k) const{
        return Point(x/k,y/k);
    }
    double rad(Point a,Point b)
    {
        Point p = *this;
        return fabs(atan2 (fabs((a-p)^(b-p)),(a-p)*(b-p)));
    }

    Point trunc(double r)
    {
        double l = len();
        if (!sgn(l)) return *this;
        r /=l ;
        return Point (x*r,y*r);
    }

    Point rotleft(){
        return Point(-y,x);
    }

    Point rotright(){
        return Point(y,-x);
    }

    Point rotate(Point p,double angle)
    {
        Point v = (*this) - p;
        double c = cos(angle) , s = sin(angle);
        return Point(p.x+v.x*c - v.y*s,p.y+v.x*s + v.y*c);
    }

};

struct Line{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e)
    {
        s = _s;
        e = _e;
    }
    bool operator == (Line v)
    {
        return (s==v.s)&&(e==v.e);
    }
    Line (Point p,double angle)
    {
        s = p;
        if (sgn(angle - pi/2) == 0)
        {
            e = (s + Point(0,1));
        }
        else
            e = (s +Point(1,tan(angle)));
    }

    Line (double a,double b,double c)
    {
        if (sgn(a)==0)
        {
            s = Point(0,-c/b);
            e = Point(1,-c/b);
        }
        else if (sgn(b)==0)
        {
            s = Point(-c/a,0);
            e = Point(-c/a,1);
        }
        else
        {
            s = Point(0,-c/b);
            s = Point(1,(-c-a)/b);
        }
    }
    void input()
    {
        s.input();
        e.input();
    }

    void adjust()
    {
        if (e<s) swap(s,e);
    }

    double length()
    {
        return s.distance(e);
    }

    double angle(){
        double k = atan2(e.y-s.y,e.x-s.x);
        if (sgn(k)<0) k+= pi;
        if (sgn(k-pi)==0) k -= pi;
        return k;
    }

    int relation(Point p)
    {
        int c = sgn((p-s)^(e-s));
        if (c<0) return 1;
        else if (c>0) return 2;
        else return 3;
    }

    bool pointonseg(Point p)
    {
        return sgn((p-s)^(e-s)) ==0 &&sgn((p-s)*(p-e))<=0;
    }

    bool parallel(Line v)
    {
        return sgn((e-s)^(v.e-v.s)) ==0;
    }

    int segcrossseg(Line v)
    {
        int d1 = sgn((e-s)^(v.s-s));
        int d2 = sgn((e-s)^(v.e-s));
        int d3 = sgn((v.e-v.s)^(s-v.s));
        int d4 = sgn((v.e-v.s)^(e-v.s));
        if ((d1^d2) ==-2 && (d3^d4)==-2 ) return 2;
        return (d1==0&& sgn((v.s-s)*(v.s-e))<=0)||
        (d2 ==0 && sgn((v.e-s)*(v.e-e))<=0)||
        (d3 ==0 && sgn((s-v.s)*(s-v.e))<=0)||
        (d4 ==0 && sgn((e-v.s)*(e-v.e))<=0);
    }

    Point crosspoint(Line v)
    {
        double a1 = (v.e-v.s)^(s-v.s);
        double a2 = (v.e-v.s)^(e-v.s);
        return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
    }

    double dispointtoline(Point p)
    {
        return fabs((p-s)^(e-s))/length();
    }

    double dispointtoseg(Point p)
    {
        if (sgn((p-s)*(e-s))<0|| sgn((p-e)*(s-e))<0)
            return min(p.distance(s),p.distance(e));
        return dispointtoline(p);
    }

    Point symmetrypoint(Point p)
    {
        Point q = lineprog(p);
        return Point(2*q.x-p.x,2*q.y-p.y);
    }

    Point lineprog(Point p)
    {
        return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()));
    }

};


struct circle{
    Point p;
    double r;
    circle() {}
    circle(Point _p , double _r)
    {
        p = _p;
        r = _r;
    }

    circle(Point a, Point b, Point c)
    {
        Line u = Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
        Line v = Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
        p = u.crosspoint(v);
        r = p.distance(a);
    }


    void input()
    {
        p.input();
        scanf("%lf",&r);
    }

    void output()
    {
        printf("%.2f,%.2f,%.2f\n",p.x,p.y,r);
    }

    bool operator == (circle v)
    {
        return (p==v.p)&&sgn(r-v.r)==0;
    }
    bool operator < (circle v)const
    {
        return ((p<v.p)||((p==v.p)&&sgn(r-v.r)<0));
    }

    double area()
    {
        return pi*r*r;
    }

    double circumference()
    {
        return 2*pi*r;
    }

    int relation(Point b)
    {
        double dst = b.distance(p);
        if (sgn(dst-r)<0) return 2;
        else if (sgn(dst-r)==0) return 1;
        return 0;
    }

    int relationseg(Line v)
    {
        double dst = v.dispointtoseg(p);
        if (sgn(dst-r) < 0) return 2;
        else if (sgn(dst -r) ==0) return 1;
        return 0;
    }

    int relationline (Line v)
    {
        double dst = v.dispointtoline(p);
        if (sgn(dst-r)<0) return 2;
        else if (sgn(dst-r)==0) return 1;
        return 0;
    }

    int relationcircle(circle v)
    {
        double d = p.distance(v.p);
        if (sgn(d-r-v.r)>0) return 5;
        if (sgn(d-r-v.r)==0) return 4;
        double l = fabs(r-v.r);
        if (sgn(d-r-v.r)<0 && sgn(d-l)>0) return 3;
        if (sgn(d-l) ==0 ) return 2;
        if (sgn(d-l)< 0) return 1;
    }

    int pointcrosscircle(circle v,Point& p1,Point& p2)
    {
        int rel = relationcircle(v);
        if (rel==1 || rel==5) return 0;
        double d = p.distance(v.p);
        double l = (d*d+r*r-v.r*v.r)/(2*d);
        double h = sqrt(r*r-l*l);
        Point tmp = p+ (v.p-p).trunc(l);
        p1 = tmp + ((v.p-p).rotleft().trunc(h));
        p2 = tmp + ((v.p-p).rotright().trunc(h));
        if (rel==2 || rel ==4)
            return 1;
        return 2;
    }

    int pointcorssline(Line v, Point &p1 , Point &p2)
    {
        if (!(*this).relationline(v))   return 0;
        Point a = v.lineprog(p);
        double d = v.dispointtoline(p);
        d = sqrt(r*r-d*d);
        if (sgn(d)==0)
        {
            p1 = a;
            p2 = a;
            return 1;
        }
        p1 = a + (v.e-v.s).trunc(d);
        p2 = a - (v.e-v.s).trunc(d);
        return 2;

    }


};

int main(){
    int T;
    scanf("%d", &T);
    for(int t = 1; t <= T; t++){
        circle cir;
        cir.input();
        Point points,pointv,pointt;
        points.input();
        pointv.input();
        vx = pointv.x;
        vy = pointv.y;
        pointt.input();

        Point pointe(points.x+vx,points.y+vy);

        ok = false;

        Line linep(points,pointe);  //构造直线


        bool flag = false;
        Point p,q;
        int res = cir.pointcorssline(linep,p,q);
        if (sgn(px-tx)==0&&sgn(py-ty)==0)
            flag = true;
        else if (res==2 && sgn(p.x - px) == sgn(vx) && sgn(p.y - py) == sgn(vy))   // 直线与圆相交
        {

            double dis1 = p.distance(points);
            double dis2 = q.distance(points);
            if(dis1 > dis2) swap(p,q);
            Line tpl(points,p);
            if (tpl.pointonseg(pointt))
            {
                flag = true;
            }
            else
            {
                Point symp = Line(cir.p,p).symmetrypoint(points);
                Line reflect(symp,p);
                double nvx = symp.x - p.x , nvy = symp.y - p.y;
                int nrels = reflect.relation(pointt);
                if (nrels==3&&sgn(pointt.x-p.x)==sgn(nvx)&&sgn(pointt.y-p.y)==sgn(nvy))
                {
                    flag = true;
                }
                else
                    flag = false;
            }

        }
        else   //相切 相离 或者未在前进方向上相交
        {
            int rels = linep.relation(pointt);
            if (rels==3&&sgn(pointt.x-points.x)==sgn(vx)&&sgn(pointt.y-points.y)==sgn(vy))// 点在直线上
            {
                flag = true;
            }
            else
                flag = false;

        }
        printf("Case #%d: ", t);
        if(flag) puts("Yes");
        else puts("No");
    }
}