Commentator problem（CF 2）

题目链接
题目大意：
给定三个圆，询问是否存在点满足该点与三个圆夹角均相等，若存在多组解返回夹角最大值。

圆外一点到两圆夹角均相等：
即 sina = sinb = r1 / d1 = r2 / d2
即 d1 / d2 = r1 / r2
引出阿波罗尼斯圆定义：
平面上相异两点A B， 且P 与A B在同一平面内， PA / PB = k (k > 0 且 k != 1), 则 P点的轨迹是圆，把他称为阿波罗尼斯圆。
注意到当k等于1时P点的运动轨迹在A，B的中垂线上，那么此题就转化为直线求交、直线与圆求交、圆与圆求交的简单问题了。
给出阿氏圆的推导结论：

设A（x1, y1）, B (x2, y2)  , k = a.r / b.r
那么有 圆心坐标为
{（k * k * x2 - x1）/ (k * k - 1), (k * k * y2 - y1) / (k * k -1)} 
半径向量为
{（x2 - x1）, (y2 - y1)} * k / (k * k - 1)

输出前记得判断点是否合法（不在任一圆内即可）

#include <bits/stdc++.h>

using namespace std;
// `计算几何模板`
const double eps = 1e-4;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 1010;
//`Compares a double to zero`
int sgn(double x){
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    else return 1;
}
//square of a double
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);
    }
    bool operator == (Point b)const{
        return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
    }
    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);
    }
    //`计算pa  和  pb 的夹角`
    //`就是求这个点看a,b 所成的夹角`
    //`测试 LightOJ1203`
    double rad(Point a,Point b){
        Point p = *this;
        return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));
    }
    //`化为长度为r的向量`
    Point trunc(double r){
        double l = len();
        if(!sgn(l))return *this;
        r /= l;
        return Point(x*r,y*r);
    }
    //`逆时针旋转90度`
    Point rotleft(){
        return Point(-y,x);
    }
    //`顺时针旋转90度`
    Point rotright(){
        return Point(y,-x);
    }
    //`绕着p点逆时针旋转angle`
    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;
    }
    //`根据一个点和倾斜角angle确定直线,0<=angle<pi`
    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)));
        }
    }
    //求线段长度
    double length(){
        return s.distance(e);
    }
    //`返回直线倾斜角 0<=angle<pi`
    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;
    }
    //`点和直线关系`
    //`1  在左侧`
    //`2  在右侧`
    //`3  在直线上`
    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;
    }
    //`直线和线段相交判断`
    //`-*this line   -v seg`
    //`2 规范相交`
    //`1 非规范相交`
    //`0 不相交`
    int linecrossseg(Line v){
        int d1 = sgn((e-s)^(v.s-s));
        int d2 = sgn((e-s)^(v.e-s));
        if((d1^d2)==-2) return 2;
        return (d1==0||d2==0);
    }
    //`两直线关系`
    //`0 平行`
    //`1 重合`
    //`2 相交`
    int linecrossline(Line v){
        if((*this).parallel(v))
            return v.relation(s)==3;
        return 2;
    }
    //`求两直线的交点`
    //`要保证两直线不平行或重合`
    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(s == e) return p.distance(s);
        if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)
            return min(p.distance(s),p.distance(e));
        return dispointtoline(p);
    }
    //`返回线段到线段的距离`
    //`前提是两线段不相交，相交距离就是0了`
    double dissegtoseg(Line v){
        return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));
    }
    //`返回点p在直线上的投影`
    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;
    }
    //输入
    void input(){
        p.input();
        scanf("%lf",&r);
    }
    //`点和圆的关系`
    //`0 圆外`
    //`1 圆上`
    //`2 圆内`
    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;
    }
    //`两圆的关系`
    //`5 相离`
    //`4 外切`
    //`3 相交`
    //`2 内切`
    //`1 内含`
    //`需要Point的distance`
    //`测试：UVA12304`
    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;
        return 0;// int 需要返回值（除条件外） 
    }
    //`求两个圆的交点，返回0表示没有交点，返回1是一个交点，2是两个交点`
    //`需要relationcircle`
    //`测试：UVA12304`
    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 pointcrossline(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;
    }
};

circle a, b, c;

Line get_l(circle a, circle b){
    Line l = {b.p ,a.p};
    double ang = l.angle();
    if(sgn(ang-pi/2)>=0) ang -= pi/2;
    else ang += pi/2;
    Point mid = (a.p+b.p)/2;
    return Line{mid, ang};
}

circle get_c(circle a, circle b){
    double k = a.r/b.r;
    circle c;
    c.p = (Point){(k*k*b.p.x-a.p.x)/(k*k-1), (k*k*b.p.y-a.p.y)/(k*k-1)};
    Point cr = {b.p.x-a.p.x, b.p.y-a.p.y};
    cr= cr*k/(k*k-1);
    c.r = cr.len();
    return c;
}

bool judge(Point k){
    if(a.relation(k) || b.relation(k) || c.relation(k)) return false;
    return true;
}

void output(Point k){
    if(judge(k)) printf("%.5lf %.5lf\n", k.x, k.y);
}

void wk(circle a, circle b, circle c){
    auto L = get_l(a, b);
    auto A = get_c(a, c), B = get_c(b, c);
    Point a1, a2;
    int t = A.pointcrosscircle(B, a1, a2);
    if(t == 0) return ;
    else if(t == 1){
        if(L.relation(a1) == 3) output(a1);
    }
    else{
        if(L.relation(a1) == 3 && L.relation(a2) == 3){
            if(judge(a1) && judge(a2)){
                if(sgn(a1.distance(a.p)-a2.distance(a.p))<0) output(a1);
                else output(a2);
            }
            else output(a1), output(a2);
        }
        else if(L.relation(a1) == 3){
            output(a1);
        }
        else if(L.relation(a2) == 3){
            output(a2);
        }
    }
}

int main(){
    a.input(), b.input(), c.input();
    if(sgn(a.r-b.r)==0 && sgn(a.r-c.r)==0){
        auto A = get_l(a, b), B = get_l(a, c), C = get_l(b, c);
        if(A.linecrossline(B) == 2){
            Point p = A.crosspoint(B);
            output(p);
        }
    }
    else if(sgn(a.r-b.r)!=0 && sgn(a.r-c.r)!=0 && sgn(b.r-c.r)!=0){
        auto A = get_c(a, b), B = get_c(a, c), C = get_c(b, c);
        Point a1, a2, c1, c2;
        int t = A.pointcrosscircle(B, a1, a2);
        if(t == 0) return 0;
        else if(t == 1){
            int t1 = A.pointcrosscircle(C, c1, c2);
            if(t1){
                if(a1 == c1 || a1 == c2) output(a1);
            }
        }
        else{
            int t1 = A.pointcrosscircle(C, c1, c2);
            if(t1 == 0) return 0;
            else if(t1 == 1){
                if(c1 == a1 || c1 == a2) output(c1);
            }
            else{
                if((a1 == c1 && a2 == c2) || (a1 == c2 && a2 == c1)){
                    if(judge(a1) && judge(a2)){
                        if(sgn(a1.distance(a.p)-a2.distance(a.p))<0) output(a1);
                        else output(a2);
                    }
                    else output(a1), output(a2);
                }
                else if(a1 == c1 || a1 == c2){
                    output(a1);
                }
                else if(a2 == c1 || a2 == c2){
                    output(a2);
                }
            }
        }
    }
    else{
        if(sgn(a.r-b.r)==0){
            wk(a, b, c);
        }
        else if(sgn(a.r-c.r)==0){
            wk(a, c, b);
        }
        else{
            wk(b, c, a);
        }
    }
    return 0;
}