[圆的反演] HDU 4773 Problem of Apollonius

反演详见ACdreamer的blog

反演的基本性质

• 不过反演中心的一条直线反演成一个过反演中心圆 反之亦然
• 不过反演中心的一个圆反演成另一个不过反演中心的圆

那题目就很显然了先把两个圆反演 变成两个圆 在求公切线 再反演回去

可能用到的几何方法

• 两圆的公切线
• 两个相交圆的交点

最后上几张图

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long double ld;

const ld eps=1e-9;
inline int sgn(ld a){
  if (fabs(a)<eps) return 0; return a<0?-1:1;
}
inline ld sqr(ld x) { return  x*x; }
struct PP{
  ld x,y;
  PP(ld x=0,ld y=0):x(x),y(y){ }
  friend PP operator +(PP a,PP b) { return PP(a.x+b.x,a.y+b.y); }
  friend PP operator -(PP a,PP b) { return PP(a.x-b.x,a.y-b.y); }
  friend PP operator *(PP a,ld b) { return PP(a.x*b,a.y*b); }
  friend PP operator /(PP a,ld b) { return PP(a.x/b,a.y/b); }
  PP rot(){ return PP(-y,x); }
  friend ld cross(PP a,PP b) { return a.x*b.y-a.y*b.x; }
  friend ld operator *(PP a,PP b) { return a.x*b.x+a.y*b.y; }
  friend ld dist(PP a,PP b){ return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y)); }
  void read(){ double _x,_y; scanf("%lf%lf",&_x,&_y); x=_x; y=_y; }
};

struct CC{
  PP c; ld r;
  CC() { }
  CC(PP _c,ld _r) { c=_c; r=_r; }
  CC trans(PP o,ld r0){
    ld oc=dist(o,c),r2=(1/(oc-r)-1/(oc+r))*r0*r0/2;
    ld oc2=r0*r0/(oc-r)-r2;
    PP c2=(c-o)/oc*oc2+o;
    return CC(c2,r2);
  }
  friend pair<PP,PP> cross(CC A,CC B){
    ld L=dist(A.c,B.c);
    ld x0=A.c.x+(B.c.x-A.c.x)*(A.r*A.r-B.r*B.r+L*L)/2/L/L;
    ld y0=A.c.y+(B.c.y-A.c.y)*(A.r*A.r-B.r*B.r+L*L)/2/L/L;
    PP E=PP(x0,y0);
    ld L1=dist(E,A.c),L2=sqrt(A.r*A.r-L1*L1);
    PP l1=E-A.c,l2=(l1/L1*L2).rot(),l3=l2.rot().rot();
    return make_pair(E+l2,E+l3);
  }
  friend bool jud(CC &A,CC &B){
    return !sgn(dist(A.c,B.c)-A.r-B.r);
  }
  void read(){ double _r; c.read(); scanf("%lf",&_r); r=_r; }
  void print() { printf("%lf %lf %lf\n",(double)c.x,(double)c.y,(double)r); }
};

const ld R=4;
PP P; CC C1,C2,C3,C4;
CC c1,c2; PP p,l,p1,p2;

int cnt;
CC ans[10];

inline void calc(){
  if (!sgn(cross(p1-P,p2-P))) return;
  PP t1=P-p1; ld tmp=t1*(p2-p1)/dist(p2,p1);
  PP t=(p2-p1)/dist(p2,p1)*tmp;
  PP T=p1+t;
  ld r=R*R/dist(P,T)/2;
  PP t2=(T-P)/dist(T,P)*r;
  PP c=P+t2;
  ans[++cnt]=CC(c,r);
  if (!(jud(ans[cnt],C1) && jud(ans[cnt],C2)))
    cnt--;
}

int main(){
  int Q;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  scanf("%d",&Q);
  while (Q--){
    cnt=0;
    C1.read(); C2.read(); P.read();
    C3=C1.trans(P,R),C4=C2.trans(P,R);
    if (!sgn(C3.r-C4.r)){
      l=C3.c-C4.c;
      l=l.rot()/dist(C3.c,C4.c)*C3.r;
      p1=C3.c+l; p2=C4.c+l;
      calc();
      p1=C3.c-l; p2=C4.c-l;
      calc();
      printf("%d\n",cnt);
      for (int i=1;i<=cnt;i++) ans[i].print();
      continue;
    }
    if (C3.r<C4.r) swap(C3,C4);
    c1=CC((C3.c+C4.c)/2,dist(C3.c,C4.c)/2);
    c2=CC(C3.c,C3.r-C4.r);
    pair<PP,PP> cro=cross(c1,c2);

    p=cro.first;
    l=(p-C3.c)/(C3.r-C4.r)*C4.r;
    p1=p+l; p2=C4.c+l;
    calc();

    p=cro.second;
    l=(p-C3.c)/(C3.r-C4.r)*C4.r;
    p1=p+l; p2=C4.c+l;
    calc();

    printf("%d\n",cnt);
    for (int i=1;i<=cnt;i++) ans[i].print();
  }
  return 0;
}