【HDU】4773 Problem of Apollonius

题意

给定相离的两个圆（圆心坐标以及半径）以及圆外的一个定点$P$，求出过点$P$的且与已知的两个圆外切的所有圆（输出总数+圆心、半径）。

分析

如果强行解方程，反正我是不会。
本题用到新姿势：圆的反演。
二维上的圆的反演通常是指定一个圆$C$为基础，其圆心$O$为反演中心，其半径$r$为反演半径。对于平面上任意一个非反演中心的点$P$，都有且仅有一个反演点$P'$，满足$OP \cdot OP' = r^2$且$P'$在$OP$射线上。对于任意一个二维上的图形，其反形就是图形上的所有点的反演点组成的。
于是可以证明：

1. $C$内的点的反演点在$C$外，$C$上的点的反演点就是自身，$C$外的点的反演点在$C$内。
2. 任意一条不过反演中心的直线，其反形为过反演中心的圆。（至于过反演中心，我只能理解为这是极限意义下有点已经无限逼近了反演中心所以就算经过）
3. 任意一个不过反演中心的圆，其反形为不过反演中心的圆，且反演中心为两圆的位似中心。

于是对于本题，交$3$个点的圆，可以看做一条直线关于$P$点反演得到的圆。由于反演点唯一，所以这条直线肯定与给出的两个圆的反演圆各相交一个点。所以就是这两个圆的反演圆的公切线！我们发现，如果是内公切线，反演成圆后会把一个圆包含，因此不符合题意。所以我们只需要考虑外公切线即可。
但是外公切线的反形圆也可能会出现把给出的两个圆包含的情况，画一下图就能发现这种情况只出现在$p$和另外两个反演圆不在公切线的同一侧。

题解

本题要注意反演半径不要开太小，否则会有精度问题。（我是设成10）

圆$(C_1, r_1)$关于圆$(P, r)$反演得到反形$(C_2, r_2)$，我们来求$r_2$。
根据反演式子可以得到：

$$\begin{align} (PC_1+r_1) \cdot (PC_2-r_2) = & r^2 \\ (PC_1-r_1) \cdot (PC_2+r_2) = & r^2 \\ \end{align}$$

消去$PC_2$可以解得： $$r_2 = \frac{r^2}{2} \left( \frac{1}{PC_1-r_1} - \frac{1}{PC_1+r_1} \right)$$ 那么再根据$C_1$的反演点是$C_2$，我们也能解出$C_2$来。

然后我们需要求切线的反演圆$(O, r_3)$了，假设切线与圆$(C_2, r_2)$交于$A$。
首先$r_3$可以由$P$到切线的距离$t$通过反演定义式求出:

$$2r_3 \cdot t = r^2 \Rightarrow r_3 = \frac{r^2}{2t}$$ 然后由于直线$OP$与$AC_1$是平行的（都与切线垂直）
所以根据相似来求出圆心坐标： $$O = P + \frac{r_3}{r_2} \overrightarrow{C_2 A}$$

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
typedef double lf;
const lf eps=1e-8;
inline int dcmp(lf x) {
	return x>-eps?x>=eps:-1;
}
struct ip {
	lf x, y;
	ip(lf _x=0, lf _y=0) : x(_x), y(_y) {}
	void scan() {
		scanf("%lf%lf", &x, &y);
	}
};
typedef ip iv;
ip operator + (ip a, iv b) {
	return ip(a.x+b.x, a.y+b.y);
}
iv operator - (ip a, ip b) {
	return iv(a.x-b.x, a.y-b.y);
}
lf operator ^ (iv a, iv b) {
	return a.x*b.y-a.y*b.x;
}
iv operator * (lf k, iv a) {
	return iv(a.x*k, a.y*k);
}
inline lf sqr(lf x) {
	return x*x;
}
lf dis(ip a, ip b) {
	return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
int onleft(ip a, ip b, iv v) {
	return dcmp(v^(a-b))==1;
}
struct ic {
	ip o;
	lf r;
	void scan() {
		o.scan();
		scanf("%lf", &r);
	}
	ip gen(lf a) {
		return ip(o.x+r*cos(a), o.y+r*sin(a));
	}
	void print() {
		printf("%.8f %.8f %.8f\n", o.x, o.y, r);
	}
}o[2], ans[2], p;
ic inv(ic a) {
	ic b;
	lf d=dis(a.o, p.o), r2=sqr(p.r),
	   k1=r2/(d-a.r), k2=r2/(d+a.r);
	b.r=0.5*(k1-k2);
	b.o=p.o+0.5*(k1+k2)/d*(a.o-p.o);
	return b;
}
ic inv(ip a, ip b) {
	ic c;
	c.r=sqr(p.r)/(abs(((b-a)^(p.o-a))/dis(a, b))*2);
	c.o=p.o+(c.r/o[0].r)*(a-o[0].o);
	return c;
}
int cal() {
	int tot=0;
	o[0]=inv(o[0]);
	o[1]=inv(o[1]);
	if(o[0].r<o[1].r) {
		swap(o[0], o[1]);
	}
	iv t=o[1].o-o[0].o;
	lf k1=atan2(t.y, t.x),
	   k2=acos((o[0].r-o[1].r)/dis(o[0].o, o[1].o));
	ip a=o[0].gen(k1+k2), b=o[1].gen(k1+k2);
	t=b-a;
	if(onleft(o[0].o, a, t)==onleft(p.o, a, t)) {
		ans[tot++]=inv(a, b);
	}
	a=o[0].gen(k1-k2), b=o[1].gen(k1-k2);
	t=b-a;
	if(onleft(o[0].o, a, t)==onleft(p.o, a, t)) {
		ans[tot++]=inv(a, b);
	}
	return tot;
}
int main() {
	int T;
	scanf("%d", &T);
	while(T--) {
		o[0].scan();
		o[1].scan();
		p.o.scan();
		p.r=10;
		int tot=cal();
		printf("%d\n", tot);
		for(int i=0; i<tot; ++i) {
			ans[i].print();
		}
	}
	return 0;
}