HDU 4773 （圆的反演）

题目：传送门

题目

题意：给定两个圆的半径和圆心，它们是相离的，在这两个圆外给定一个点p，输出所有过点p且与已知的两个

圆外切的圆。

思路

根据圆的反演的一些性质，可以将过点 p 的圆与已知两个圆外切转化为过点 p 的直线，与已知两个圆外切。

圆的反演的一些性质：

1、过反演中心的圆，反形(经过反演之后的图形)为不过反演中心的直线。

2、不过反演中心的直线，反形为过反演中心的圆。

3、不过反演中心的圆，反形也为不过反演中心的圆，并且反演中心为这两个互为反形的圆的位似中心。

4、相切两圆的反象仍相切，若切点恰是反演中心，则其反象为两平行线。

5、相离的两圆反演(反演中心不在圆上)后仍然相离；两圆相切，若反演中心在某圆上，则为反形为相切的直线与圆。

我们以 p 为反演中心，半径自设，然后对已知的两个圆反演，然后求两圆的外切线，将外切线反演为圆即为所求。

圆的反演：参考

#include <bits/stdc++.h>
#define LL long long
#define ULL unsigned long long
#define UI unsigned int
#define mem(i, j) memset(i, j, sizeof(i))
#define rep(i, j, k) for(int i = j; i <= k; i++)
#define dep(i, j, k) for(int i = k; i >= j; i--)
#define pb push_back
#define make make_pair
#define INF 0x3f3f3f3f
#define inf LLONG_MAX
#define PI acos(-1)
#define fir first
#define sec second
#define lb(x) ((x) & (-(x)))
#define dbg(x) cout<<#x<<" = "<<x<<endl;
using namespace std;

const int N = 1e6 + 5;
const double eps = 1e-8;

struct Point {
    double x, y;
    Point(double a = 1.0, double b = 1.0) : x(a), y(b) {}
    Point operator + (const Point &a) { return Point(x+a.x, y+a.y); }
    Point operator - (const Point &a) { return Point(x-a.x, y-a.y); }
    Point operator * (const double a) { return Point(a*x, a*y); }
    Point operator / (const double a) { return Point(a/x, a/y); }
    void Input() { scanf("%lf %lf", &x, &y); }
    void Output() { printf("%.8f %.8f\n", x, y); }
};

struct Circle {
    Point o;
    double r;
    Circle(Point a = Point(), double b = 1.0) : o(a), r(b) { }
    Point getPoint(double alpha) { return o + Point(r*cos(alpha), r*sin(alpha)); }
    void Input() { o.Input(); scanf("%lf", &r); }
    void Output() { printf("%.8f %.8f %.8f\n", o.x, o.y, r); }
};

double dis(Point a, Point b) { return sqrt((a.x-b.x) * (a.x-b.x) + (a.y-b.y) * (a.y-b.y)); }
double Cross(Point a, Point b, Point c) { return (b.x-a.x)*(c.y-a.y) - (c.x-a.x)*(b.y-a.y); }
double Cross(Point a, Point b) { return a.x*b.y - a.y*b.x; }
double Dot(Point a, Point b, Point c) { return (b.x-a.x)*(c.x-a.x) + (b.y-a.y)*(c.y-a.y); }
int dcmp(double x) { return (x > eps) - (x < -eps); }

Point Point_Inver(Circle c0,Point P){
    Point OP = P - c0.o;
    double len = dis(c0.o,P);
    len = len*len;
    return c0.o + OP*( c0.r * c0.r / len );
}
Circle Circle_Inver(Circle c0,Circle a){
    Circle res;
    Point OA = a.o - c0.o;
    double len = dis(a.o,c0.o);
    Point up = c0.o + OA * ( ( len + a.r) / len );
    Point down = c0.o + OA *( (len - a.r) / len );
    up = Point_Inver(c0,up);
    down = Point_Inver(c0,down);
    res.o = (up+down) * 0.5;
    res.r = dis(up,down) * 0.5;
    return res;
}
Circle Line_Inver(Circle c0,Point a,Point b){
    Circle res = Circle();
    double d = fabs(Cross(a,c0.o,b) / dis(a,b));
    res.r = c0.r * c0.r / (2.0 * d);
    double len = Dot(a,b,c0.o) / dis(a,b);
    Point AB = b - a;
    Point c = a + AB * (len/dis(a,b));
    Point CO = c - c0.o;
    res.o = c0.o + CO * (res.r/d);
    //double len = dis(a,c[1].o);
    //res.o = c0.o + (a-c[1].o) * (res.r/len);
    return res;
}

Circle c[20];
Circle c0;
int tot = 0;

void solve() {

    c[0].Input(); c[1].Input();
    c0.o.Input();
    c0.r = 5.00;

    Point p = c0.o;

    c[0] = Circle_Inver(c0, c[0]);
    c[1] = Circle_Inver(c0, c[1]);

    if(c[1].r > c[0].r) swap(c[0], c[1]);

    Point v = c[1].o - c[0].o;

    double alpha = atan2(v.y, v.x);
    double d = dis(c[0].o, c[1].o);
    double beta = acos((c[0].r - c[1].r) / d);
    int k = 2;
    Point a = c[0].getPoint(alpha + beta);
    Point b = c[1].getPoint(alpha + beta);

    if(dcmp(Cross(a, c[0].o, b)) == dcmp(Cross(a, p, b)) &&
       dcmp(Cross(a, c[1].o, b)) == dcmp(Cross(a, p, b)))
        c[++k] = Line_Inver(c0, a, b);

    a = c[0].getPoint(alpha - beta);
    b = c[1].getPoint(alpha - beta);

    if(dcmp(Cross(a, c[0].o, b)) == dcmp(Cross(a, p, b)) &&
       dcmp(Cross(a, c[1].o, b)) == dcmp(Cross(a, p, b)))
        c[++k] = Line_Inver(c0, a, b);

    printf("%d\n", k - 2);

    rep(i, 3, k) c[i].Output();

}


int main() {

    int _; scanf("%d", &_);
    while(_--) solve();

//    solve();

    return 0;
}