HDU4773 Problem of Apollonius（反演变换）

HDU4773 Problem of Apollonius（反演变换）

题意：

给定两个相离的圆，一个点P，求出过点P并且和两个圆都外切的所有圆。

思路:

以P为反演中心将两个圆O1,O2分别反演，得到O1',O2'，则和O1,O2相切的圆反演后是一条不经过P点且与O1'、O2‘相切的直线。显然这条直线就是O1'、O2’的公切线。因为要求和O1、O2外切的圆，因此反演后的直线必定使O1'、O2'、P在圆的同一侧。筛选出符合条件的公切线再反演回去就是所求的圆了。

代码：

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const double PI=acos(-1.0);
int sgn(double x) {
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
struct Point {
    double x,y;
    Point() {}
    Point(double _x,double _y) {
        x = _x;
        y = _y;
    }
    Point operator -(const Point &b)const {
        return Point(x - b.x,y - b.y);
    }
    Point operator +(const Point &b)const {
        return Point(x + b.x,y +b.y);
    }
    double operator ^(const Point &b)const {
        return x*b.y - y*b.x;
    }
    double operator *(const Point &b)const {
        return x*b.x + y*b.y;
    }
    Point operator *(double b)const {
        return Point(x*b,y*b);
    }
    Point Rotate(double rad){ //逆时针旋转
        return Point(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));
    }
    double len(){//Vector
        return sqrt(x*x+y*y);
    }
    void stdd(){//Vector
        double le=len();
        x/=le;
        y/=le;
    }
    void input(){
        scanf("%lf%lf",&x,&y);
    }
    void output(){
        printf("(%.8f,%.8f) ",x,y);
    }
};
struct Circle {
    Point o;
    double r;
    Circle(){}
    void input(){
        o.input();
        scanf("%lf",&r);
    }
    void output(){
        printf("%.8f %.8f %.8f
",o.x,o.y,r);
    }
};
double xmult(Point p0,Point p1,Point p2) { //p0p1 X p0p2
    return (p1-p0)^(p2-p0);
}
double dist(Point a,Point b) {
    return sqrt( (b - a)*(b - a) );
}
double getArea(Point a,Point b,Point c){//三角形面积S
    return fabs(xmult(a,b,c))/2.0;
}
double distLP(Point a,Point b,Point p){//p到ab的距离
    return getArea(a,b,p)*2/dist(a,b);
}
Circle CtC(Circle C,Point p,double r){//圆到圆的反演
    Circle res;
    double t = dist(C.o,p);
    double x = r*r / (t - C.r);
    double y = r*r / (t + C.r);
    res.r = (x - y) / 2.0;
    double s = (x + y) / 2.0;
    res.o = p + (C.o - p) * (s / t);
    return res;
}
Circle LtC(Point a,Point b,Point p,double r){//直线到圆的反演
    double d=distLP(a,b,p);
    d=r*r/d;
    Circle c;
    c.r=d/2;
    Point v1;
    if(xmult(a,b,p)>0)
        v1=(a-b).Rotate(PI/2);
    else
        v1=(b-a).Rotate(PI/2);
    v1.stdd();
    c.o=p+v1*c.r;
    return c;
}
Circle c[10];
Point p;
bool check(Point a,Point b,Point p){
    if(sgn(xmult(a,b,p))==sgn(xmult(a,b,c[0].o))&&sgn(xmult(a,b,p))==sgn(xmult(a,b,c[1].o)))
        return 1;
    else return 0;
}
vector<Circle>ans;
void solve(){
    ans.clear();
    double r=1;
    c[0]=CtC(c[0],p,r);
    c[1]=CtC(c[1],p,r);
    if(c[0].r>c[1].r)swap(c[0],c[1]);
    Point v1=c[0].o-c[1].o;
    double d=dist(c[0].o,c[1].o);
    Point a,b;
    double alpha=acos((c[1].r-c[0].r)/d);
    v1.stdd();
    a=c[1].o+v1.Rotate(-alpha)*c[1].r;
    b=c[0].o+v1.Rotate(-alpha)*c[0].r;
    if(check(a,b,p))ans.push_back(LtC(a,b,p,r));
    a=c[1].o+v1.Rotate(alpha)*c[1].r;
    b=c[0].o+v1.Rotate(alpha)*c[0].r;
    if(check(a,b,p))ans.push_back(LtC(a,b,p,r));
}
int main () {
    int T;
    scanf("%d",&T);
    while(T--){
        c[0].input();c[1].input();
        scanf("%lf%lf",&p.x,&p.y);
        solve();
        printf("%d
",ans.size());
        for(int i=0;i<ans.size();i++){
            ans[i].output();
        }
    }
}