计算几何（ZOJ 3783，Coins Game）

一开始RE，不知道咋回事，后来发现原来是assert(0)搞的（前面写的后面给忘了），说明我没有找到任何一个与移动硬币相切的固定硬币。

当时觉得应该是精度的问题，但是试了各种eps的值依然找不到。

我判定相切的办法是求两个圆交点的个数是不是一。后来觉得没必要这样判定，有更好的办法，那就是直接判定半径和是否等于圆心距，然后就AC了。

经过测试后发现错误原因是相切的时候函数判断有两个交点。理论上来说判断点是否重合的函数应该会更严格才对，所以误差应该是在计算两点坐标的时候就被放大了，一旦涉及到角度的问题一般误差都会被放大很多。这个误差似乎没有什么好的办法可以避免，所以说对于计算几何的题目来说，每个判断都应该在误差没被放大的时候用误差最小的判断方法来判断，否则就可能出现错误。尝试了各种eps的值依然不行，我也不清楚为什么会这样，可能是因为当eps足够宽容的时候，另一些判断又被掩盖了。所以说有时候调整eps是没有用，或者说很难成功的。我们更应该把每个细节都处理得尽量精确。

虽然说被卡了，但是数据还是很水的，就算我去掉了一些增加鲁棒性的代码，也可以AC。如果去掉那段代码的话，我自己是能找到一些简单数据把代码跑成死循环的。

具体思路就是圆心走过的距离除以圆周长就是转过的圈数，乘以360度就是转过的角度数。

然后模拟滚动一遍计算圆心走过的距离就好了。

一些注意的地方，angle(Vector v)函数可以返回向量v与x轴之间的夹角，逆时针为正，顺时针为负，范围是(-PI,PI]。

实现方法是调用atan2()函数，这个函数可以返回方位角（与现实中的方位角定义有所不同，唯一区别是值减少了PI），即向量v与y轴之间的夹角，顺时针为正，逆时针为负，范围是(-PI,PI]，所以要调用atan2(v.y.v.x)才对。

因为我们要顺时针旋转，所以极角要从大到小排序。（其实逆时针模拟也一样，相应的极角就从小到大排序）

为了保证能转一圈找到下一个点，所以每个点都要复制一个极角小2*PI的点。（如果逆时针，就大2*PI）

代码

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1010;
const double PI = acos(-1);

struct Point
{
    double x,y;
    double ang;
    int bel;
    Point(double x=0,double y=0,double ang=0,int bel=0):x(x),y(y),ang(ang),bel(bel){}
    void Read()
    {
        scanf("%lf %lf",&x,&y);
        ang=bel=0;
    }
    bool operator < (const Point& rhs) const
    {
        return ang>rhs.ang;
    }
};

typedef Point Vector;

const double eps = 1e-10;
int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return x>0?1:-1;
}

bool operator == (const Point& a,const Point& b)
{
    return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}

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 t)
{
    return Vector(A.x*t,A.y*t);
}
Vector operator / (Vector A,double t)
{
    return Vector(A.x/t,A.y/t);
}

double Dot(Vector A,Vector B)
{
    return A.x*B.x+A.y*B.y;
}

double Len(Vector A)
{
    return sqrt(Dot(A,A));
}

struct Circle
{
    Point c;
    double r;
    vector<Point>vec;
    Circle(Point c=Point(),double r=0):c(c),r(r){}
    Point point(double a)
    {
        //cout<<c.x+cos(a)*r<<" "<<c.y+sin(a)*r<<endl;
        return Point(c.x+cos(a)*r,c.y+sin(a)*r);
    }
    void Read()
    {
        c.Read();
        scanf("%lf",&r);
        vec.clear();
    }
};

double angle(Vector v){return atan2(v.y,v.x);}

int gcci(Circle C1,Circle C2,vector<Point>&sol)
{
    double d=Len(C1.c-C2.c);
    if(dcmp(d)==0)
    {
        if(dcmp(C1.r-C2.r)==0) return -1;
        return 0;
    }
    if(dcmp(C1.r+C2.r-d)<0) return 0;
    if(dcmp(fabs(C1.r-C2.r)-d)>0) return 0;

    double a = angle(C2.c-C1.c);
    double da = acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
    Point p1 = C1.point(a-da);
    Point p2 = C1.point(a+da);

    sol.push_back(p1);
    if(p1==p2) return 1;
    sol.push_back(p2);
    return 2;
}

int N;
Circle C[maxn];
Circle O,S;
int kkk;

void GetS()
{
    vector<Point>sol;
    for(int i=0;i<N;i++) if(dcmp(O.r+C[i].r-Len(O.c-C[i].c))==0)
    {
        S=C[i];
        kkk=i;
        return;
    }
    assert(0);
}

void pre()
{
    for(int i=0;i<N;i++) C[i].r+=O.r;
    for(int i=0;i<N;i++) for(int j=i+1;j<N;j++)
    {
        vector<Point>sol;
        int num=gcci(C[i],C[j],sol);
        if(num==2)
        {
            //puts("----------");
            //printf("%d %d:\n",i,j);
            //cout<<sol[0].x<<" "<<sol[0].y<<endl;
            //cout<<sol[1].x<<" "<<sol[1].y<<endl;
            //puts("----------");
            sol[0].bel=j;
            sol[1].bel=j;
            C[i].vec.push_back(sol[0]);
            C[i].vec.push_back(sol[1]);
            sol[0].bel=i;
            sol[1].bel=i;
            C[j].vec.push_back(sol[0]);
            C[j].vec.push_back(sol[1]);
        }
        else if(num==1)
        {
            sol[0].bel=j;
            C[i].vec.push_back(sol[0]);
            sol[0].bel=i;
            C[j].vec.push_back(sol[0]);
        }
    }
    O.c.bel=-1;
    C[kkk].vec.push_back(O.c);
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<(int)C[i].vec.size();j++) C[i].vec[j].ang=angle(C[i].vec[j]-C[i].c);
        sort(C[i].vec.begin(),C[i].vec.end());
        int sz=C[i].vec.size();
        for(int j=0;j<sz;j++)
        {
            Point pt = C[i].vec[j];
            pt.ang=C[i].vec[j].ang-2*PI;
            C[i].vec.push_back(pt);
        }
    }
}

int r[maxn<<2];
vector<Point>sb;
int cmp1(int i,int j){return sb[i]<sb[j];}

void solve()
{
    scanf("%d",&N);
    for(int i=0;i<N;i++) C[i].Read();
    O.Read();
    GetS();
    pre();
    double len=0;
    Circle nowc = C[kkk];
    Point nowp = O.c;
    while(1)
    {
        nowp.ang = angle(nowp-nowc.c);
        int k = upper_bound(nowc.vec.begin(),nowc.vec.end(),nowp)-nowc.vec.begin();
        sb.clear();
        for(int i=k;i<(int)nowc.vec.size();i++)
            if(nowc.vec[i]==nowc.vec[k]) sb.push_back(nowc.vec[i]);
            else break;
        if(sb.size()>1)
        {
            for(int i=0;i<(int)sb.size();i++)
            {
                sb[i].ang=angle(C[sb[i].bel].c-nowc.c);
                r[i]=i;
            }
            sort(r,r+sb.size(),cmp1);
            k+=r[0];
        }
        //puts("----------");
        //for(int i=0;i<(int)nowc.vec.size();i++) printf("(%lf %lf) ",nowc.vec[i].x,nowc.vec[i].y);
        //puts("\n----------");

        len+=nowc.r*(nowp.ang-nowc.vec[k].ang);
        nowp=nowc.vec[k];
        if(nowp.bel==-1) break;
        nowc=C[nowp.bel];
    }
    printf("%.10lf\n",len/PI*180/O.r);
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--) solve();
    return 0;
}