hdu 4773 Problem of Apollonius

题解：圆的反演，测试一下模版。

他人详细讲解

#include"bits/stdc++.h"

using namespace std;
const double eps = 1e-12;
struct Point {
    double x, y;
    Point() {}
    Point(double x,double y):x(x),y(y) {}
    void print(){
        printf("%.8f %.8f\n",x,y);
    }
};
typedef Point Vector;
int dcmp(double x) { //返回x的正负
    if(fabs(x)<eps) return 0;
    return x<0?-1:1;
}
Vector operator-(Vector A,Vector B) {return Vector(A.x - B.x, A.y - B.y);}
Vector operator+(Vector A,Vector B) {return Vector(A.x + B.x, A.y + B.y);}
Vector operator*(Vector A,double p) {return Vector(A.x*p, A.y*p);}
Vector operator/(Vector A,double p) {return Vector(A.x/p, A.y/p);}
bool operator<(const Point&a,const Point&b) {return a.x<b.x||(a.x==b.x&&a.y<b.y);}
bool operator==(const Point&a,const Point&b) {return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;}
double Dot(Vector A,Vector B) { //点积
    return A.x*B.x+A.y*B.y;//如果改成整形记得加LL
}
double Cross(Vector A,Vector B) { //叉积
    return A.x*B.y-A.y*B.x;//如果改成整形记得加LL
}
//向量长度
double Length(Vector A) {
    return sqrt(Dot(A,A));
}
//返回A的逆时针旋转90度的单位法向量
Vector Normal(Vector A) {
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}

struct Circle {
    Point c;
    double r;
    Circle() {}
    Circle(Point c,double r):c(c),r(r) {}
    Point getpoint(double rad) {
        return Point(c.x+cos(rad)*r,c.y+sin(rad)*r);
    }
    void print(){
        printf("%.8f %.8f %.8f\n",c.x, c.y, r);
    }
};
struct Line {
    Point p;
    Vector v;//方向向量，左边为半平面
    double ang;//极角
    Line() {}
    Line(Point p,Vector v):p(p),v(v) {
        ang=atan2(v.y,v.x);
    }
    bool operator<(const Line &L)const { //极角排序
        return ang<L.ang;
    }
    Point point(double t) {
        return p + v*t;
    }
    Line move(double d) { //向左边平移d单位
        return Line(p + Normal(v)*d, v);
    }
    void print(){
        p.print();
        v.print();
    }
};
Point GetLineIntersection(const Line& a, const Line& b)
{
    Vector u = a.p-b.p;
    double t = Cross(b.v, u) / Cross(a.v, b.v);
    return a.p+a.v*t;
}

Circle c[2];
vector<Circle> ans;
//求圆u关于c的反演圆
Circle InvCircletoCircle(Circle C,Circle u) {
    Circle T;
    double t = Length(C.c-u.c);
    double x = 1.0 / (t - u.r);
    double y = 1.0 / (t + u.r);
    T.r = (x - y)*C.r*C.r/ 2.0;
    double s = (x + y)*C.r*C.r/ 2.0;
    T.c = C.c + (u.c - C.c) * (s / t);
    return T;
}
//求直线u关于c的反演圆
//先确保直线u不经过反演圆心，不然反演过后就是它本身
Circle InvLinetoCircle(Circle C,Line u) {
    Circle T;
    Point w=GetLineIntersection(Line(C.c,Normal(u.v)),Line(u.p, u.v));//垂足
    double dis = Length(w-C.c);
    double t = C.r*C.r/ dis;
    Point p=C.c+(w-C.c)/dis*t;
    T.r=t/2;
    T.c=(p+C.c)/2;
    return T;
}
//求圆u关于c的反演直线
//先保证u经过C的圆心
Line InvCircletoLine(Circle C,Circle u) {
    Line T;
    double t=C.r*C.r/(2*u.r);
    Point p=(u.c-C.c)/Length(u.c-C.c)*t+C.c;//垂足
    T.p=p;
    T.v=Normal(u.c-C.c);
    return T;
}

void solve(Point p)
{
    Circle C = Circle(p,10.0);
    for(int i = 0; i < 2; i++) c[i] = InvCircletoCircle(C,c[i]);
    if(c[0].r < c[1].r) swap(c[0],c[1]);
    Point temp = c[1].c - c[0].c;
    double base = atan2(temp.y,temp.x);
    double ang = acos( (c[0].r - c[1].r) / Length(c[0].c - c[1].c));

    Point p1 = c[0].getpoint(base+ang);
    Point p2 = c[1].getpoint(base+ang);

    if(dcmp(Cross(p1-c[0].c,p2-c[0].c)) == dcmp(Cross(p1-p,p2-p)) )
        ans.push_back(InvLinetoCircle(C,Line(p1,p2-p1)));
    p1 = c[0].getpoint(base-ang);
    p2 = c[1].getpoint(base-ang);
    if(dcmp(Cross(p1-c[0].c,p2-c[0].c)) == dcmp(Cross(p1-p,p2-p)) )
        ans.push_back(InvLinetoCircle(C, InvCircletoLine(C, InvLinetoCircle(C,Line(p1,p2-p1)) ) ) );
}

int tot = 0;
int main()
{
#ifdef LOCAL
    freopen("in.txt","r",stdin);
#endif // LOCAL
    int T;
    cin>>T;
    while(T--){
        Point p;
        scanf("%lf %lf %lf %lf %lf %lf %lf %lf",&c[0].c.x, &c[0].c.y, &c[0].r, &c[1].c.x, &c[1].c.y, &c[1].r, &p.x, &p.y);
        tot = 0;
        ans.clear();
        solve(p);
        printf("%d\n",ans.size());
        for(int i = 0; i < ans.size(); i++) ans[i].print();
    }
    return 0;
}