最小圆覆盖利用随机增量法，是一种看似暴力但实际上很高效的算法。
概述
最小圆覆盖算法用来求：给定 $n$ 个点求覆盖所有点的最小的圆。
实现
首先我们对这 $n$ 个点随机打乱，可以利用 $\text{C++}$ 中的 $\text{random_shuffle}$ 函数。
首先令当前最小圆的圆心为点 $1$，半径为 $0$；接下来分成 $3$ 重循环。
- 枚举点 $i\in [1,n]$，如果点 $i$ 不在圆中则把圆心定为点 $i$ 并进入下层循环。
- 枚举点 $j\in[1,i)$，如果点 $j$ 不在圆中则用点 $i,j$ 构造圆并进入下层循环。
- 枚举点 $k\in[1,j)$，如果点 $k$ 不在圆中则用点 $i,j,k$ 构造圆（显然三点定圆），圆心使用原方程随便推一下就可以计算了。
时间复杂度：$\mathcal O(n)$，详见下方证明。
复杂度
上述过程看起来复杂度为 $\mathcal O(n^3)$，接下来进行严格分析。
由于这 $n$ 个点的最小覆盖圆只需要 $3$ 个点就可以确定了，因此每个点参与确定最小圆的概率为 $\frac{3}{n}$。所以第 $1,2$ 层循环中在第 $i$ 个点进入下一层循环的概率约为 $\frac{3}{i}$ 和 $\frac{2}{i}$。
设该算法的三个循环从外到内的复杂度分别为 $T_1(n)$；则有：
$$ \begin{aligned} T_3(n)&=\mathcal O(n) \\ T_2(n)&=\mathcal O(n)+\sum_{i=1}^n\frac{i-2}{i}\cdot \mathcal O(1)+\frac{2}{i}\cdot T_3(i)=\mathcal O(n) \\ T_1(n)&=\mathcal O(n)+\sum_{i=1}^n\frac{i-3}{i}\cdot \mathcal O(1)+\frac{3}{i}\cdot T_2(i)=\mathcal O(n) \end{aligned} $$
故整个算法的复杂度为 $\mathcal O(n)$！
代码
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <algorithm>
const int N = 1e6 + 5;
const double eps = 1e-12;
int n;
double r;
struct Point {
double x, y;
Point(double _x = 0, double _y = 0) {
x = _x, y = _y;
}
} p[N], o;
double sqr(double x) {
return x * x;
}
double dis(Point a, Point b) {
return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
bool cmp(double a, double b) {
return fabs(a - b) < eps;
}
bool inCir(Point a, Point o, double r) {
double d = dis(a, o);
return d < r || cmp(d, r);
}
void calc(Point a, Point b, Point c, Point &o, double &r) {
double a1 = 2 * (a.x - b.x), b1 = 2 * (a.y - b.y), c1 = sqr(a.x) + sqr(a.y) - sqr(b.x) - sqr(b.y);
double a2 = 2 * (a.x - c.x), b2 = 2 * (a.y - c.y), c2 = sqr(a.x) + sqr(a.y) - sqr(c.x) - sqr(c.y);
o.x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
o.y = (c1 * a2 - c2 * a1) / (b1 * a2 - b2 * a1);
r = dis(a, o);
}
void getCir(Point *p, int n, Point &o, double &r) {
srand(time(0) + *(new char));
std::random_shuffle(p + 1, p + n + 1);
o = p[1], r = 0;
for (int i = 1; i <= n; i++) {
if (inCir(p[i], o, r)) continue;
o = p[i], r = 0;
for (int j = 1; j < i; j++) {
if (inCir(p[j], o, r)) continue;
o.x = (p[i].x + p[j].x) / 2;
o.y = (p[i].y + p[j].y) / 2;
r = dis(o, p[i]);
for (int k = 1; k < j; k++) {
if (inCir(p[k], o, r)) continue;
calc(p[i], p[j], p[k], o, r);
}
}
}
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
scanf("%lf%lf", &p[i].x, &p[i].y);
}
getCir(p, n, o, r);
printf("%.10lf %.10lf %.10lf\n", o.x, o.y, r);
return 0;
}