bzoj1337 最小圆覆盖

Description

---

给出平面上N个点,N<=10^5.请求出一个半径最小的圆覆盖住所有的点

Solution

---

新东西当然是要资瓷的，计算几何都不太记得了，何况是这么神奇的算法
显然一个最优的圆一定过至少两点，那么我们默认1作为初始圆心。
每次加入第i个点的时候如果不能覆盖全部说明第i个点一定在圆上，返回去枚举另一个点做一个新的圆重新判断
如果由两个点做出的圆仍然不满足全覆盖那么找第三个点，使得三点确定一个圆，这个圆就是前i个点的答案了

为了保证复杂度需要随机乱排n个点。由于i个点中加入一个点使得在圆外的几率大致是3/i的，因此时间复杂度期望是O（n）的

我大概是弱智选手了，推半天初中内容求直线方程交点，gg

Code

---

#include <cstdio>
#include <cmath>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
const int N=200005;
struct pos {double x,y;}p[N],cir;
double r=0;
double sqr(double x) {return x*x;}
double get_dis(pos a,pos b) {return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}
pos get_mid(pos a,pos b) {return (pos){(a.x+b.x)*0.5,(a.y+b.y)*0.5};}
pos get_center(pos a,pos b,pos c) {
    pos cen1=get_mid(a,b);
    pos cen2=get_mid(a,c);
    pos ret;
    double k1,k2,b1,b2;
    if (a.y==b.y) {
        k2=-1.0/((c.y-a.y)/(c.x-a.x));
        b2=cen2.y-k2*cen2.x;
        ret=(pos){cen1.x,cen1.x*k2+b2};
    } else if (a.y==c.y) {
        k1=-1.0/((b.y-a.y)/(b.x-a.x));
        b1=cen1.y-k1*cen1.x;
        ret=(pos){cen2.x,cen2.x*k1+b1};
    } else {
        k1=-1.0/((b.y-a.y)/(b.x-a.x));
        b1=cen1.y-k1*cen1.x;
        k2=-1.0/((c.y-a.y)/(c.x-a.x));
        b2=cen2.y-k2*cen2.x;
        ret.x=(b1-b2)/(k2-k1); ret.y=k1*ret.x+b1;
    }
    return ret;
}
int main(void) {
    int n; scanf("%d",&n);
    rep(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y);
    std:: random_shuffle(p+1,p+n+1);
    cir=p[1];
    rep(i,2,n) {
        double d=get_dis(p[i],cir);
        if (d<=r) continue;
        r=0; cir=p[i];
        rep(j,1,i-1) {
            if (get_dis(p[j],cir)<=r)  continue;
            r=get_dis(p[j],p[i])*0.5;
            cir=get_mid(p[j],p[i]);
            rep(k,1,j-1) {
                if (get_dis(p[k],cir)<=r) continue;
                cir=get_center(p[i],p[j],p[k]);
                r=get_dis(cir,p[i]);
            }
        }
    }
    printf("%.3lf\n", r);
    return 0;
}