【HDU4773】Problem of Apollonius（圆的反演）

传送门

---

题解：

题意：给两个相离的圆和圆外一点 $P$，求所有过 $P$ 且与两圆外切的圆。

如果直接考虑列方程，只能列出两个方程，而且是平方和开根的形式，不好解。

由于要求的圆要过一个定点，考虑反演。

反演的定义请自行参见百度百科。

反演是平面上点到点的一个映射，除反演中心外所有点形成一一对应，且该映射的逆映射就是其本身。

反演有几个性质

1. 经过反演中心的直线反形就是其自身
2. 不经过反演中心的直线反形是过反演中心的圆
3. 经过反演中心的圆反形是不过反演中心的直线。
4. 不经过反演中心的圆，反形也是一个不经过反演中心的圆，且两个圆关于反演中心位似。
5. 除反演中心外，平面上的所有点有唯一的反演点。

根据性质 5 反演前相切的圆的反形仍然相切，不管它的反形是直线还是圆。

（为了防止杠精，直线是半径为无穷大的圆，不明白这个概念的请放下你手里的键盘，冷静）

所以我们要求一个过某个点的满足某个条件的圆，直接以这个点为反演中心做反演变换，计算满足条件的直线，然后再反演回去就行了。

对于这道题，要求与原来的两个圆外切。画一下图或者感性理解一下就是要求两个圆心和反演中心三个点在所求直线的同一侧。

求反形就各凭本事了。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define db double
#define cs const

using std::cerr;
using std::cout;

cs db eps=1e-10;

inline int sign(db x){return x<-eps?-1:(x>eps?1:0);}

struct Pnt{
	db x,y;Pnt(){}Pnt(cs db &_x,cs db &_y):x(_x),y(_y){}
	friend Pnt operator+(cs Pnt &a,cs Pnt &b){return Pnt(a.x+b.x,a.y+b.y);}
	friend Pnt operator-(cs Pnt &a,cs Pnt &b){return Pnt(a.x-b.x,a.y-b.y);}
	friend Pnt operator*(cs Pnt &a,cs db &b){return Pnt(a.x*b,a.y*b);}
	friend Pnt operator/(cs Pnt &a,cs db &b){return Pnt(a.x/b,a.y/b);}
	friend db operator*(cs Pnt &a,cs Pnt &b){return a.x*b.y-b.x*a.y;}
	inline db len()cs{return sqrt(x*x+y*y);}
};

struct Cir{
	Pnt o;db r;Cir(){}
	inline Pnt rot(db a)cs{
		return Pnt(o.x+r*cos(a),o.y+r*sin(a));
	}
};

Cir P,C1,C2;

inline db area(cs Pnt &a,cs Pnt &b,cs Pnt &c){
	return (b-a)*(c-a);
}
inline db dis(cs Pnt &a,cs Pnt &b){
	return (a-b).len();
}
inline db dis(cs Pnt &a,cs Pnt &b,cs Pnt &c){
	return fabs(area(a,b,c))/dis(a,b);
}
inline Pnt inter(db a1,db b1,db c1,db a2,db b2,db c2){
	if(a1==0){db y=-c1/b1;return Pnt((-c2-b2*y)/a2,y);}
	db y=(a2*c1/a1-c2)/(b2-b1*a2/a1);return Pnt(-(c1+b1*y)/a1,y);
}

inline void inv(Cir &C){
	db d=dis(C.o,P.o),s=P.r*P.r/(d-C.r),t=P.r*P.r/(d+C.r);
	C.r=(s-t)/2;C.o=P.o+(C.o-P.o)/d*(t+C.r);
}
inline Cir inv(cs Pnt &a,cs Pnt &b){
	db a1=b.y-a.y,b1=a.x-b.x,c1=a.y*b.x-a.x*b.y;
	Cir C;C.o=inter(a1,b1,c1,b1,-a1,a1*P.o.y-b1*P.o.x);
	C.r=P.r*P.r/dis(a,b,P.o)/2;db d=dis(C.o,P.o);
	C.o=P.o+(C.o-P.o)*C.r/d;return C;
}

inline bool sameside(cs Pnt &a,cs Pnt &b,cs Pnt &st,cs Pnt &ed){
	return sign(area(st,ed,a))==sign(area(st,ed,b));
}

int ct;
Pnt st[3],ed[3];

void calc_tan(cs Cir &a,cs Cir &b){
	ct=0;db dlt=acos((a.r-b.r)/dis(a.o,b.o));
	db rad=atan2(b.o.y-a.o.y,b.o.x-a.o.x);
	st[ct]=a.rot(rad-dlt),ed[ct]=b.rot(rad-dlt);
	if(sameside(P.o,a.o,st[ct],ed[ct]))++ct;
	st[ct]=a.rot(rad+dlt),ed[ct]=b.rot(rad+dlt);
	if(sameside(P.o,a.o,st[ct],ed[ct]))++ct;
}

void Main(){
	int T;scanf("%d",&T);
	while(T--){
		scanf("%lf%lf%lf",&C1.o.x,&C1.o.y,&C1.r);
		scanf("%lf%lf%lf",&C2.o.x,&C2.o.y,&C2.r);
		scanf("%lf%lf",&P.o.x,&P.o.y);P.r=10;
		inv(C1),inv(C2);calc_tan(C1,C2);
		printf("%d\n",ct);
		for(int re i=0;i<ct;++i){
			Cir t=inv(st[i],ed[i]);
			printf("%.6f %.6f %.6f\n",t.o.x,t.o.y,t.r);
		}
	}
}

void file(){
#ifdef zxyoi
	freopen("apollonius.in","r",stdin);
#endif
}
signed main(){file();Main();return 0;}