HDU 4773 Problem of Apollonius 圆的反演

传送门
坐标系上给定两个不相交的圆，求所有过一定点且与这两个圆都外切的圆。

反演变换

如果把以那个定点为原点进行反演变换，那么过定点的圆就变成了直线，不过定点的圆还是圆，问题就转化成了求那两个圆的公切线。求出所有公切线再反演回圆再判是否都是外切。

#include <bits/stdc++.h>
using namespace std;
#define prt(k) cerr<<#k" = "<<k<<endl
typedef long long LL;
const int inf = 0x3f3f3f3f;
const double PI = acos(-1.0);
const double eps = 1e-8;
int cmp(double x)
{
    return x < -eps ? -1 : x > eps;
}
struct P
{
    double x, y, r;
    P() { x = 0, y = 0, r = 1; }
    P(double xx, double yy, double rr =  1)
    {
        x = xx, y = yy, r = rr;
    }
    P getP(double a)
    {
        return P(x + cos(a)*r, y + sin(a)*r);
    }
    void in() { scanf("%lf%lf%lf", &x, &y, &r); }
    P operator * (double k) { return P(x*k, y*k, r*k); }
    P operator / (double k) {
        assert(abs(k) > 1e-8);
        return P(x/k, y/k);
    }
    void out()
    {
        printf("%.10f %.10f %.10f\n", x, y, r);
    }
};
double pf(double x) { return x * x; }
P operator-(P a, P  b)
{
    return P(a.x - b.x, a.y - b.y);
}
P operator+(P a, P  b)
{
    return P(a.x + b.x, a.y + b.y);
}
double abs(P p) { return sqrt(pf(p.x) + pf(p.y)); }
double dis(P a, P b) { return abs(a - b); }
typedef P Circle ;
P a, b, o;
double dot(P a, P b) { return a.x*b.x + a.y*b.y; }
P cuizhu(P p, P s, P t) ///点到直线垂足
{
    double r = dot(t-s, p-s)/dot(t-s, t-s);
    return s + (t-s)*r;
}
const double r = 1;
/// 圆的反演（结果是一个圆），反演中心O为原点，反演满足OC * OC' = r*r  这里 r = 1
/// 圆不经过反演中心
Circle inv(Circle a)
{
    double d1  =  abs(a);
    double r1 = a.r;
    double r2 = r1 / (d1 * d1 - r1 * r1) * r * r;
    double d2 = (r1 * r2 + r * r) / d1;
    double x = a.x / d1 * d2;
    double y = a.y / d1 * d2;
    P p = a * (r2 / r1);
    p.r = r2;
    return P(x, y, r2);
}
Circle inv(P s, P t) /// 直线的反演，直线不经过反演中心
{
    P O(0,0);
    P p = cuizhu(O, s, t);
    double d = dis(O, p);
    p = p * (0.5 / d / d);
    p.r = dis(p, O);
    return p;
}
/// 圆的切线，返回切线条数，见白书 267 页
int getTan(Circle A, Circle B, P a[], P b[])
{
    int cnt = 0;
    if (A.r < B.r) { swap(A, B); swap(a, b); }
    double d2 = pf(dis(A, B));
    double rdiff = (A.r - B.r);
    double rsum = A.r + B.r;
    if (d2 < pf(rdiff)) return 0;
    double base = atan2(B.y - A.y, B.x - A.x);
    if (d2 == 0 && A.r == B.r)  return -1;
    if (d2 == pf(rdiff)) {
        a[cnt] = A.getP(base); b[cnt++] = B.getP(base);
        return cnt;
    }
    double ang = acos((A.r - B.r) / sqrt(d2));
    a[cnt] = A.getP(base+ang); b[cnt++] = B.getP(base + ang);
    a[cnt] = A.getP(base - ang); b[cnt++] = B.getP(base - ang);
    if (d2 == rsum * rsum) {
        a[cnt] = A.getP(base); b[cnt++] = B.getP(PI + base);
    }
    else if (d2 > rsum * rsum) {
        double ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.getP(base + ang); b[cnt++] = B.getP(base + ang);
        a[cnt] = A.getP(base - ang); b[cnt++] = B.getP(base - ang);
    }
    return cnt;
}
bool XiangQie(P a, P b)
{
    return cmp(dis(a, b) - a.r - b.r) == 0;
}
bool in(Circle C, P p)
{
    return cmp(dis(C, p) - C.r) == 0;
}
Circle A, B;
bool check(Circle ans)
{
    return XiangQie(A, ans) && XiangQie(B, ans) && in(ans, o);
}
int main()
{
    int re; scanf("%d", &re);
    while (re--) {
        a.in(); b.in();
        A = a, B = b;
        scanf("%lf%lf", &o.x, &o.y);
        a.x -= o.x; a.y -= o.y;
        b.x -= o.x; b.y -= o.y;
        a = inv(a); b = inv(b);
        P p1[55], p2[55];
        int cnt = getTan(a, b, p1, p2);
        vector<Circle> v;
        for (int i=0;i<cnt;i++) {
            Circle ans = inv(p1[i], p2[i]);
            ans.x += o.x; ans.y += o.y;
            if (check(ans)) v.push_back(ans);
        }
        printf("%d\n", v.size());
        for (auto C : v) C.out();
    }
}