最小圆覆盖（随机增量||模拟退火）

最小圆覆盖

https://www.lydsy.com/JudgeOnline/problem.php?id=1337

Time Limit: 1 Sec  Memory Limit: 64 MB
Submit: 1464  Solved: 737
[Submit][Status][Discuss]

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

 

ac1() 模拟退火，ac2()随机增量法

#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
const double eps=1e-8;
const double INF=1e20;
const double PI=acos(-1.0);


///pick定理：面积=内部格点数+（边上格点数/2）-1

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

inline double sqr(double x){return x*x;}
inline double hypot(double a,double b){return sqrt(a*a+b*b);}
struct Point{
    double x,y;
    Point(){}
    Point(double _x,double _y){
        x=_x;
        y=_y;
    }
    void input(){
        scanf("%lf %lf",&x,&y);
    }
    Point operator - (const Point &b)const{
        return Point(x-b.x,y-b.y);
    }
    //叉积
    double operator ^ (const Point &b)const{
        return x*b.y-y*b.x;
    }
    //点积
    double operator * (const Point &b)const{
        return x*b.x+y*b.y;
    }
    //返回长度
    double len(){
        return hypot(x,y);
    }
    //返回两点的距离
    double distance(Point p){
        return hypot(x-p.x,y-p.y);
    }
    Point operator + (const Point &b)const{
        return Point(x+b.x,y+b.y);
    }
    Point operator * (const double &k)const{
        return Point(x*k,y*k);
    }
    Point operator / (const double &k)const{
        return Point(x/k,y/k);
    }

    //逆时针转90度
    Point rotleft(){
        return Point(-y,x);
    }
};

struct Line{
    Point s,e;
    Line(){}
    Line(Point _s,Point _e){
        s=_s;
        e=_e;
    }
    //根据一个点和倾斜角angle确定直线，0<=angle<pi
    Line(Point p,double angle){
        s=p;
        if(sgn(angle-PI/2)==0){
            e=(s+Point(0,1));
        }
        else{
            e=(s+Point(1,tan(angle)));
        }
    }
    //ax+by+c=0;
    Line(double a,double b,double c){
        if(sgn(a)==0){
            s=Point(0,-c/b);
            e=Point(1,-c/b);
        }
        else if(sgn(b)==0){
            s=Point(-c/a,0);
            e=Point(-c/a,1);
        }
        else{
            s=Point(0,-c/b);
            e=Point(1,(-c-a)/b);
        }
    }

    //求两直线的交点
    //要保证两直线不平行或重合
    Point crosspoint(Line v){
        double a1=(v.e-v.s)^(s-v.s);
        double a2=(v.e-v.s)^(e-v.s);
        return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
    }
};

struct circle{
    Point p;
    double r;
    circle(){}
    circle(Point _p,double _r){
        p=_p;
        r=_r;
    }

    circle(double x,double y,double _r){
        p=Point(x,y);
        r=_r;
    }

    circle(Point a,Point b,Point c){///三角形的外接圆
        Line u=Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
        Line v=Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
        p=u.crosspoint(v);
        r=p.distance(a);
    }


};

Point p[100005];
int n;

void ac2(){
    random_shuffle(p+1,p+n+1);
    circle c(p[1],0);
    for(int i=1;i<=n;i++){
        if(sgn((c.p-p[i]).len()-c.r)>0){
            c.p=p[i],c.r=0;
            for(int j=1;j<i;j++){
                if(sgn((c.p-p[j]).len()-c.r)>0){
                    c.p=(p[i]+p[j])/2;
                    c.r=(c.p-p[j]).len();
                    for(int k=1;k<j;k++){
                        if(sgn((c.p-p[k]).len()-c.r)>0){
                            c=circle(p[i],p[j],p[k]);
                        }
                    }
                }
            }
        }
    }
    printf("%f\n%f %f\n",c.r,c.p.x,c.p.y);
}

double dist(Point a,Point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

void ac1(){

    double ans=1e99;
    Point tmp;
    tmp.x=tmp.y=0;
    int s=1;
    double step=10000;
    double esp=0.000001;
    while(step>esp){
        for(int i=1;i<=n;i++){
            if(dist(tmp,p[s])<dist(tmp,p[i])) s=i;
        }
        double Dist=dist(tmp,p[s]);
        ans=min(ans,Dist);
        tmp.x+=(p[s].x-tmp.x)/Dist*step;
        tmp.y+=(p[s].y-tmp.y)/Dist*step;
        step*=0.98;
    }
    printf("%.3f\n",ans);
}
int main(){

    scanf("%d",&n);
        for(int i=1;i<=n;i++){
            p[i].input();
        }
        ac1();

    return 0;
}