HDU - 4773 Problem of Apollonius 圆的反演

题目链接点这里

圆的反演主要有3条性质

//1.不过反演中心的圆经过反演变换仍然是一个不过反演中心的圆.

//2.不过反演中心的直线经过反演变换是一个经过反演中心的圆.

//3.反演变换不改变图形的相切性.

然后这道题就解决了。。

//注意精度，直接用rad计算圆弧中点，
#include<iostream>
#include<cstdio>
#include<math.h>
#include<algorithm>
#include<map>
#include<set>
#include<bitset>
#include<stack>
#include<queue>
#include<string.h>
#include<string>
#include<cstring>
#include<vector>
#include<time.h>
#include<stdlib.h>
using namespace std;
#define INF 0x3f3f3f3f
#define INFLL 0x3f3f3f3f3f3f3f3f
#define FIN freopen("input.txt","r",stdin);
#define mem(x,y) memset(x,y,sizeof(x));
typedef unsigned long long ULL;
typedef long long LL;
#define fuck(x) cout<<x<<endl;
const int  MX=333;
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
typedef pair<pair<int,int>,int> PIII;
typedef pair<int,int> PII;
const double eps=1e-8;
const double PI=acos(-1);
struct Point {
    double x, y;
    Point() {}
    Point(double x,double y):x(x),y(y) {}
};
typedef Point Vector;
int dcmp(double x) { //返回x的正负
    if(fabs(x)<eps)return 0;
    return x<0?-1:1;
}
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);
}
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;//如果改成整形记得加LL
}
double Cross(Vector A,Vector B) { //叉积
    return A.x*B.y-A.y*B.x;//如果改成整形记得加LL
}
//向量长度
double Length(Vector A) {
    return sqrt(Dot(A,A));
}
//2个向量之间的夹角
double Angle(Vector A,Vector B) {
    return acos(Dot(A,B)/Length(A)/Length(B));
}
//向量的极角
double angle(Vector v) {
    return atan2(v.y,v.x);
}
//将A向量逆时针旋转rad
Vector Rotate( Vector A,double rad) {
    return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
//返回A的逆时针旋转90度的单位法向量
Vector Normal(Vector A) {
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}
//计算2条直线P+tv和Q+tw的交点，请先确保不是平行(v！=w)
Point GetLineIntersection(Point P,Vector v,Point Q,Vector w) {
    Vector u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}
//P到直线AB的距离
double DistanceToLine(Point P,Point A,Point B) {
    Vector v1=B-A,v2=P-A;
    return fabs(Cross(v1,v2))/Length(v1);
}
//园
struct Circle {
    Point c;
    double r;
    Circle() {}
    Circle(Point c,double r):c(c),r(r) {}
    Point getpoint(double rad) {
        return Point(c.x+cos(rad)*r,c.y+sin(rad)*r);
    }
};
//有向线段
struct Line {
    Point p;
    Vector v;//方向向量，左边为半平面
    double ang;//极角
    Line() {}
    Line(Point p,Vector v):p(p),v(v) {
        ang=atan2(v.y,v.x);
    }
    bool operator<(const Line &L)const { //极角排序
        return ang<L.ang;
    }
    Point point(double t) {
        return p + v*t;
    }
    Line move(double d) { //向左边平移d单位
        return Line(p + Normal(v)*d, v);
    }
};
//返回切线条数-1表示无穷多条切线
//a[i],b[i]分别是第i条切线在圆A，圆B上的切点
int getCircleTangents(Circle A,Circle B,Point *a,Point *b) {
    int cnt=0;//切线条数
    if(A.r<B.r) {
        swap(A,B);
        swap(a,b);
    }
    double d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);//圆心距
    double rdiff=A.r-B.r;
    double rsum=A.r+B.r;
    if(dcmp(d2-rdiff*rdiff)<0) return 0;//内含
    double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);//求出圆心连线的极角
    if(dcmp(d2)==0&&dcmp(A.r-B.r)==0) return -1;//两圆重合
    if(dcmp(d2-rdiff*rdiff)==0) {  //内切 一条外公切线

        a[cnt]=A.getpoint(base);
        b[cnt]=B.getpoint(base);
        cnt++;
        return 1;
    }
    //有两条外公切线
    double ang=acos((A.r-B.r)/sqrt(d2));//求出切线与圆心连线的夹角
    a[cnt]=A.getpoint(base+ang);
    b[cnt]=B.getpoint(base+ang);
    cnt++;
    a[cnt]=A.getpoint(base-ang);
    b[cnt]=B.getpoint(base-ang);
    cnt++;
    if(d2==rsum*rsum) { //两圆外切
        a[cnt]=A.getpoint(base);
        b[cnt]=B.getpoint(base+PI);
        cnt++;
    } else if(d2>rsum*rsum) {
        double ang=acos((A.r+B.r)/sqrt(d2));//求出内公切线和圆心连线的夹角
        a[cnt]=A.getpoint(base+ang);
        b[cnt]=B.getpoint(PI+base+ang);
        cnt++;
        a[cnt]=A.getpoint(base-ang);
        b[cnt]=B.getpoint(PI+base-ang);
        cnt++;
    }
    return cnt;
}
//求圆u关于c的反演圆
Circle InvCircletoCircle(Circle C,Circle u) {
    Circle T;
    double t = Length(C.c-u.c);
    double x = 1.0 / (t - u.r);
    double y = 1.0 / (t + u.r);
    T.r = (x - y)*C.r*C.r/ 2.0;
    double s = (x + y)*C.r*C.r/ 2.0;
    T.c = C.c + (u.c - C.c) * (s / t);
    return T;
}
//求直线u关于c的反演圆
Circle InvLinetoCircle(Circle C,Line u) {
    Circle T;
    Point w=GetLineIntersection(C.c,Normal(u.v),u.p,u.v);//垂足
    double dis = Length(w-C.c);
    double t = C.r*C.r/ dis;
    Point p=C.c+(w-C.c)/dis*t;
    T.r=t/2;
    T.c=(p+C.c)/2;
    return T;
}
//求圆u关于c的反演直线
//先保证u经过C的圆心
Line InvCircletoLine(Circle C,Circle u) {
    Line T;
    double t=C.r*C.r/(2*u.r);
    Point p=(u.c-C.c)/Length(u.c-C.c)*t+C.c;//垂足
    T.p=p;
    T.v=Normal(u.c-C.c);
    return T;
}
int main() {
    //FIN;
    int T;
    cin>>T;
    while(T--) {
        Circle a,b,inva,invb,C;
        scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&a.c.x,&a.c.y,&a.r,&b.c.x,&b.c.y,&b.r,&C.c.x,&C.c.y);
        C.r=3;
        inva=InvCircletoCircle(C,a);
        invb=InvCircletoCircle(C,b);
        //cout<<inva.c.x<<" "<<inva.c.y<<" "<<inva.r<<endl;
        //cout<<invb.c.x<<" "<<invb.c.y<<" "<<invb.r<<endl;
        Point pa[10],pb[10];
        int t=getCircleTangents(inva,invb,pa,pb);
        vector<Circle> ans;
        //cout<<t<<endl;
        for(int i=0; i<t; i++) {
            int t1=dcmp(Cross(pa[i]-pb[i],inva.c-pb[i]));
            int t2=dcmp(Cross(pa[i]-pb[i],invb.c-pb[i]));
            int t3=dcmp(Cross(pa[i]-pb[i],C.c-pb[i]));
            //cout<<t1<<t2<<t3<<endl;
            if(t1==t2&&t2==t3) {
                ans.push_back(InvLinetoCircle(C,Line(pa[i],pa[i]-pb[i])));
            }
        }
        cout<<ans.size()<<endl;
        for(auto i:ans) {
            printf("%.8f %.8f %.8f\n",i.c.x,i.c.y,i.r);
        }


    }
    return 0;
}