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 就是反演之后圆可能很大，角度处理细节要注意。

 

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

#define op operator
#define db double
#define cn const
#define cp const P&
#define rt return
using namespace std;
cn db pi = acos(-1.0);
cn db eps = 1e-9;
const int N = 1007;

inline int sig(db x) {rt (x>eps) - (x<-eps);}

int n;
struct P{
    db x, y;
    P(db _x= 0, db _y =0) : x(_x), y(_y) {}
    void in() {scanf("%lf %lf", &x, &y); }
    P op-(cp a)cn {rt P(x-a.x, y-a.y); }
    P op+(cp a)cn {rt P(x+a.x, y+a.y); }
    P op*(db a)cn {rt P(x*a, y*a);}
    db dis() {rt sqrt(x*x + y*y);}
};
P p;

struct CCL{
    P o;
    db r;
    CCL() {}
    CCL(P _o, db _r) : o(_o), r(_r) {}
    void in() {o.in(), scanf("%lf", &r); }
}c[N];

CCL Inverse(CCL ci) {
    CCL T;
    db t = (ci.o - p).dis();
    db x = 97. / (t - ci.r);
    db y = 97. / (t + ci.r);
    T.r = (x - y) / 2.0;
    db s = (x + y) / 2.0;
    T.o = p + (ci.o - p) * (s / t);
    rt T;
}

struct L{
    db alpha;
    int t;
    L(db a=0, int t = 0) : alpha(a), t(t) {}
    bool op<(cn L& a)cn {
        if(sig(alpha - a.alpha)) rt sig(alpha - a.alpha) < 0;
        rt t > a.t;
    }
};
L s[N<<3];
int ct;

void add_ang(db a1, db a2) {
    s[ct++] = L(a1, 1), s[ct++] = L(a2, -1);
    if(sig(a2 - pi) > 0) s[ct++] = L(a1-pi-pi, 1), s[ct++] = L(a2-pi-pi, -1);
    if(sig(a1 + pi) < 0) s[ct++] = L(a1+pi+pi, 1), s[ct++] = L(a2+pi+pi, -1);
}

void qie(CCL A, CCL B) {
    db d = (A.o-B.o).dis();
    bool f = 0;
    db rdiff = A.r - B.r, rsum = A.r + B.r;
    if(sig(d + A.r - B.r) <= 0) {
        s[ct++] = L(-pi-pi, 1), s[ct++] = L(pi+pi, -1); return ;
    }
    if(sig(d + B.r - A.r) < 0) return ;

    db base = atan2(B.o.y - A.o.y, B.o.x - A.o.x);
    db alpha = acos(rdiff / d);

    if(sig(d - rsum) <= 0) {
        add_ang(base - alpha, base + alpha);
    } else {
        db ang = acos(rsum / d);
        add_ang(base - alpha, base - ang);
        add_ang(base + ang, base + alpha);
    }
}

void solve() {
    for(int i = 0; i < n; ++i) c[i] = Inverse(c[i]);

    int mx = 0;
    bool flag = 0;
    for(int i = 0; i < n; ++i) {
        ct = 0, flag = 0;
        for(int j = 0; j < n; ++j) {
            if(i == j) continue;
            qie(c[i], c[j]);
        }
        sort(s, s+ct);

        int mv = 0;
        for(int j = 0; j < ct; ++j) {
            mv += s[j].t;
            mx = max(mx, mv);
        }
    }
    printf("%d
", mx+1);
}

int main()
{
    while(scanf("%d", &n) != EOF) {
        p.in();
        for(int i = 0; i < n; ++i) c[i].in();
        solve();
    }
    return 0;
}