HDU 4773 Problem of Apollonius——圆反演

题面

　　HDU4773

解析

　　  大概是圆反演的模板吧。

　　以点$P(x3, y3)$为反演中心，任意长为反演半径，将两个已知圆反演，设反演后的圆为$A'$, $B'$，所求圆反演后为一条直线，根据题目中的要求，该直线为两圆的外公切线。因此我们只需要求出两圆的外公切线即可。

　　然后会发现WA了，因为题目中还有一个要求，所求圆要外切于两圆，即反演变换后反演中心$P$和$A'$的圆心要在同侧。

　　还有一个我一开始做错了的地方，原来的圆心$O$反演后就不是新的圆心了！！！可以连接$PO$，求其与圆的两个交点，两个交点的反演点中点才是新的圆心。

　代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const double eps = 1e-11, Pi = 3.14159265358979323;

int T, cnt, num;
double R;

struct node{
    double fir, sec;
    node(){}
    node(double x, double y) {
        fir = x;
        sec = y;
    }
}p, c[5], d[5];

node operator + (node x, node y)
{
    return node(x.fir + y.fir, x.sec + y.sec);
}

node operator - (node x, node y)
{
    return node(x.fir - y.fir, x.sec - y.sec);
}

node operator * (node x, double y)
{
    return node(x.fir * y, x.sec * y);
}

node operator / (node x, double y)
{
    return node(x.fir / y, x.sec / y);
}

//Rotate vector counterclockwise by alpha
node rotate(node x, double y)//rad
{
    return node(x.fir * cos(y) - x.sec * sin(y), x.fir * sin(y) + x.sec * cos(y));
}

//Cross product
double crp(node x, node y)
{
    return x.fir * y.sec - x.sec * y.fir;
}

//Dot product
double dtp(node x, node y)
{
    return x.fir * y.fir + x.sec * y.sec;
}

//Distance
double dist(node x)
{
    return sqrt(dtp(x, x));
}

//Find the intersection of two straight lines
//x, y are the vector starting points
//u, v are the direction vectors
node itpt(node x, node u, node y, node v)
{
    node t = x - y;
    double rat = crp(v, t) / crp(u, v);
    return x + u * rat;
}

struct circ{
    node o;
    double r;
}a, b, aa, bb, ans[5];

//Find the common tangent of two circles
void tagt(circ x, circ y)
{
    if(x.r < y.r)    swap(x, y);
    double len = dist(x.o - y.o);

    //External common tangent
    double alp = acos((x.r - y.r) / len);
    node tmp = y.o - x.o, t = rotate(tmp, alp);
    c[++cnt] = x.o + t / len * x.r;
    tmp = x.o - y.o; t = rotate(tmp, Pi + alp);
    d[cnt] = y.o + t / len * y.r;
    tmp = y.o - x.o; t = rotate(tmp, 2 * Pi - alp);
    c[++cnt] = x.o + t / len * x.r;
    tmp = x.o - y.o; t = rotate(tmp, Pi - alp);
    d[cnt] = y.o + t / len * y.r;

    //Internal common tangent
    /*alp = acos((x.r + y.r) / len);
    tmp = y.o - x.o; t = rotate(tmp, alp);
    c[++cnt] = x.o + t / len * x.r;
    tmp = x.o - y.o; t = rotate(tmp, alp);
    d[cnt] = y.o + t / len * y.r;
    tmp = y.o - x.o; t = rotate(tmp, 2 * Pi - alp);
    c[++cnt] = x.o + t / len * x.r;
    tmp = x.o - y.o; t = rotate(tmp, 2 * Pi - alp);
    d[cnt] = y.o + t / len * y.r;*/
}

int check(double x)
{
    return (x > eps) - (x < -eps);
}

void out(double x)
{
    printf("%.8f ", fabs(x) < 0.000000005? 0: x);
}

int main()
{
    scanf("%d", &T);
    R = 1.0;
    double len1, len2;
    node tmp, t;
    while(T --)
    {
        scanf("%lf%lf%lf%lf%lf%lf%lf%lf", &a.o.fir, &a.o.sec, &a.r, &b.o.fir, &b.o.sec, &b.r, &p.fir, &p.sec);
        len1 = dist(a.o - p);
        tmp = p + (a.o - p) / len1 * (len1 - a.r);
        len2 = dist(tmp - p);
        tmp = p + (tmp - p) / len2 * (R / len2);
        t = p + (a.o - p) / len1 * (len1 + a.r);
        len2 = dist(t - p);
        t = p + (t - p) / len2 * (R / len2);
        aa.o = (tmp + t) / 2;
        aa.r = dist(aa.o - tmp);

        len1 = dist(b.o - p);
        tmp = p + (b.o - p) / len1 * (len1 - b.r);
        len2 = dist(tmp - p);
        tmp = p + (tmp - p) / len2 * (R / len2);
        t = p + (b.o - p) / len1 * (len1 + b.r);
        len2 = dist(t - p);
        t = p + (t - p) / len2 * (R / len2);
        bb.o = (tmp + t) / 2;
        bb.r = dist(bb.o - tmp);


        cnt = num = 0;
        tagt(aa, bb);    
        for(int i = 1; i <= cnt; ++i)
        {
            if(fabs(crp(c[i] - p, d[i] - p)) < eps || check(crp(d[i] - c[i], p - c[i])) !=  check(crp(d[i] - c[i], aa.o - c[i])))    continue;
            ++ num;
            tmp = c[i] - d[i];
            tmp = rotate(tmp, Pi / 2);
            tmp = itpt(p, tmp, d[i], c[i] - d[i]);
            len1 = dist(tmp - p);
            tmp = p + (tmp - p) / len1 * (R / len1);
            ans[num].o = (p + tmp) / 2;
            ans[num].r = dist(ans[num].o - p);
        }    
        printf("%d
", num);
        for(int i = 1; i <= num; ++i)
        {
            out(ans[i].o.fir);
            out(ans[i].o.sec);
            printf("%.8f
", ans[i].r);
        }
    }
    return 0;
}