【luogu P1742】最小圆覆盖（随机增量法）（计算几何）

最小圆覆盖

题目链接：luogu P1742

题目大意

给你二维平面上的一些点，要你找到一个半径最小的圆覆盖所有的点。

思路

这个有随机增量法和模拟退火。
（其实怎么看都是模拟退火好写，但随机增量法至少能保证绝对正确，故趁着二维还好就写了）
（不过三维和多维的模拟退火也不一定能过）

首先有一个引理：如果 \(p\) 不在集合 \(S\) 的最小圆覆盖上，那么 \(p\) 一定在集合 \(S\cup p\) 的最小圆覆盖上。
（这个在高维也适用）

所以就有一个方法：
我们随便打乱点的顺序，然后按顺序依次选一个点，此时的一个点只能确定圆心在自己，半径为 \(0\) 的圆。
然后接着再前面的点中选一个点，那两个点就是圆心在这两个点西那段的中点，半径直接随便选一个点跟中点算一下即可（后面也一样）。
然后再在这两个的前面选一个点，那这个时候就是找这个三角形的外心，可以用几何的方法，也可以用高斯消元（见某球形什么的题）。

然后要注意的是因为上面的引理，所以你是要找到不在当前圆里面的才需要搞这些操作，三个循环都是。

接着要证明它的复杂度其实是 \(O(n)\)。（非常玄学的证明方法）
因为我们打乱了它的顺序，所以它是随机的，就是一个随机的增量。
所以第 \(i\) 个点不在那个圆里面的概率是 \(\dfrac{3}{i}\)。

所以要修改的次数就是：
\(\sum\limits_{i=1}^n(\dfrac{3}{i}\times\sum\limits_{j=1}^n(\dfrac{3}{j}\times j))=9n\)
所以就是 \(O(n)\) 的了。

代码

#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>

using namespace std;

struct node {
	long double x, y;
}a[100001], ans;
int n;
long double r2, f[12][11], fc[12][11];

node operator +(node x, node y) {
	return (node){x.x + y.x, x.y + y.y};
}

node operator -(node x, node y) {
	return (node){x.x - y.x, x.y - y.y};
}

node operator /(node x, long double y) {
	return (node){x.x / y, x.y / y};
}

long double dis2(node x) {
	return x.x * x.x + x.y * x.y;
}

long double Abs(long double x) {
	return x < 0 ? -x : x;
}

void gaosi(int n) {
	for (int i = 1; i <= n; i++) {
		int maxx = i;
		for (int j = i + 1; j <= n; j++)
			if (Abs(fc[j][i]) > Abs(fc[maxx][i]))
				maxx = i;
		swap(fc[i], fc[maxx]);
		for (int j = 1; j <= n; j++) {
			if (i == j) continue;
			long double tmp = fc[j][i] / fc[i][i];
			for (int k = i + 1; k <= n + 1; k++) {
				fc[j][k] -= fc[i][k] * tmp;
			}
		}
	}
}

node get_cir(node x, node y, node z, int n) {
	f[1][1] = x.x; f[1][2] = x.y;
	f[2][1] = y.x; f[2][2] = y.y;
	f[3][1] = z.x; f[3][2] = z.y;
	
	memset(fc, 0, sizeof(fc));
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++) {
			fc[i][j] = 2 * (f[i][j] - f[i + 1][j]);
			fc[i][n + 1] += f[i][j] * f[i][j] - f[i + 1][j] * f[i + 1][j];
		}
	
	gaosi(n);
	
	node re;
	re.x = fc[1][n + 1] / fc[1][1];
	re.y = fc[2][n + 1] / fc[2][2];
	return re;
}

int main() {
	srand(19491001);
	
	scanf("%d", &n);
	for (int i = 1; i <= n; i++) {
		scanf("%Lf %Lf", &a[i].x, &a[i].y);
		rand();
	}
	
	random_shuffle(a + 1, a + n + 1);
	
	for (int i = 1; i <= n; i++) {
		if (dis2(a[i] - ans) > r2) {
			ans = a[i]; r2 = 0;
			for (int j = 1; j < i; j++) {
				if (dis2(a[j] - ans) > r2) {
					ans = (a[i] + a[j]) / 2; r2 = dis2(a[j] - ans);
					for (int k = 1; k < j; k++) {
						if (dis2(a[k] - ans) > r2) {
							ans = get_cir(a[i], a[j], a[k], 2); r2 = dis2(a[k] - ans);
						}
					}
				}
			}
		}
	}
	
	printf("%.10Lf\n%.10Lf %.10Lf", sqrt(r2), ans.x, ans.y);
	
	return 0;
}