1. 程式人生 > >bzoj1336&1337 最小圓覆蓋

bzoj1336&1337 最小圓覆蓋

1336: [Balkan2002]Alien最小圓覆蓋

Time Limit: 1 Sec  Memory Limit: 162 MBSec  Special Judge
Submit: 1473  Solved: 648
[Submit][Status][Discuss]

Description

給出N個點,讓你畫一個最小的包含所有點的圓。

Input

先給出點的個數N,2<=N<=100000,再給出座標Xi,Yi.(-10000.0<=xi,yi<=10000.0)

Output

輸出圓的半徑,及圓心的座標

Sample Input

6
8.0 9.0
4.0 7.5
1.0 2.0
5.1 8.7
9.0 2.0
4.5 1.0


Sample Output

5.00
5.00 5.00

經典模型,隨機增量法求最小圓覆蓋

假如已經求出了前i-1個點的最小圓覆蓋,假如第i個點後,如果不在所求的圓內,那麼這個點一定在新圓的邊界上。這樣我們只要列舉另外兩個點就可以了,因為三點確定一個圓。

列舉方法為:從1-(i-1)列舉一個不在圓內的點j,再從1-(j-1)列舉一個不在圓內的點k,那麼i j k三點就可以確定一個圓。(注意三點共線的特殊情況)

這樣的列舉方法一定可以保證列舉到所有的點對,從而保證所有點在圓內且圓最小。(傻傻的我想這個想了好久……)

這樣的列舉複雜度看似是O(n^3),但在順序隨機的情況下可以證明覆雜度是O(n)。


#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#define F(i,j,n) for(int i=j;i<=n;i++)
#define D(i,j,n) for(int i=j;i>=n;i--)
#define ll long long
#define maxn 100005
#define eps 1e-8
using namespace std;
int n;
double r;
struct P
{
	double x,y;
	friend P operator +(P a,P b){return (P){a.x+b.x,a.y+b.y};}
	friend P operator -(P a,P b){return (P){a.x-b.x,a.y-b.y};}
	friend P operator /(P a,double b){return (P){a.x/b,a.y/b};}
}p[maxn],t;
inline double dis(P a,P b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
inline P rev(P a)
{
	return (P){a.y,-a.x};
}
inline P centre(P a,P b,P c)
{
	if (fabs((a.x-b.x)*(a.y-c.y)-(a.x-c.x)*(a.y-b.y))<=eps)
	{
		if (a.x<b.x) swap(a,b);
		if (a.x<c.x) swap(a,c);
		if (b.x<c.x) swap(b,c);
		return (a+c)/2;
	}
	P x1=(a+b)/2,x2=x1+rev(a-b),y1=(a+c)/2,y2=y1+rev(a-c);
	if (fabs(y1.x-y2.x)<eps) swap(x1,y1),swap(x2,y2);
	double k2=(y1.y-y2.y)/(y1.x-y2.x),b2=y2.y-y2.x*k2;
	if (fabs(x1.x-x2.x)<eps) return (P){x1.x,k2*x1.x+b2};
	double k1=(x1.y-x2.y)/(x1.x-x2.x),b1=x2.y-x2.x*k1;
	double x=(b2-b1)/(k1-k2);
	return (P){x,k1*x+b1};
}
int main()
{
	scanf("%d",&n);
	F(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y);
	random_shuffle(p+1,p+n+1);
	t=p[1];r=0.0;
	F(i,2,n) if (dis(t,p[i])>r+eps)
	{
		t=(p[1]+p[i])/2;r=dis(p[i],t);
		F(j,2,i-1) if (dis(t,p[j])>r+eps)
		{
			t=(p[i]+p[j])/2;r=dis(p[i],t);
			F(k,1,j-1) if (dis(t,p[k])>r+eps)
			{
				t=centre(p[i],p[j],p[k]);
				r=dis(p[i],t);
			}
		}
	}
	printf("%.10lf\n%.10lf %.10lf\n",r,t.x,t.y);
	return 0;
}