2013杭州region D HDU 4773

Problem of Apollonius

题意：给定两个圆，告诉半径和圆心，它们是相离的，在这两个圆外给定一个点p，求符合条件：过点p的圆且与已知的两个圆外切的所有圆的总数和它们的圆心坐标和半径。

链接:http://acm.hdu.edu.cn/showproblem.php?pid=4773

直接联立方程求解十分复杂，用圆的反演求解此题（几乎是模板题）

圆的反演：由于反演之后相切关系保持不变，长用作相切圆的处理

好懂的反演圆的介绍
反形圆求法
动图理解
$c_1$为原始圆圆心，$c_o$为反演中心， $c_2$为反型圆圆心
其中反形圆的圆心 $c = c_o + \frac{OC_2}{OC_1} (c_1 - c_o)$
又有$\frac{OC_2}{OC_1} = \frac{r_2}{r_1}$

注意：不仅要求的是公切线，为了防止出现内切的情况，我们要求的是“外公切线”，并且要使得P点和反形圆的圆心在外公切线的同一侧。

代码：

#include<bits/stdc++.h>
using namespace std;

int t;
const double eps = 1e-9;
const double pi = acos(-1);
int sign(double x) {
    if (fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;
}
struct point {
    double x, y;
    point(double x = 0, double y = 0) :x(x), y(y) {}
    point operator + (const point &rhs) const {
        return point(x + rhs.x, y + rhs.y);
    }
    point operator - (const point &rhs) const {
        return point(x - rhs.x, y - rhs.y);
    }
    point operator * (const double &rhs) const {
        return point(x * rhs, y * rhs);
    }
    point operator / (const double &rhs) const {
        return point(x / rhs, y / rhs);
    }
    void read() {
        scanf("%lf%lf", &x, &y);
    }
};
double dot(point a, point b) {
    return a.x * b.x + a.y * b.y;
}
double len(point v) {
    return sqrt(dot(v, v));
}
double cross(point a, point b) {
    return a.x * b.y - a.y * b.x;
}
double cross(point A, point B, point C) {
    return (B.x-A.x)*(C.y-A.y) - (B.y-A.y) * (C.x-A.x);
}
struct circle {
    point c;
    double r;
    circle() {}
    circle(point c, double r) :c(c), r(r) {}
    point get_point(double ang) {
        return point(c.x + cos(ang) * r, c.y + sin(ang) * r);
    }
    void read() {
        c.read();
        scanf("%lf", &r);
    }
};

circle get_cir(point p1, point p2, point p3) {
    double x1 = p1.x, x2 = p2.x, x3 = p3.x;
    double y1 = p1.y, y2 = p2.y, y3 = p3.y;
    double a = x1 - x2, b = y1 - y2, c = x1 - x3, d = y1 - y3;
    double e = ((x1 * x1 - x2 * x2) + (y1 * y1 - y2 * y2)) / 2;
    double f = ((x1 * x1 - x3 * x3) + (y1 * y1 - y3 * y3)) / 2;
    double delta = b * c - a * d;
    point pc = point(-(d * e - b * f) / delta, -(a * f - c * e) / delta);
    return circle(pc, len(pc - p1));
}
circle get_inv(circle c) {
    double r = (1 / (len(c.c) - c.r) - 1 / (len(c.c) + c.r)) / 2;
    point p = c.c * (r / c.r);
    return circle(p, r);
}
circle get_inv(point a, point b) {
    point p1 = a * (1.0 / (len(a) * len(a)));
    point p2 = b * (1.0 / (len(b) * len(b)));
    return get_cir(p1, p2, point(0, 0));
}
point p, pa[10], pb[10];
circle c1, c2;
int get_tangents(circle A, circle B, point* a, point* b) {
    int cnt = 0;
    if (A.r < B.r) {
        swap(A, B);
        swap(a, b);
    }
    double d2 = (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
    double rdiff = A.r - B.r;
    double rsum = A.r + B.r;

    double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    double ang = acos((A.r - B.r) / sqrt(d2));
    a[cnt] = A.get_point(base + ang);
    b[cnt] = B.get_point(base + ang);
    if (sign(cross(a[cnt], A.c, b[cnt])) == sign(cross(a[cnt], point(), b[cnt])) &&
        sign(cross(a[cnt], A.c, b[cnt])) == sign(cross(a[cnt], point(), b[cnt])))
    cnt++;
    a[cnt] = A.get_point(base - ang);
    b[cnt] = B.get_point(base - ang);
    if (sign(cross(a[cnt], A.c, b[cnt])) == sign(cross(a[cnt], point(), b[cnt])) &&
        sign(cross(a[cnt], A.c, b[cnt])) == sign(cross(a[cnt], point(), b[cnt])))
    cnt++;
    return cnt;
}

int main() {
    scanf("%d", &t);
    while (t--) {
        c1.read();
        c2.read();
        p.read();
        c1.c = c1.c - p;
        c2.c = c2.c - p;
        circle invc1 = get_inv(c1);
        circle invc2 = get_inv(c2);
        int flag = get_tangents(invc1, invc2, pa, pb);
        printf("%d\n", flag);
        for (int i = 0; i < flag; i++) {
            circle c = get_inv(pa[i], pb[i]);
            c.c = c.c + p;
            printf("%.9f %.9f %.9f\n", c.c.x, c.c.y, c.r);
        }
    }
}