bzoj1337 最小圆覆盖

Description

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

Input

第一行给出数字N，现在N行，每行两个实数x,y表示其坐标.

Output

输出最小半径，输出保留三位小数.

Sample Input

4

1 0

0 1

0 -1

-1 0

Sample Output

1.000

分析：
经典算法之最小圆覆盖

#include<bits/stdc++.h>

using namespace std;

const double eps=1e-8;
const double Pi=acos(-1.0);
struct node{
    double x,y;
    node (double xx=0,double yy=0) {
        x=xx; y=yy;
    }
};
node p[100005];

node operator +(const node &a,const node &b) {return node(a.x+b.x,a.y+b.y);}
node operator -(const node &a,const node &b) {return node(a.x-b.x,a.y-b.y);}
node operator *(const node &a,const double &b) {return node(a.x*b,a.y*b);}
node operator /(const node &a,const double &b) {return node(a.x/b,a.y/b);}

double Dot(node a,node b) {return a.x*b.x+a.y*b.y;}
double Cross(node a,node b) {return a.x*b.y-a.y*b.x;}
double lenth(node a) {return sqrt(Dot(a,a));}

int n;
double r;
node c;

int dcmp(double x) {
    if (fabs(x)<eps) return 0;
    else if (x>0) return 1;
    else return -1;
}

node rotate(node a,double t) {
    return node(a.x*cos(t)-a.y*sin(t),a.x*sin(t)+a.y*cos(t));
}

node jiao(node P,node v,node Q,node w) {
    node u=P-Q;
    double t=Cross(w,u)/Cross(v,w);
    return P+v*t;
}

node get_c(node a,node b,node c) {
    node P=(a+b)/2.0;
    node Q=(a+c)/2.0;
    node v=rotate(b-a,Pi/2.0),w=rotate(c-a,Pi/2.0);
    if (dcmp(lenth(Cross(v,w)))==0) {
        if (dcmp(lenth(a-b)+lenth(b-c)-lenth(a-c))==0)
            return (a+c)/2.0;
        if (dcmp(lenth(a-c)+lenth(c-b)-lenth(a-b))==0)
            return (a+b)/2.0;
        if (dcmp(lenth(b-a)+lenth(a-c)-lenth(b-c))==0)
            return (b+c)/2.0;
    }
    return jiao(P,v,Q,w);
}

void min_circular() {
    random_shuffle(p+1,p+1+n);
    c=p[1]; r=0;
    for (int i=2;i<=n;i++) 
        if (dcmp(lenth(p[i]-c)-r)>0) {
            c=p[i],r=0;
            for (int j=1;j<i;j++)
                if (dcmp(lenth(p[j]-c)-r)>0) {
                    c=(p[i]+p[j])/2.0;
                    r=lenth(p[i]-c);
                    for (int k=1;k<j;k++) 
                        if (dcmp(lenth(p[k]-c)-r)>0) {
                            c=get_c(p[i],p[j],p[k]);
                            r=lenth(p[i]-c);
                        }
                }
        }
}

int main() {
    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
    min_circular();
    printf("%0.3lf\n",r);
    return 0;
}