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;
}