随机增量法

随机增量法

来自随机增量法 - OI Wiki

想法来自于洛谷一位未知人士，提问：最小覆盖圆问题（高中物理实验中的一个问题）。

其次，随机增量算法是随机和增量法的儿子。

增量法（Incremental Algorithm）

引入

随机增量算法是有关计算几何的重要算法，具有时间复杂度低（一般是线性的），且对理论知识要求不高，所以具有很广的应用范围。

类似于第一数学归纳法以后有时间学习一下，其本质是将一个问题化为规模刚好小一层的子问题。解决了子问题后加入当前对象。以下是递归式：

\[T(n) = T(n - 1) + g(n) \]

可以看出，其递归式十分的简洁，可以应用于很多几何题目中。

增量法加上随机化，可以避免很多坏的情况出现。

最小覆盖圆问题

题意

在一个平面上有 \(n\) 个点，求一个最小的圆，能够覆盖所有点

分析

假设，圆 \(O\) 是前 \(i - 1\) 个点的最小覆盖圆，那么，加入第 \(i\) 个点，如果在圆内或边上，则什么都不做。否则，新得到的最小覆盖圆肯定经过第 \(i\) 个点。

然后，以第 \(i\) 个点为基础（半径为 \(0\)），重复上述过程，加入第 \(j\) 个点，若第 \(j\) 个点在圆外，则最小覆盖圆必经过第 \(j\) 个点。

重复上述步骤。（因为需要三个点来确定这个最小覆盖圆，所以重复三次）

当我们遍历完所有点后，所得到的圆就是覆盖所有点的最小覆盖圆。

总结

1. 一开始只有一个点，构成一个半径为 \(0\) 的圆，考虑加入第 \(i\) 个点。
2. 如果这个新加的点在之前的最小圆中则跳过不管。
3. 否则我们尝试以 \(j \in [1,i - 1]\) 与 \(i\) 作为直径作圆。如果这个 \(j\) 已经在当前最小圆中，我们就什么都不做。
4. 否则就说明 \(k \in [1,j - 1]\) 有可能不在圆中，所以我们枚举这些点，然后以 \(i, j, k\) 的外接圆为新圆。

复杂度的讨论

时间复杂度：\(O(n)\)

空间复杂度：\(O(n)\)

关于时间复杂度的证明

看完上述分析，乍一看，复杂度是 \(O(n^3)\)，但是如果我们加入随机排列，期望的复杂度就会降到 \(O(n)\)。显然三个操作分别是 \(O(n)\) 的复杂度，最主要的是，会不会进入下一次循环。

我们考虑 \(1\) 到 \(2\) 的过程，三点定一个圆，所以一个点是答案的概率为 \(3 \over i\)。而 \(2\) 到 \(3\) 的过程是类似的，最后 \(3\) 的复杂度，应该是 \(O(j)\) 的。分析得，过程 \(2\) 的复杂度是 \(\sum_{j = 1} ^ i \space O(j) \times {3 \over j} = O(j)\)，所以，易得过程 \(3\) 的复杂度是 \(\sum_{i = 1} ^ n \space O(i) \times {3 \over i} = O(n)\)。所以总结，其时间复杂度为 \(O(n)\)。

一些细节

• 比较半径来判断点是否在圆内
• 两点中点作为圆心作圆
• 三点外接圆圆心是三角形的外心

代码

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
struct node {
	double x, y;
	friend node operator + (node a, node b) {return node{a.x + b.x, a.y + b.y};}
	friend node operator - (node a, node b) {return node{a.x - b.x, a.y - b.y};}
	friend node operator * (node a, double t) {return node{a.x * t, a.y * t};}
	friend node operator / (node a, double t) {return node{a.x / t, a.y / t};}
} p[100005], ans, now; int n; double ansr, nowr;
node rotate_90(node a) {return node{a.y, -a.x};}
node mid(node a, node b) {return node{(a.x + b.x) / 2, (a.y + b.y) / 2};}
double dist(node a, node b) {return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) *(a.y - b.y));}
bool check(node a, double r, node b) { 
	if (r - dist(a, b) >= -eps) return true; return false;
} double cross(node a, node b) {return a.x * b.y - a.y * b.x;}
node get_intersection(node p1, node v1, node p2, node v2) {
	double t = cross((p2 - p1), v2) / (cross(v1, v2)); return p1 + v1 * t;
} node get_circle(node a, node b, node c) {return get_intersection(mid(a, b), rotate_90(b - a), mid(a, c), rotate_90(c - a));}
mt19937 rnd(chrono::system_clock::now().time_since_epoch().count());
signed main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
    scanf("%lf%lf", &p[i].x, &p[i].y);
	shuffle(p + 1, p + 1 + n, rnd);
	ans = p[1]; ansr = 0;
	for (int i = 2; i <= n; i++) {
		if (check(ans, ansr, p[i])) continue;
		now = p[i]; nowr = 0;
		for (int j = 1; j <= i - 1; j++) {
			if (check(now, nowr, p[j])) continue;
			now = mid(p[i], p[j]);
			nowr = dist(p[i], p[j]) / 2;
			for (int k = 1; k <= j - 1; k++) {
				if (check(now, nowr, p[k])) continue;
				now = get_circle(p[i], p[j], p[k]);
				nowr = dist(now, p[i]);
			}
		}
		ans = now; ansr = nowr;
	} printf("%.10lf\n", ansr);
	printf("%.10lf %.10lf\n", ans.x, ans.y);
	return 0;
}