题意：给出两个相离的圆O1,O2和圆外一点P，求构造这样的圆：同时与两个圆相外切，且经过点P，输出圆的圆心和半径
分析：画图很容易看出这样的圆要么存在一个，要么存在两个：此题直接解方程是不容易的，先看看反演的定义：已知一圆C，圆心为O，半径为r，如果P与P’在过圆心O的直线上，且，则称P与P'关于O互为反演。
反演的性质：
首先设出反演圆心O和反演半径R
1、圆外一点P与圆内一点P‘会一一对应的反演OP*OP'=R*R
2、经过O的圆，反演后成为不经过O的一条直线
3、不经过O的圆，反演后成为另一个圆，且圆心并不对应
4、不经过O的直线反演后成为一个经过O的圆
5、过O的直线反演后不变
具体参考该http://blog.csdn.net/acdreamers/article/details/16966369
注意：反演中心圆的半径要设大一点，否则会有精度问题
分析：画图很容易看出这样的圆要么存在一个，要么存在两个：此题直接解方程是不容易的，先看看反演的定义：已知一圆C，圆心为O，半径为r，如果P与P’在过圆心O的直线上，且，则称P与P'关于O互为反演。
反演的性质：
首先设出反演圆心O和反演半径R
1、圆外一点P与圆内一点P‘会一一对应的反演OP*OP'=R*R
2、经过O的圆，反演后成为不经过O的一条直线
3、不经过O的圆，反演后成为另一个圆，且圆心并不对应
4、不经过O的直线反演后成为一个经过O的圆
5、过O的直线反演后不变
具体参考该http://blog.csdn.net/acdreamers/article/details/16966369
注意：反演中心圆的半径要设大一点，否则会有精度问题
代码：
#include"stdio.h"
#include"string.h"
#include"stdlib.h"
#include"queue"
#include"algorithm"
#include"string.h"
#include"string"
#include"math.h"
#include"vector"
#include"stack"
#include"map"
#define eps 1e-4
#define inf 0x3f3f3f3f
#define M 609
#define PI acos(-1.0)
using namespace std;
struct node//二维点坐标
{
double x,y;
node(){};
node(double xx,double yy)
{
x=xx;
y=yy;
}
node operator+(node b)
{
return node(x+b.x,y+b.y);
}
node operator-(node b)
{
return node(x-b.x,y-b.y);
}
node operator*(double b)
{
return node(x*b,y*b);
}
double operator*(node b)
{
return x*b.x+y*b.y;
}
double operator^(node b)
{
return x*b.y-y*b.x;
}
};
struct Circle//定义圆心和半径
{
node center;
double r;
}p[M];
struct Line//定义直线一般式的三个参数ABC
{
double A,B,C;
};
struct Line2//定义两条直线
{
Line s,e;
};
double len(node a)//求向量的长度
{
return sqrt(a*a);
}
double dis(node a,node b)//求两个点的距离
{
return len(a-b);
}
double cross(node a,node b,node c)//求叉乘
{
return (b-a)^(c-a);
}
double dot(node a,node b,node c)//求点乘
{
return (b-a)*(c-a);
}
Circle InverCircle(node p,double R,Circle c)//已知反演中心和反演半径，求圆c的反形圆
{
Circle ret;
double pc=dis(p,c.center);
ret.r=R*R/2*(1/(pc-c.r)-1/(pc+c.r));
double pret=R*R/(pc+c.r)+ret.r;
ret.center.x=p.x+(pret/pc)*(c.center.x-p.x);
ret.center.y=p.y+(pret/pc)*(c.center.y-p.y);
return ret;
}
node InterPoint(node p,Line L)//过直线外一点与直线L的垂足
{
node ret;
ret.y=(L.A*L.A*p.y-L.A*L.B*p.x-L.B*L.C)/(L.A*L.A+L.B*L.B);
ret.x=(L.B*L.B*p.x-L.A*L.B*p.y-L.A*L.C)/(L.A*L.A+L.B*L.B);
return ret;
}
Circle InverLine(node p,double R,Line L)//已知不过反演中心的直线，求其反形圆（反形圆过反演中心）
{
Circle ret;
node q=InterPoint(p,L);
double l1=dis(p,q);
double l2=R*R/l1;
ret.r=l2/2;
ret.center.x=p.x+ret.r/l1*(q.x-p.x);//利用相似三角形求解
ret.center.y=p.y+ret.r/l1*(q.y-p.y);
return ret;
}
Line line(node a,node b)//已知两个不同的点，求过这两个点的直线
{
Line l;
l.A=b.y-a.y;
l.B=a.x-b.x;
l.C=b.x*a.y-a.x*b.y;
return l;
}
double disLL(Line L1,Line L2)//两条平行线间垂直距离
{
return fabs(L1.C-L2.C)/sqrt(L1.A*L1.A+L1.B*L1.B);
}
double disPL(node p,Line L)//求点到直线的距离
{
return fabs(L.A*p.x+L.B*p.y+L.C)/sqrt(L.A*L.A+L.B*L.B);
}
node Ratate(node a,double rad)//向量逆时针旋转rad的角度
{
return node(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
Line PointCutCircle(node p,Circle c,int clock)//过圆外一点p且与圆的切线,clock代表两个不同的切线
{
Line ret;
node she;
she=c.center-p;
double pc=dis(c.center,p);
double rad=asin(c.r/pc);
if(clock==-1)//逆时针旋转
{
node she1=Ratate(she,rad);
ret.A=she1.y;
ret.B=-she1.x;
ret.C=she1.x*p.y-p.x*she1.y;
}
if(clock==1)//顺时针旋转
{
node she1=Ratate(she,-rad);
ret.A=she1.y;
ret.B=-she1.x;
ret.C=she1.x*p.y-p.x*she1.y;
}
return ret;
}
Line2 CircleCutCircle(Circle O1,Circle O2)//求两个相离的圆的两条外公切线
{
Line2 L;
Line l1,l2,L1;
if(fabs(O1.r-O2.r)<eps)//当量圆半径相同时
{
L1=line(O1.center,O2.center);
l1.A=L1.A;
l1.B=L1.B;
l1.C=L1.C+fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
l2.A=L1.A;
l2.B=L1.B;
l2.C=L1.C-fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
L.s=l1;
L.e=l2;
}
else//当两个圆的半径不同时
{
if(O1.r>O2.r)
swap(O1,O2);
Circle O;
O.center=O2.center;
O.r=O2.r-O1.r;
L1=PointCutCircle(O1.center,O,-1);
l1.A=L1.A;
l1.B=L1.B;
l1.C=L1.C+fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
l2.A=L1.A;
l2.B=L1.B;
l2.C=L1.C-fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
if(fabs(disPL(O1.center,l1)-O1.r)<eps&&fabs(disPL(O2.center,l1)-O2.r)<eps)
L.s=l1;
if(fabs(disPL(O1.center,l2)-O1.r)<eps&&fabs(disPL(O2.center,l2)-O2.r)<eps)
L.s=l2;
L1=PointCutCircle(O1.center,O,1);
l1.A=L1.A;
l1.B=L1.B;
l1.C=L1.C+fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
l2.A=L1.A;
l2.B=L1.B;
l2.C=L1.C-fabs(O1.r)*sqrt(L1.A*L1.A+L1.B*L1.B);
if(fabs(disPL(O1.center,l1)-O1.r)<eps&&fabs(disPL(O2.center,l1)-O2.r)<eps)
L.e=l1;
if(fabs(disPL(O1.center,l2)-O1.r)<eps&&fabs(disPL(O2.center,l2)-O2.r)<eps)
L.e=l2;
}
return L;
}
int main()
{
int T,cnt;
scanf("%d",&T);
while(T--)
{
Circle C1,C2,C3,C4,C,ans[4];
Line2 L;
node P;
scanf("%lf%lf%lf%lf%lf%lf%lf%lf",&C1.center.x,&C1.center.y,&C1.r,&C2.center.x,&C2.center.y,&C2.r,&P.x,&P.y);
C3=InverCircle(P,30.0,C1);
C4=InverCircle(P,30.0,C2);
L=CircleCutCircle(C3,C4); cnt=0;
C=InverLine(P,30.0,L.s);
if(fabs(dis(C.center,C1.center)-C1.r-C.r)<eps&&fabs(dis(C.center,C2.center)-C2.r-C.r)<eps)
ans[cnt++]=C; C=InverLine(P,30.0,L.e);
if(fabs(dis(C.center,C1.center)-C1.r-C.r)<eps&&fabs(dis(C.center,C2.center)-C2.r-C.r)<eps)
ans[cnt++]=C; printf("%d\n",cnt);
for(int i=0;i<cnt;i++)
printf("%.8lf %.8lf %.8lf\n",ans[i].center.x,ans[i].center.y,ans[i].r);
}
}