hdu 4773 圆的反演

第一次接触反演算法。

通过反演圆以求得反演后的直线。

圆的相切即为直线和圆的相切，切点是关键。

值得注意的是，不过中心的直线反演后得到不过中心的圆，圆圆反演后得到另一个圆，因为中途可以得到一条直线，所以可以简化计算。

反演半径不可以随意选取，过大会导致精度问题。

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <string>
#include <iostream>
#include <algorithm>
using namespace std;

#define MEM(a,b) memset(a,b,sizeof(a))
#define REP(i,n) for(int i=1;i<=(n);++i)
#define REV(i,n) for(int i=(n);i>=1;--i)
#define FOR(i,a,b) for(int i=(a);i<=(b);++i)
#define RFOR(i,a,b) for(int i=(a);i>=(b);--i)
#define getmid(l,r) ((l) + ((r) - (l)) / 2)
#define MP(a,b) make_pair(a,b)

typedef long long ll;
typedef pair<int,int> pii;
const int INF = (1 << 30) - 1;
const double eps = 1e-10;

double add(double a,double b) //有误差的浮点加法
{
    if(abs(a + b) < eps * (abs(a) + abs(b))) return 0;
    return a + b;
}

struct Point
{
    double x,y;
    Point(double tx = 0,double ty = 0) : x(tx),y(ty) {}
    Point operator + (Point p)
    {
        return Point(add(x,p.x),add(y,p.y));
    }
    Point operator - (Point p)
    {
        return Point(add(x,-p.x),add(y,-p.y));
    }
    Point operator * (double d)
    {
        return Point(x * d,y * d);
    }
    Point operator / (double d)
    {
        return Point(x / d,y / d);
    }
    Point Move(double a,double d)//从x正方向开始，逆时针
    {
        return Point(x + d * cos(a),y + d * sin(a));
    }
    void Read()
    {
        scanf("%lf%lf",&x,&y);
    }
};

struct Circle
{
    Point o;
    double r;
    Circle(double tx = 0,double ty = 0,double tr = 0) : o(tx,ty),r(tr) {}
    void Read()
    {
        o.Read();
        scanf("%lf",&r);
    }
    void Out()
    {
        printf("%.8f %.8f %.8f\n",o.x,o.y,r);
    }
};

int Sign(double x) //判断x的正负
{
    return (x > eps) - (x < -eps);
}

double Cross(Point a,Point b,Point c) //a为顶点的叉积
{
    return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
}

double Dis(Point a,Point b)//距离
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

Circle c[5];
Point P;
int T,tot;
double R;

Circle Inver(Circle c1) //圆反演圆过程
{
    Circle res;
    double oc1 = Dis(P,c1.o);
    double k1 = 1.0 / (oc1 - c1.r);
    double k2 = 1.0 / (oc1 + c1.r);
    res.r = 0.5 * (k1 - k2) * R * R;
    double oc2 = 0.5 * (k1 + k2) * R * R;
    res.o = P + (c1.o - P) * oc2 / oc1;
    return res;
}

void Mark(Point a,Point b) //记下答案圆，直线反演回圆
{
    ++tot;
    double t = fabs(Cross(a,P,b) / Dis(a,b)); //求出P点到直线的距离
    c[tot].r = R * R / (2.0 * t);
    double d = Dis(a,c[1].o);
    c[tot].o = P + (a - c[1].o) * (c[tot].r / d); //因为向量(a,c[1].o)与公切线垂直，所以可以利用其长度相似出P到c[tot]圆心的距离
}

void Solve()
{
    REP(i,2) c[i] = Inver(c[i]); //将已知两圆反演
    if(c[1].r < c[2].r) swap(c[1],c[2]); //c[1]圆的半径较大
    Point tmp = c[2].o - c[1].o;//基础应该转角多少
    double a1 = atan2(tmp.y,tmp.x); //atan2的经典应用，求出直线的倾斜角！
    double a2 = acos((c[1].r - c[2].r) / Dis(c[1].o,c[2].o));//两圆心应该倾斜多少
    Point P1 = c[1].o.Move(a1 + a2,c[1].r);//计算切点
    Point P2 = c[2].o.Move(a1 + a2,c[2].r);//计算切点
    if(Sign(Cross(P1,c[1].o,P2)) == Sign(Cross(P1,P,P2))) Mark(P1,P2); //保证P与c[1]圆心在公切线同侧
    P1 = c[1].o.Move(a1 - a2,c[1].r);//同样要考虑公切线在下（上）方的情况，计算切点
    P2 = c[2].o.Move(a1 - a2,c[2].r);//计算切点
    if(Sign(Cross(P1,c[1].o,P2)) == Sign(Cross(P1,P,P2))) Mark(P1,P2); //保证P与c[1]圆心在公切线同侧
}

int main()
{
    R = 10.0; //自定义的反演半径，不可过大
    scanf("%d",&T);
    REP(tt,T)
    {
        tot = 2;
        REP(i,2) c[i].Read();
        P.Read();
        Solve();
        printf("%d\n",tot - 2);
        FOR(i,3,tot) c[i].Out();
    }
    return 0;
}