HOJ 13102 Super Shuttle （圆的反演变换）

HOJ 13102 Super Shuttle

链接：http://49.123.82.55/online/?action=problem&type=show&id=13102

题意：给定一个点 p 和 n 个圆，做某个经过点 p的圆穿过尽可能多的圆，问可穿过最多的圆是多少。

思路：圆的反演变换：

　　给出反演极O和反演幂k>0，作点A的反演点A′。

　　　　令k=r^2，作出反演基圆⊙O（r），

　　　　1）若点A在⊙O（r）外，则过点A作圆的切线(两条),两个切点相连与OA连线交点就是点A′。

　　　　2）若点A在⊙O（r）内，则把上述过程逆过来：连结OA，过点A作直线垂直于OA，直线与⊙O（r）的交点处的切线的交点就是点A′。

　　　　3）若点A在⊙O（r）上，反演点A′就是点A自身。

推荐看 ACdreamers 的博客。

　　我们取反演点就是点P。如果你看懂了反演，就可以知道，我们所做的圆经过反演变成一条直线（经过点P），而其它圆还是圆，问题就变成了一条直线穿过尽可能多的圆。当中的tricks 就是反演之后圆可能很大，角度处理细节要注意。

代码：

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <cmath>
  5 #include <algorithm>
  6 
  7 #define op operator
  8 #define db double
  9 #define cn const
 10 #define cp const P&
 11 #define rt return
 12 using namespace std;
 13 cn db pi = acos(-1.0);
 14 cn db eps = 1e-9;
 15 const int N = 1007;
 16 
 17 inline int sig(db x) {rt (x>eps) - (x<-eps);}
 18 
 19 int n;
 20 struct P{
 21     db x, y;
 22     P(db _x= 0, db _y =0) : x(_x), y(_y) {}
 23     void in() {scanf("%lf %lf", &x, &y); }
 24     P op-(cp a)cn {rt P(x-a.x, y-a.y); }
 25     P op+(cp a)cn {rt P(x+a.x, y+a.y); }
 26     P op*(db a)cn {rt P(x*a, y*a);}
 27     db dis() {rt sqrt(x*x + y*y);}
 28 };
 29 P p;
 30 
 31 struct CCL{
 32     P o;
 33     db r;
 34     CCL() {}
 35     CCL(P _o, db _r) : o(_o), r(_r) {}
 36     void in() {o.in(), scanf("%lf", &r); }
 37 }c[N];
 38 
 39 CCL Inverse(CCL ci) {
 40     CCL T;
 41     db t = (ci.o - p).dis();
 42     db x = 97. / (t - ci.r);
 43     db y = 97. / (t + ci.r);
 44     T.r = (x - y) / 2.0;
 45     db s = (x + y) / 2.0;
 46     T.o = p + (ci.o - p) * (s / t);
 47     rt T;
 48 }
 49 
 50 struct L{
 51     db alpha;
 52     int t;
 53     L(db a=0, int t = 0) : alpha(a), t(t) {}
 54     bool op<(cn L& a)cn {
 55         if(sig(alpha - a.alpha)) rt sig(alpha - a.alpha) < 0;
 56         rt t > a.t;
 57     }
 58 };
 59 L s[N<<3];
 60 int ct;
 61 
 62 void add_ang(db a1, db a2) {
 63     s[ct++] = L(a1, 1), s[ct++] = L(a2, -1);
 64     if(sig(a2 - pi) > 0) s[ct++] = L(a1-pi-pi, 1), s[ct++] = L(a2-pi-pi, -1);
 65     if(sig(a1 + pi) < 0) s[ct++] = L(a1+pi+pi, 1), s[ct++] = L(a2+pi+pi, -1);
 66 }
 67 
 68 void qie(CCL A, CCL B) {
 69     db d = (A.o-B.o).dis();
 70     bool f = 0;
 71     db rdiff = A.r - B.r, rsum = A.r + B.r;
 72     if(sig(d + A.r - B.r) <= 0) {
 73         s[ct++] = L(-pi-pi, 1), s[ct++] = L(pi+pi, -1); return ;
 74     }
 75     if(sig(d + B.r - A.r) < 0) return ;
 76 
 77     db base = atan2(B.o.y - A.o.y, B.o.x - A.o.x);
 78     db alpha = acos(rdiff / d);
 79 
 80     if(sig(d - rsum) <= 0) {
 81         add_ang(base - alpha, base + alpha);
 82     } else {
 83         db ang = acos(rsum / d);
 84         add_ang(base - alpha, base - ang);
 85         add_ang(base + ang, base + alpha);
 86     }
 87 }
 88 
 89 void solve() {
 90     for(int i = 0; i < n; ++i) c[i] = Inverse(c[i]);
 91 
 92     int mx = 0;
 93     bool flag = 0;
 94     for(int i = 0; i < n; ++i) {
 95         ct = 0, flag = 0;
 96         for(int j = 0; j < n; ++j) {
 97             if(i == j) continue;
 98             qie(c[i], c[j]);
 99         }
100         sort(s, s+ct);
101 
102         int mv = 0;
103         for(int j = 0; j < ct; ++j) {
104             mv += s[j].t;
105             mx = max(mx, mv);
106         }
107     }
108     printf("%d
", mx+1);
109 }
110 
111 int main()
112 {
113     while(scanf("%d", &n) != EOF) {
114         p.in();
115         for(int i = 0; i < n; ++i) c[i].in();
116         solve();
117     }
118     return 0;
119 }