随机增量法：bzoj 1336 && bzoj 1337 最小圆覆盖

1337: 最小圆覆盖

Time Limit: 1 Sec   Memory Limit: 64 MB
Submit: 1170   Solved: 573
[ Submit][ Status][ Discuss]

Description

给出平面上N个点,N<=10^5.请求出一个半径最小的圆覆盖住所有的点

Input

第一行给出数字N，现在N行，每行两个实数x,y表示其坐标.

Output

输出最小半径，输出保留三位小数.

Sample Input

4

1 0

0 1

0 -1

-1 0

Sample Output

1.000

题目：求出n个点的最小圆覆盖

随机增量法：

①将数组重新打乱（random_shuffle()函数）

②假设已经求出了前i个点的最小圆覆盖为Ci，检查第i+1个点是否在圆内，如果在，则C(i+1)=Ci并继续判定下一个点，不在的话，就令C(i+1)的圆心为第1个点和第i+1个点连线的中点，半径为连线长度的一半，并执行步骤③

③重新检查前i个点是否在圆C(i+1)上，如果某个点j不在，就用类似步骤②的方法继续修改圆心和半径使得j点在圆C(i+1)上，并执行步骤④

④再来一次一模一样的检查过程，不过这次只要有点k不在圆上，就令圆C(i+1)的圆心为三角形(i, j, k)的外心并扩增圆C(i+1)的半径

搞定，复杂度O(n^3)不过平均复杂度只有O(n)，证明略，应该很容易找到证明

#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
#define eps 1e-8
typedef struct Point
{
	double x;
	double y;
}Point;
Point s[100005];
Point operator + (const Point a, const Point b)
{  
	Point k;
	k.x = a.x+b.x, k.y = a.y+b.y;
	return k;
}
Point operator - (const Point a, const Point b)
{
	Point k;
	k.x = a.x-b.x, k.y = a.y-b.y;
	return k;
}
Point operator / (const Point a, double x)
{
	Point k;
	k.x = a.x/x, k.y = a.y/x;
	return k;
}
Point operator * (const Point a, double x)
{
	Point k;
	k.x = a.x*x, k.y = a.y*x;
	return k;
}

double Dist(Point a, Point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point Waixin(Point a, Point b, Point c)		//求三角形外接圆圆心（外心）
{
	Point ans;
	double a1, a2, b1, b2, c1, c2;
	a1 = 2*(b.x-a.x), b1 = 2*(b.y-a.y); c1 = (b.x*b.x)-(a.x*a.x)+(b.y*b.y)-(a.y*a.y);
	a2 = 2*(c.x-a.x), b2 = 2*(c.y-a.y); c2 = (c.x*c.x)-(a.x*a.x)+(c.y*c.y)-(a.y*a.y);
	if(fabs(a1)<eps)
		ans.y = c1/b1, ans.x = (c2-ans.y*b2)/a2;
	else if(fabs(b1)<eps)
		ans.x = c1/a1, ans.y = (c2-ans.x*a2)/b2;
	else
		ans.x = (c2*b1-c1*b2)/(a2*b1-a1*b2), ans.y = (c2*a1-c1*a2)/(b2*a1-b1*a2);
	return ans;
}
int main(void)
{
	Point P;
	double R;
	int n, i, j, k;
	scanf("%d", &n);
	for(i=1;i<=n;i++)
		scanf("%lf%lf", &s[i].x, &s[i].y);
	random_shuffle(s+1, s+n+1);
	P = s[1], R = 0;
	for(i=2;i<=n;i++)
	{
		if(Dist(P, s[i])-R>eps)
		{
			P = (s[1]+s[i])/2;
			R = Dist(P, s[i]);
			for(j=2;j<=i-1;j++)
			{
				if(Dist(P, s[j])-R>eps)
				{
					P = (s[i]+s[j])/2;
					R = Dist(P, s[i]);
					for(k=1;k<=j-1;k++)
					{
						if(Dist(P, s[k])-R>eps)
						{
							P = Waixin(s[i], s[j], s[k]);
							R = Dist(P, s[i]);
						}
					}
				}
			}
		}
	}
	printf("%.3f\n", R);
	return 0;
}