1. 程式人生 > 其它 >[模板] 最小圓覆蓋

[模板] 最小圓覆蓋

一、題目

點此看題

二、解法

定理:如果 \(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);
}