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