最小圆覆盖

最小圆覆盖

期望复杂度 (mathcal O(n)) 求解面积最小的覆盖 (n) 个点的圆，采用的是随机增量法。随机增量法类似归纳法，但由于随机性能够保证复杂度优秀。

算法思路：

• 假设我们找到了能够包含前 (i-1) 个点的最小圆，如果第 (i) 个点在圆中，则前 (i) 个点的最小圆不变
• 否则最小圆一定经过点 (i)，接着我们继续采用这样的思路找出包含前 (i-1) 个点的圆，且满足圆经过点 (i)，这一步最终可以归纳到所有点的最小圆
• 假设满足前 (j-1) 个点在过 (i) 的最小圆上，考虑点 (j)，如果在圆中，则圆不变，否则一定经过点 (j)，这一步最终可以归纳到前 (i) 个点的最小圆
• 此时对于前 (j) 个点加上 (i) 号点，我们已经知道必定经过 (i) 和 (j)，只需要枚举前 (j-1) 个点做为圆上第三个点。同样地，对于前 (k-1) 个点和 (i,j) 的最小圆，如果 (k) 在里面则圆不变，否则一定经过 (i,j,k)，更换圆即可，这里最终能得到前 (j) 个点和 (i) 号点的最小圆。

所以我们这样做

• 第一层循环判断点 (i) 是否在现有的圆中（初始时圆为空）
• 如果不在，接着第二层循环枚举第二个点 (j(j<i))
• 以这两个点为直径的圆判断点 (k(k<j)) 是否在园里，如果不在，选取以 (i,j,k) 构成的圆

当然，请注意三点共线的情况。模板题没有考虑。

分析下复杂度。显然最里面一层循环复杂度是 (mathcal O(j))。

第二层循环当且仅当 (j) 不在前 (j-1) 和 (i) 构成的最小圆中，会进入第三层。这里期望是 (mathcal O(dfrac3j))，因为仅仅三个点（或两个点）确定最小圆，(j) 只有是这三个点才会扩大圆，进入第三层。所以第二层复杂度是 (mathcal O(i))。

同样的道理，第一层复杂度为 (mathcal O(n))。

#include <bits/stdc++.h>
using std::swap;
typedef double DB;
const int N = 1000005;
const DB eps = 1e-6;
int n;
struct Point {
	DB x, y;
	Point(DB x = 0, DB y = 0) : x(x), y(y) {}
	Point operator + (Point a) { return Point(x + a.x, y + a.y); }
	Point operator - (Point a) { return Point(x - a.x, y - a.y); }
	friend Point operator * (DB k, Point a) { return Point(k * a.x, k * a.y); }
	operator DB() { return sqrt(x * x + y * y); }
} a[N], O;
DB R;
Point work(Point P1, Point P2, Point P3) {
	DB a, b, c, d, e, f, x, y;
	a = P2.x - P1.x, b = P2.y - P1.y, c = P3.x - P1.x, d = P3.y - P1.y;
	e = P1.x * P1.x + P1.y * P1.y - P2.x * P2.x - P2.y * P2.y; e /= 2;
	f = P1.x * P1.x + P1.y * P1.y - P3.x * P3.x - P3.y * P3.y; f /= 2;
	x = (b * f - d * e) / (a * d - b * c), y = (c * e - a * f) / (a * d - b * c);
	return Point(x, y);
}
int main() {
	srand(time(0));
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
		scanf("%lf%lf", &a[i].x, &a[i].y);
	for (int i = 1; i <= n; i++)
		swap(a[i], a[rand() % n + 1]);
	for (int i = 1; i <= n; i++)
		if ((DB)(a[i] - O) > R + eps) {
			O = a[i], R = 0;
			for (int j = 1; j < i; j++)
				if ((DB)(a[j] - O) > R + eps) {
					O = 0.5 * (a[i] + a[j]), R = (DB)(a[i] - O);
					for (int k = 1; k < j; k++)
						if ((DB)(a[k] - O) > R + eps) O = work(a[i], a[j], a[k]), R = (DB)(a[i] - O);
				}
		}
	printf("%.2lf %.2lf %.2lf
", O.x, O.y, R);
	return 0;
}