hdu4773 Problem of Apollonius【反演变换】

题目大意：

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

解题思路：

这个blog解释的挺不错:

http://blog.csdn.net/acdreamers/article/details/16966369。

不过本蒟蒻还是理一下自己的思路。

首先了解一下反演的定义：

设反演中心为O，反演半径为R，那么若有两点P,P′连线过O点，且OP∗OP′=R2，则称P和P′关于O互为反演

这个定义可推广到任意图形上。

然后就是反演的性质了。

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