hdu4773 Problem of Apollonius【反演变换】

题目大意：

给两个相离的圆和一个不在圆上的点O，问经过该点且与两个圆都外切的圆有哪些。

解题思路：

这个blog解释的挺不错:
http://blog.csdn.net/acdreamers/article/details/16966369。
不过本蒟蒻还是理一下自己的思路。

首先了解一下反演的定义：
设反演中心为$O$，反演半径为$R$，那么若有两点$P,P'$连线过$O$点，且$OP*OP'=R^2$，则称$P$和$P'$关于 O <script type="math/tex" id="MathJax-Element-1648">O</script>互为反演
这个定义可推广到任意图形上。

然后就是反演的性质了。
1、不经过O的直线反演后成为一个经过O的圆
2、经过O的圆，反演后成为不经过O的一条直线
2、不经过O的圆，反演后成为另一个圆，且两个圆关于O位似
4、过O的直线反演后不变
5、反演不改变相切性

那再来看看本题。

由于直接求圆难度很大，但求两个圆的公切线较简单（其实也挺麻烦，全靠画图yy），所以我们可以先求两个圆关于O点的反演（反演半径随便），再求反演后两个圆的公切线，再将公切线关于O反演回来，根据性质1和5，就得到与原来两个圆相切的圆了。

注意题目要求的是都外切，所以所求切线需在反演后两圆圆心及O点同侧（画画图就知道啦，其他情况的切线分别对应都内切及一内一外的情况）。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
#include<ctime>
#include<vector>
#include<set>
using namespace std;

int getint()
{
    int i=0,f=1;char c;
    for(c=getchar();(c!='-')&&(c<'0'||c>'9');c=getchar());
    if(c=='-')f=-1,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}


const double R=1,eps=1e-11;
struct point
{
    double x,y;
    point(){}
    point(double _x,double _y):
        x(_x),y(_y){}
    inline friend point operator + (const point &a,const point &b)
    {return point(a.x+b.x,a.y+b.y);}
    inline friend point operator - (const point &a,const point &b)
    {return point(a.x-b.x,a.y-b.y);}
    inline friend point operator * (const point &a,double b)
    {return point(a.x*b,a.y*b);}
    inline friend double operator * (const point &a,const point &b)
    {return a.x*b.x+a.y*b.y;}
    inline friend double cross(const point &a,const point &b)
    {return a.x*b.y-a.y*b.x;} 
    inline  point n(){return point(-y,x);}
    inline double dis(){return sqrt(x*x+y*y);}
}O;
struct circle
{
    point O;double r;
}C1,C2,C3,C4,C5,C6;
struct line
{
    point st,ed;
}l1,l2;
int bz1,bz2;

double Pow(double x)
{
    return x*x;
}

void getline()
{
    if(C3.r==C4.r)
    {
        point t=(C4.O-C3.O).n();
        t=t*(1/t.dis());
        l1.st=C3.O+t*C3.r,l1.ed=C4.O+t*C4.r;
        l2.st=C3.O-t*C3.r,l2.ed=C4.O-t*C4.r;
    }
    else
    {
        l1.st=l2.st=C3.O+(C3.O-C4.O)*(C3.r/(C4.r-C3.r));
        point t=C3.O-l1.st;
        double Sin=C3.r/t.dis(),d=C3.r*Sin;
        t=t*(1/t.dis());
        l1.ed=l2.ed=C3.O-t*d;
        t=t.n(),d=sqrt(Pow(C3.r)-Pow(d));
        l1.ed=l1.ed+t*d,l2.ed=l2.ed-t*d;
    }
}

point getcross(line a,line b)
{
    double x=cross(a.ed-a.st,b.st-a.st);
    double y=cross(b.ed-a.st,a.ed-a.st);
    point t=b.st+(b.ed-b.st)*(x/(x+y));
    return t;
}

void solve(line a,circle &C,int &bz)
{
    if(abs(cross(a.ed-O,a.st-O))<eps)return;
    line b;
    b.st=O;b.ed=O+((a.ed-a.st).n());
    point p=getcross(a,b);
    if((O-p)*(C3.O-p)<0)return;
    bz=1;
    point t=p-O;t=t*(1/t.dis());
    C.r=Pow(R)/(p-O).dis()*0.5;
    C.O=O+t*C.r;
}

int main()
{
    //freopen("lx.in","r",stdin);
    //freopen("lx.out","w",stdout);
    int T=getint();
    while(T--)
    {
        bz1=bz2=0;
        C1.O.x=getint(),C1.O.y=getint(),C1.r=getint();
        C2.O.x=getint(),C2.O.y=getint(),C2.r=getint();
        O.x=getint(),O.y=getint();
        C3.r=Pow(R)*C1.r/(Pow((O-C1.O).dis())-Pow(C1.r));
        C3.O=O+(C1.O-O)*((Pow(R)/((C1.O-O).dis()+C1.r)+C3.r)/(C1.O-O).dis());
        C4.r=Pow(R)*C2.r/(Pow((O-C2.O).dis())-Pow(C2.r));
        C4.O=O+(C2.O-O)*((Pow(R)/((C2.O-O).dis()+C2.r)+C4.r)/(C2.O-O).dis());
        if(C3.r>C4.r)swap(C3,C4);
        getline();
        solve(l1,C5,bz1);
        solve(l2,C6,bz2);
        printf("%d\n",bz1+bz2);
        if(bz1)printf("%.8f %.8f %.8f\n",C5.O.x,C5.O.y,C5.r);
        if(bz2)printf("%.8f %.8f %.8f\n",C6.O.x,C6.O.y,C6.r);
    }
    return 0;
}