[模板] 最小圆覆盖

一、题目

点此看题

二、解法

定理：如果 \(p\) 不在集合 \(S\) 的最小圆覆盖中，那么它一定在集合 \(S\cup\{p\}\) 的最小圆覆盖上。

多次运用此定理即可，我们先把所有点 \(\tt random\_shuffle\) 一遍，然后维护前 \(i\) 个点的最小圆覆盖。如果 \(i\) 不在前 \(i-1\) 个点的最小圆覆盖中，那么它一定在圆上，所以我们以 \(i\) 作为圆的定点之一，暂时以它为圆心考虑让前 \(i-1\) 个点在新的圆中。

\(j\) 从 \(1\) 开始扫，如果 \(j\) 在圆内则跳过；如果 \(j\) 不在圆内，那么它一定是圆的定点之一（考虑 \([1,j]\cup {i}\) 这些点的情况下），这样两点可以确定一个半径最小的圆。再用 \(k\) 从 \(1\) 开始扫，如果 \(k\) 在圆内则跳过；如果 \(k\) 不在圆内，那么它一定是圆的定点之一，这样可以三点定圆，写一个求线段交点即可解决。

我们管它叫随机增量法，看上去时间复杂度是 \(O(n^3)\) 的，实际上每个点只有 \(\frac{3}{|S|}\) 的概率作为圆上定点，才能进入下一重循环。我们设每重循环的时间复杂度分别为 \(T_1(n),T_2(n),T_3(n)\)，可以列出方程：

\[T_3(n)=O(n) \]

\[T_2(n)=O(n)+\sum_{i=1}^n\frac{3}{i}T_3(i)=O(n) \]

\[T_1(n)=O(n)+\sum_{i=1}^n\frac{3}{i}T_2(i)=O(n) \]

所以总时间复杂度 \(O(n)\)，但是注意要 \(\tt srand(time(0))\) 不然会被卡。

#include <cstdio>
#include <algorithm>
#include <ctime>
#include <cmath>
using namespace std;
#define db double
const int M = 100005;
int read()
{
	int x=0,f=1;char c;
	while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
	while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
	return x*f;
}
int n;
struct point
{
	db x,y;
	point(db X=0,db Y=0) : x(X) , y(Y) {}
	point operator + (db t) {return point(x+t,y+t);}
	point operator - (db t) {return point(x-t,y-t);}
	point operator * (db t) {return point(x*t,y*t);}
	point operator / (db t) {return point(x/t,y/t);}
	point operator + (point r) {return point(x+r.x,y+r.y);}
	point operator - (point r) {return point(x-r.x,y-r.y);}
	point rotate() {return point(y,-x);}
	db len() {return x*x+y*y;}
}p[M];
db cdt(point a,point b) {return a.x*b.x+a.y*b.y;}
db crs(point a,point b) {return a.x*b.y-a.y*b.x;}
point get(point a,point l0,point b,point l1)
{
	db k=crs(b-a,l1)/crs(l0,l1);
	return a+l0*k;
}
point cir(point a,point b,point c)
{
	return get((a+b)/2,(b-a).rotate(),(b+c)/2,(c-b).rotate());
}
signed main()
{
	n=read();srand(time(0));
	for(int i=1;i<=n;i++)
		scanf("%lf %lf",&p[i].x,&p[i].y);
	random_shuffle(p+1,p+1+n);
	point o;db r=0;
	for(int i=1;i<=n;i++)
	{
		if((p[i]-o).len()>r)
		{
			o=p[i];r=0;
			for(int j=1;j<i;j++)
				if((o-p[j]).len()>r)
				{
					o=(p[i]+p[j])/2;r=(o-p[j]).len();
					for(int k=1;k<j;k++)
						if((o-p[k]).len()>r)
						{
							o=cir(p[i],p[j],p[k]);
							r=(o-p[k]).len();
						}
				}
		}
	}
	printf("%.10f\n%.10f %.10f\n",sqrt(r),o.x,o.y);
}