前言
最小圆覆盖和最小矩形覆盖其实是有一定区别的。在矩形中有两组平行直线,在圆中则没有。
随机增量法
因为复杂度是在期望意义下计算的,所以无论题目给出什么样的数据,我们都要随机一下,这就是“随机”的含义。
我们考虑递推的过程。
假设我们现在已经求出来前 (i-1) 个点的最小覆盖圆。如果第 (i) 个点在这个圆内,那么这个点不产生影响。如果第 (i) 个点在圆外,我们考虑重构这个最小覆盖圆。
令 (i) 号点 (P_i) 为圆心,(r=0) 为半径,枚举 (j:1 ext{ to }i-1),重新进行递推。
假设我们现在已经求出来 (P_i) 和第 (1sim j-1) 个点的最小覆盖圆。如果第 (j) 个点在这个圆内,那么这个点不产生影响。如果第 (j) 个点在圆外,那么同样重构这个最小覆盖圆。但是此时已经钦定了 (i) 号点在圆上,所以设 (i) 号点和 (j) 号点的中点为 (M)。
令 (M) 点为圆心,(r=frac{left|P_iP_j ight|}2) 为半径,枚举 (k:1 ext{ to }j-1),进一步递推。
现在我们已经知道了 (P_i,P_j) 一定在圆上。假设我们已经求出来 (P_i,P_j) 和 (1sim k-1) 号点的最小覆盖圆。如果 (k) 不在圆内,直接更新答案(此时是过 (P_i) 和 (1sim j) 号的点的最小覆盖圆)为 (P_i,P_j,P_k) 的外接圆。
过 (P_i) 和 (1sim j) 号点的最小覆盖圆找到了,我们下一步就可以做 (1sim j+1) 了,返回第 5 段,进一步进行。
综上所述本题一共用了 3 层迭代,最终通过递推的方式找出了最小覆盖圆。
复杂度分析
但是这样的三重循环看样子是 (Oleft(n^3 ight)) 的啊。
我们现在回到一开始的随机操作,它的意义是让数据尽可能达到期望意义。
对于过 (P_i,P_j) 的圆的枚举,它全面地涉及到了 (1sim j-1) 中的每个点,至少循环了一遍,这一层复杂度是 (O(j))。
对于过 (P_i) 的圆的枚举,它要枚举 (1sim i-1) 中的那些点,考虑到三点定圆,覆盖到前 (i) 个点的最小覆盖圆也是由三个点定下来的,因此在前 (i-1) 号点中,期望有两个点做出贡献。因此能迭代进下一层的只有 (Oleft(sum_{j=1}^ifrac 2j ight)) 个。
同理,对于全局最小圆覆盖,在前 (i) 号点中,期望有三个点能做出贡献。能迭代进下一层的只有 (Oleft(sum_{i=1}^nfrac 3i ight)) 个。
总复杂度为 (Oleft(sum_{i=1}^nfrac 3icdot sum_{j=1}^ifrac 2jcdot j ight)=O(n))。
代码
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define db double
#define eps 1e-12
struct pnt
{
db x,y;
pnt(db x,db y)
{
this->x=x;
this->y=y;
}
pnt(){}
friend pnt operator +(pnt a,pnt b)
{return pnt(a.x+b.x,a.y+b.y);}
friend pnt operator -(pnt a,pnt b)
{return pnt(a.x-b.x,a.y-b.y);}
friend db operator *(pnt a,pnt b)
{return a.x*b.y-a.y*b.x;}
db dis()
{return sqrt(x*x+y*y);}
pnt times(db k)
{return pnt(k*x,k*y);}
}p[100100];
pnt its(pnt a,pnt b,pnt c,pnt d)//ab和cd交点
{
db u=(c-a)*(d-a);
db D=u+(d-b)*(c-b);
return a+(b-a).times(u/D);
}
int main()
{
int n;
scanf("%d",&n);
for(int i=1;i<=n;++i)
scanf("%lf%lf",&p[i].x,&p[i].y);
std::random_shuffle(p+1,p+1+n);
pnt c(0,0);
db r=0.0;
for(int i=1;i<=n;++i)
{
if((p[i]-c).dis()<r+eps)
continue;
c=p[i];
r=0.0;
for(int j=1;j<i;++j)
{
if((p[j]-c).dis()<r+eps)
continue;
r=(p[j]-p[i]).dis()/2;
c=p[i]+(p[j]-p[i]).times(.5);
for(int k=1;k<j;++k)
{
if((p[k]-c).dis()<r+eps)
continue;
pnt a=p[j]-p[i];
pnt f=p[i]+a.times(.5);
pnt b=p[k]-p[j];
pnt g=p[j]+b.times(.5);
c=its(f,f+pnt(-a.y,a.x),g,g+pnt(-b.y,b.x));
r=(c-p[i]).dis();
}
}
}
printf("%.10lf
%.10lf %.10lf
",r,c.x,c.y);
return 0;
}