sgu 414 Orthogonal Circles 圆的正交

题意：给定n个圆，问能否画出一个圆，使得这个圆与所有圆正交。若不存在输出-1，若个数大于一个就输出-2，若刚好一个则输出圆心坐标以及半径。正交的意思是两个圆相交，且两个圆的圆心与交点的连线正好是另一个圆的切线，即两圆心与交点的连线互相垂直。

题解：

拿两个圆来说存在公式r1^2+r^2=|OA|^2,r2^2+r^2=|OB|^2,如下图，其中A,B是圆心，r1和r2是半径，O是所求正交圆的圆心，r是半径。两式化简下，r1^2-r2^2=d1^2-d2^2。过O作AB的垂线交AB于P，由于令|AP|=a,|BP|=b,|OP|=c，则d1^2-d2^2=(c^2+a^2)-(c^2+b^2)=a^2-b^2=r1^2-r2^2,是一个固定值，所以所有符合条件的正交圆圆心都在直线OP上，我们称这条线为正交线。

当A和B相交时，OP直线上在圆内的点就不符合，因为在园内的圆心不可能作这个圆的正交圆（正交圆的圆心都是一个圆上的两条切线的交点，圆不可能使切线交点在圆内）。

当A和B内含时，也有类似的情况，只不过P点是在AB的延长线。

特殊情况：两圆重合的时候，我们可以只计算一个圆即可。当圆同心（这里的同心不包括重合）的时候，不可能存在正交圆，因为d1^2-d2^2=r1^2-r2^2=0和d1=d2矛盾。

之后推到n个圆，先要排除掉同心圆的情况，所以需要排序下。然后找出各个圆之间的正交线（如果把重合的圆算在内，也只要n-1条，因为正交圆的性质可以传递，ABC三个圆，如果AB有点正交，BC有点正交，两个的正交线有交点，那么这个交点也符合AC的情况），遍历所有的正交线，所有的线都重合，那么线上的所有点都是正交点，输出-2；如果有平行的线出现，或交点不唯一，那么说明没有一个符合点存在，输出-1；如果所有线交一个点，那么这个点符合所有圆的正交，在输出结果前还需要判断下点是否在圆内。

注意：这道题死在了浮点精度上面，一开始改错误，过的测试数据一个个多起来，最后我卡在了第7组上。看了n遍之后还是无解。。最后试了下改大浮点误差eps，竟然过了。这样的浮点问题很少遇到，不清楚怎么改，浮点运算太多了。。哎

代码：

#include <cstdio>#include <cstdlib>#include <cmath>#include <cstring>#include <iostream>#include <algorithm>#include <map>#include <set>#include <queue>using namespace std;//基础点和向量运算struct Point{    double x,y;    Point(double x=0,double y=0):x(x),y(y){}};typedef Point Vector;Vector operator + (Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}Vector operator - (Vector A,Vector 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);}bool operator <(const Point& a, const Point& b){    return a.x<b.x||(a.x==b.x&&a.y<b.y);}const double eps=1e-4;int dcmp(double x)//判断正负，或者等于0{    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;}double Dot(Vector A, Vector B){return A.x*B.x+A.y*B.y;}//点积double Length(Vector A){return sqrt(Dot(A,A));}//OA长double Angle(Vector A,Vector B){return acos(Dot(A,B)/Length(A)/Length(B));}//OA和OB的夹角double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}//叉积Vector Rotate(Vector A,double rad)//rad为弧度，旋转rad度{    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));}//点和直线//P+tv表示一条直线，P为点，tv为方向向量Point GetLineIntersection(Point P,Vector v,Point Q,Vector w)//求直线交点，确保存在交点，即Cross(v,w)非0{    Vector u=P-Q;    double t=Cross(w,u)/Cross(v,w);    return P+v*t;}//圆struct Circle{//定义圆    Point c;    double r;    Circle(){}    Circle(Point c,double r):c(c),r(r){}    Point point(double a){//根据圆心角计算圆上的坐标        return Point(c.x+cos(a)*r,c.y+sin(a)*r);    }};struct Line{//定义线    Point p,v;    Line(){}    Line(Point p,Point v):p(p),v(v){}    Point point(double a){        return p+(v-p)*a;    }};const int maxn=1e5+10;const double pi=atan(1.0)*4;int vis[maxn];Circle cir[maxn];Point p[maxn];Line l[maxn];Line get_line(Circle a,Circle b)//获得正交线{    Line AB(a.c,b.c);    Point P,Q;    Vector R;    double d=Length(a.c-b.c);    //printf("%f\n",d);    double t=(d*d+a.r*a.r-b.r*b.r)/(2*d*d);    //printf("%f %f**\n",d,t);    P=AB.point(t);    R=a.c-b.c;    Q=Point(-R.y,R.x)+P;    //printf("P:%f %f\n",P.x,P.y);    return Line(P,Q);}bool judge(Point P,int n)//判点是否在圆内{    int i;    for(i=0;i<n;i++)        if(dcmp(Length(cir[i].c-P)-cir[i].r)<=0)return false;    return true;}int cmp(Circle a,Circle b)//排序比较{    if(a.c==b.c)return a.r<b.r;    else if(dcmp(a.c.x-b.c.x)==0)return a.c.y<b.c.y;    else return a.c.x<b.c.x;}int main(){    int n;    while(scanf("%d",&n)!=EOF)    {        double a,b,c;        int i,j,k;        for(i=0;i<n;i++)        {            scanf("%lf%lf%lf",&a,&b,&c);            cir[i]=Circle(Point(a,b),c);        }        if(n==1){printf("-2\n");continue;}        int m=0,t=0;        sort(cir,cir+n,cmp);//排序，用于查找重合以及同心圆        //找出所有正交线        for(i=0;i<n-1;i++)        {            if(cir[i].c==cir[i+1].c&&dcmp(cir[i].r-cir[i+1].r)==0)continue;//重合的圆可以跳过            else if(cir[i].c==cir[i+1].c)break;//同心圆（非重合）就输出-1            l[m++]=get_line(cir[i],cir[i+1]);        }        if(i<n-1){printf("-1\n");continue;}        //找出所有正交线的交点        for(i=0;i<m-1;i++)        {            if(dcmp(Cross(l[i].p-l[i].v,l[i+1].p-l[i+1].v))==0)//重合无视，平行无解            {                if(dcmp(Cross(l[i+1].p-l[i].p,l[i+1].p-l[i].v))==0                   &&dcmp(Cross(l[i+1].v-l[i].p,l[i+1].v-l[i].v))==0)continue;                break;            }            p[t++]=GetLineIntersection(l[i].p,l[i].p-l[i].v,l[i+1].p,l[i+1].p-l[i+1].v);        }        if(i<m-1){printf("-1\n");continue;}//平行的时候        if(t==0){printf("-2\n");continue;}//全部重合的时候        int flag=0;        for(i=0;i<t-1;i++)        {            if(p[i]==p[i+1])continue;            flag=1;//相交点不同        }        if(!flag)        {            if(judge(p[0],n))//判断点是否在圆内，在圆就是不符合的点                printf("%.10f %.10f %.10f\n",p[0].x,p[0].y,sqrt(pow(Length(p[0]-cir[0].c),2)-cir[0].r*cir[0].r));            else                printf("-1\n");        }        else printf("-1\n");    }    return 0;}