最小圆覆盖

HDU3007

P1742 最小圆覆盖

P2533 [AHOI2012]信号塔

这三道题都是最小圆覆盖问题，只不过数组或者输出形式有点变化而已

随机增量法是一个可以在 期望 O(n)时间内求出最小圆覆盖的算法，首先它的算法流程是这样的

枚举第一个点 i，若不在目前圆内，设它为圆心

枚举第二个点 j，若不在当前圆内，设当前圆为以 i,j 为直径的圆

枚举第三个点 k，若不在当前圆内，设当前圆为 i,j,k 的外接圆

正确性

显然最优解一定是两个点为直径的圆或者一个三角形的外接圆，否则肯定能缩的更小。那么这么枚举的正确性是比较显然的了

时间复杂度

这是一个重点，这么做看似是O(n^3) 的，不过对于 随机顺序 的点，是可以 期望 O(n) 的。下面考虑证明：

显然，最后一层循环枚举从 1~j，只要进入循环就一定要跑完，所以是O(j) 的

考虑倒数第二层循环，什么情况下会进入第三层循环呢？仅当 j 不在前 j-1 个点形成的圆中，考虑 j 个点形成的圆是由三个点确定的，那么第 j 个 (最后一个点) 若是三个点之一，则需要扩大圆，否则不需要进入第三层循环，这个概率是3/j 的，所以第二层的复杂度是O(i) 的同理，第一层的复杂度就是O(n)的了

当然，这是基于随机数据的期望复杂度，所以我们一般需要手动打乱排列

附上代码：

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

using namespace std;

const int M=1e5+5;
const double eps=1e-10;

double dcmp(double x)
{
    if(fabs(x)<eps){
        return 0;
    }
    return x;
}

int n;
struct point{
    double x,y;
    point(){}
    point(double _x,double _y){x=_x,y=_y;}
    point operator + (point p){
        return point(x+p.x,y+p.y);
    }
    point operator - (point p){
        return point(x-p.x,y-p.y);
    }
    point operator * (double k){
        return point(x*k,y*k);
    }
    point operator / (double k){
        return point(x/k,y/k);
    }
}pp[M];

double cross(point a,point b)
{
    return a.x*b.y-a.y*b.x;
}

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

point getmid(point a,point b)
{
    return point((a.x+b.x)/2,(a.y+b.y)/2);
}

point rotate(point a)
{
    return point(-a.y,a.x);
}

struct circle{
    point o;
    double r;
    circle(){}
    circle(point a,double r):o(a),r(r){}
};

circle circlein(point t1,point t2,point t3)
{
    point v1=t2-t1,v2=t1-t3;
    v1=rotate(v1),v2=rotate(v2);
    point p1=getmid(t1,t2),p2=getmid(t1,t3);
    point u=p1-p2;
    double t=cross(v2,u)/cross(v1,v2);
    point oo=p1+v1*t;
    return circle(oo,dcmp(dis(oo,t1)));
}

void work()
{
    circle ans;
    ans.o=pp[1];ans.r=0;
    for(int i=2;i<=n;i++){
        if(dis(pp[i],ans.o)>ans.r+eps){
            ans=circle(pp[i],0);
            for(int j=1;j<i;j++){
                if(dis(pp[j],ans.o)>ans.r+eps){
                    ans.o=getmid(pp[i],pp[j]);
                    ans.r=dis(pp[j],pp[i])/2;
                    for(int k=1;k<j;k++){
                        if(dis(pp[k],ans.o)>ans.r+eps){
                            ans=circlein(pp[i],pp[j],pp[k]);
                        }
                    }
                }
            }
        }
    }
    //printf("%.2lf %.2lf %.2lf\n",ans.o.x,ans.o.y,ans.r);
    printf("%.10lf\n%.10lf %.10lf\n",ans.r,ans.o.x,ans.o.y);
}

int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&pp[i].x,&pp[i].y);
    }
    random_shuffle(pp+1,pp+n+1);
    work();
    return 0;
}