最小圆覆盖

最小圆覆盖

期望复杂度 \(\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\n", O.x, O.y, R);
	return 0;
}