UVa 10095 Saving the Planet

题目描述

在科幻电影中，我们常常看到人类离开地球寻找新的栖息地。但在现实中，地球是我们唯一的家园，保护地球免受各种威胁是我们的神圣职责。现在，你将承担这一责任。

日本、澳大利亚、美国和俄罗斯的科学家发现了一个令人担忧的情况：一个巨大的未知形状和尺寸的物体正径直飞向地球。科学家们探测到一些来自分散位置的信号，并发现这些分散位置（可视为点）之间的相对距离保持不变。他们意识到，这些信号来自放置在该透明行星内部及表面的放射性元素。

科学家们正在制造武器来摧毁这颗半隐形的行星，但武器的威力取决于行星的大小。如果武器威力略小于摧毁行星所需的值，是可以接受的；但如果威力过大，地球本身将面临危险。

假设所有放射性元素都在行星内部，且行星为球体。你需要根据给定的放射性元素位置（三维空间中的点），确定可能的最小半径以及采样时球心的坐标。

输入格式

输入包含多组数据。每组数据第一行是一个整数 $n$ （ $n<10001$ ），随后 $n$ 行，每行是放射性元素的坐标 $(x,y,z)$。 $n=0$ 表示输入结束。

输出格式

每组数据输出一行四个浮点数：第一个是行星半径，后三个是球心的坐标。每个数保留四位小数，用空格分隔。假设 $\pi =3.14159265359$（虽然题目中给出了 $\pi$ ，但实际计算半径时并不需要用到它）。

样例输入

10
0.05 0.01 10.08
1.21 0.71 0.74
0.13 4.23 13.60
1.61 3.48 0.86
1.58 1.86 1.14
1.63 5.26 0.76
0.35 1.19 0.97
5.31 0.38 0.43
2.00 0.82 204.27
1.65 0.64 0.65
0

样例输出

101.9337 3.6550 0.6000 102.3500

问题分析

本题的核心问题是：给定三维空间中的 $n$ 个点，求能包含所有这些点的最小球体的半径和球心坐标。

这是一个经典的计算几何问题，称为 最小包围球（$\text{Minimum Enclosing Sphere}$）或 最小外接球 问题。它在许多领域都有应用，如碰撞检测、聚类分析、机器学习中的 $\text{Support Vector Machine}$ 等。

问题的数学形式

设点集为 $P=\{p_1,p_2,\dots ,p_n\}\subset {\mathbb{R}}^3$，我们需要找到球心 $c\in {\mathbb{R}}^3$ 和半径 $r\in \mathbb{R}$，使得：

等价地，我们需要最小化 $r$，使得所有点都在以 $c$ 为球心、$r$ 为半径的球内。

解的数学性质

最小包围球具有以下重要性质：

• 唯一性：对于任意点集，最小包围球是唯一的。
• 支撑点：最小包围球的边界上至少包含 $2$ 个点，最多包含 $4$ 个点（在三维空间中）。
• $\text{Karush-Kuhn-Tucker}$ 条件：球心 $c$ 位于其支撑点所构成的凸包内部（或边界上）。

这些性质为我们设计算法提供了理论基础。

解题思路

朴素方法

最直接的方法是枚举所有可能的球，但这是不可行的：可能的球心位置是连续的，枚举空间无限大。

常见算法比较

对于本题， $n$ 可达 $10000$，必须使用 $O(n)$ 或 $O(n\log n)$ 的算法，且需要保证解的精确性（保留四位小数）。

$\text{Welzl}$ 算法原理

$\text{Welzl}$ 算法是由 $\text{Emo Welzl}$ 于 $1991$ 年提出的随机增量算法，用于计算点集的最小包围球（也可推广到最小包围圆等问题）。其核心思想是 递归分治 + 随机化。

基本思路

设函数 $\text{welzl}(P,B)$ 表示：在给定点集 $P$ 和边界点集 $B$（这些边界点一定在最终球的边界上）的情况下，计算包含 $P$ 且以 $B$ 中所有点为边界点的最小球。

算法递归过程如下：

时间复杂度分析

• 随机化保证了点被选为边界点的概率很低。
• 期望情况下，每层的递归深度为 $O(\log n)$，且边界点集大小不超过 $4$。
• 总体期望时间复杂度为 $O(n)$。

几何子问题

在递归过程中，我们需要解决几个小规模几何问题：

• $1$ 个点定球：球心即为该点，半径为 $0$。
• $2$ 个点定球：球心为两点中点，半径为半距。
• $3$ 个点定球：计算三点的外接圆。需要注意三点共线的情况。
• $4$ 个点定球：计算四点的外接球。需要注意四点共面或共线的情况。

算法优化的关键

递归版本的 $\text{Welzl}$ 算法存在以下性能问题：

本题的通过代码采用了 非递归迭代版本 的 $\text{Welzl}$ 算法，具体优化如下：

• 使用三层循环代替递归，避免函数调用开销。
• 使用固定大小的数组存储边界点。
• 增量式构建球体：遇到球外点时，只用该点之前的点重新构建。
• 内层循环尽早检查点是否在球内，减少不必要的计算。

算法步骤详解

步骤 1：预处理

将输入的点集随机打乱。随机化保证了算法的时间复杂度期望为 $O(n)$，避免最坏情况的出现。

步骤 2：初始化

以第一个点作为初始球，半径为 $0$。

步骤 3：增量扫描

遍历所有点（ $i=1\dots n-1$ ）：

• 如果当前点 $p_i$ 在当前球内，继续下一个点。

• 否则，需要重新计算包含前 $i$ 个点的最小球：

  步骤 3.1：以 $p_i$ 为球心，半径为 $0$，重置边界点集。

  步骤 3.2：遍历 $j=0\dots i-1$ 的点 $p_j$：

  • 如果 $p_j$ 在当前球内，继续。
  • 否则，以 $p_i$ 和 $p_j$ 为直径构造球，边界点集设为 $\{p_i,p_j\}$。

  步骤 3.3：遍历 $k=0\dots j-1$ 的点 $p_k$：

  • 如果 $p_k$ 在当前球内，继续。
  • 否则，以 $p_i,p_j,p_k$ 三点构造外接球，边界点集设为 $\{p_i,p_j,p_k\}$。

  步骤 3.4：遍历 $l=0\dots k-1$ 的点 $p_l$：

  • 如果 $p_l$ 在当前球内，继续。
  • 否则，以 $p_i,p_j,p_k,p_l$ 四点构造外接球，边界点集设为 $\{p_i,p_j,p_k,p_l\}$。

步骤 4：几何计算

两点定球

给定点 $a$ 和 $b$，球心 $c=\left(\frac{x_a+x_b}{2},\frac{y_a+y_b}{2},\frac{z_a+z_b}{2}\right)$，半径 $r=\frac{\Vert a-b\Vert }{2}$。

三点定球

给定不共线的三点 $a,b,c$，求外接圆（球心在三点的平面上）。

设向量 $u=b-a$， $v=c-a$。我们需要找到参数 $u,v$ 使得：

$$c=a+uu+vv$$

且满足 $\Vert c-a\Vert =\Vert c-b\Vert =\Vert c-c\Vert$（即到三点等距）。

由 $\Vert c-a{\Vert }^2=\Vert c-b{\Vert }^2$ 可得：

$$\Vert uu+vv{\Vert }^2=\Vert (u-1)u+vv{\Vert }^2$$

展开并化简：

$$-2u\Vert u{\Vert }^2+\Vert u{\Vert }^2-2v(u\cdot v)=0$$

即：

$$2u\Vert u{\Vert }^2+2v(u\cdot v)=\Vert u{\Vert }^2$$

同理，由 $\Vert c-a{\Vert }^2=\Vert c-c{\Vert }^2$ 可得：

$$2u(u\cdot v)+2v\Vert v{\Vert }^2=\Vert v{\Vert }^2$$

解这个线性方程组即可得到 $u$ 和 $v$。

如果三点共线（ $\Vert u\times v\Vert \approx 0$ ），则退化为两点定球情况，取最远点对。

四点定球

给定四点 $a,b,c,d$，求外接球。设球心 $c$ 满足：

$$\Vert c-a{\Vert }^2=\Vert c-b{\Vert }^2=\Vert c-c{\Vert }^2=\Vert c-d{\Vert }^2$$

由 $\Vert c-a{\Vert }^2=\Vert c-b{\Vert }^2$ 可得：

$$2(b_x-a_x)x+2(b_y-a_y)y+2(b_z-a_z)z=b_x^2+b_y^2+b_z^2-a_x^2-a_y^2-a_z^2$$

类似地可以得到另外两个线性方程（对于 $c$ 和 $d$）。求解这个 $3\times 3$ 线性方程组即可得到球心坐标，半径即为球心到任意点的距离。

如果四点共面或共线，则退化为三点定球或两点定球的情况。

代码实现

// Saving the Planet
// UVa ID: 10095
// Verdict: Accepted
// Submission Date: 2026-05-28
// UVa Run Time: 0.090s
//
// 版权所有（C）2026，邱秋。metaphysis # yeah dot net

#include <bits/stdc++.h>
using namespace std;

const double EPS = 1e-8;

struct Point {
double x, y, z;
Point() : x(0), y(0), z(0) {}
Point(double x_, double y_, double z_) : x(x_), y(y_), z(z_) {}
};

double dist2(const Point& a, const Point& b) {
double dx = a.x – b.x, dy = a.y – b.y, dz = a.z – b.z;
return dx * dx + dy * dy + dz * dz;
}

double dist(const Point& a, const Point& b) {
return sqrt(dist2(a, b));
}

struct Sphere {
Point c;
double r;
Sphere() : c(Point()), r(–1) {}
Sphere(const Point& c_, double r_) : c(c_), r(r_) {}

bool contains(const Point& p) const {
return dist2(p, c) <= r * r + EPS;
}
};

// 由边界点构造最小球
Sphere sphereFromBoundary(Point* boundary, int m) {
if (m == 0) return Sphere(Point(0, 0, 0), –1);
if (m == 1) return Sphere(boundary[0], 0);
if (m == 2) {
Point c((boundary[0].x + boundary[1].x) / 2,
(boundary[0].y + boundary[1].y) / 2,
(boundary[0].z + boundary[1].z) / 2);
return Sphere(c, dist(boundary[0], boundary[1]) / 2);
}
if (m == 3) {
// 三点外接圆
double ax = boundary[1].x – boundary[0].x;
double ay = boundary[1].y – boundary[0].y;
double az = boundary[1].z – boundary[0].z;
double bx = boundary[2].x – boundary[0].x;
double by = boundary[2].y – boundary[0].y;
double bz = boundary[2].z – boundary[0].z;

double d1 = ax * ax + ay * ay + az * az;
double d2 = bx * bx + by * by + bz * bz;
double d3 = ax * bx + ay * by + az * bz;
double denom = d1 * d2 – d3 * d3;

if (fabs(denom) < EPS) {
// 共线，取最远两点
double d12 = dist2(boundary[0], boundary[1]);
double d13 = dist2(boundary[0], boundary[2]);
double d23 = dist2(boundary[1], boundary[2]);
if (d12 >= d13 && d12 >= d23) {
Point c((boundary[0].x + boundary[1].x) / 2,
(boundary[0].y + boundary[1].y) / 2,
(boundary[0].z + boundary[1].z) / 2);
return Sphere(c, sqrt(d12) / 2);
}
if (d13 >= d12 && d13 >= d23) {
Point c((boundary[0].x + boundary[2].x) / 2,
(boundary[0].y + boundary[2].y) / 2,
(boundary[0].z + boundary[2].z) / 2);
return Sphere(c, sqrt(d13) / 2);
}
Point c((boundary[1].x + boundary[2].x) / 2,
(boundary[1].y + boundary[2].y) / 2,
(boundary[1].z + boundary[2].z) / 2);
return Sphere(c, sqrt(d23) / 2);
}

double u = (d2 * d1 – d3 * d2) / (2 * denom);
double v = (d1 * d2 – d3 * d1) / (2 * denom);
Point c(boundary[0].x + u * ax + v * bx,
boundary[0].y + u * ay + v * by,
boundary[0].z + u * az + v * bz);
return Sphere(c, dist(c, boundary[0]));
}
// m == 4: 四点外接球
double a11 = boundary[1].x – boundary[0].x;
double a12 = boundary[1].y – boundary[0].y;
double a13 = boundary[1].z – boundary[0].z;
double a21 = boundary[2].x – boundary[0].x;
double a22 = boundary[2].y – boundary[0].y;
double a23 = boundary[2].z – boundary[0].z;
double a31 = boundary[3].x – boundary[0].x;
double a32 = boundary[3].y – boundary[0].y;
double a33 = boundary[3].z – boundary[0].z;

double b1 = (boundary[1].x * boundary[1].x + boundary[1].y * boundary[1].y + boundary[1].z * boundary[1].z –
boundary[0].x * boundary[0].x – boundary[0].y * boundary[0].y – boundary[0].z * boundary[0].z) / 2;
double b2 = (boundary[2].x * boundary[2].x + boundary[2].y * boundary[2].y + boundary[2].z * boundary[2].z –
boundary[0].x * boundary[0].x – boundary[0].y * boundary[0].y – boundary[0].z * boundary[0].z) / 2;
double b3 = (boundary[3].x * boundary[3].x + boundary[3].y * boundary[3].y + boundary[3].z * boundary[3].z –
boundary[0].x * boundary[0].x – boundary[0].y * boundary[0].y – boundary[0].z * boundary[0].z) / 2;

double det = a11 * (a22 * a33 – a23 * a32) – a12 * (a21 * a33 – a23 * a31) + a13 * (a21 * a32 – a22 * a31);

if (fabs(det) < EPS) {
// 退化，用三点球
return sphereFromBoundary(boundary, 3);
}

double x0 = (b1 * (a22 * a33 – a23 * a32) – a12 * (b2 * a33 – a23 * b3) + a13 * (b2 * a32 – a22 * b3)) / det;
double y0 = (a11 * (b2 * a33 – a23 * b3) – b1 * (a21 * a33 – a23 * a31) + a13 * (a21 * b3 – b2 * a31)) / det;
double z0 = (a11 * (a22 * b3 – b2 * a32) – a12 * (a21 * b3 – b2 * a31) + b1 * (a21 * a32 – a22 * a31)) / det;

Point c(boundary[0].x + x0, boundary[0].y + y0, boundary[0].z + z0);
return Sphere(c, dist(c, boundary[0]));
}

// Welzl 算法 – 非递归实现
Sphere welzl(vector<Point>& pts) {
int n = pts.size();
if (n == 0) return Sphere(Point(0, 0, 0), 0);
if (n == 1) return Sphere(pts[0], 0);

// 随机打乱
for (int i = n – 1; i > 0; —i) {
int j = rand() % (i + 1);
swap(pts[i], pts[j]);
}

Sphere s = sphereFromBoundary(&pts[0], 1);
s.c = pts[0];
s.r = 0;

Point boundary[4];
int bcnt = 1;
boundary[0] = pts[0];

for (int i = 1; i < n; ++i) {
if (s.contains(pts[i])) continue;

// 当前点在球外，需要重新计算
bcnt = 1;
boundary[0] = pts[i];
s = sphereFromBoundary(boundary, 1);

for (int j = 0; j < i; ++j) {
if (s.contains(pts[j])) continue;

bcnt = 2;
boundary[1] = pts[j];
s = sphereFromBoundary(boundary, 2);

for (int k = 0; k < j; ++k) {
if (s.contains(pts[k])) continue;

bcnt = 3;
boundary[2] = pts[k];
s = sphereFromBoundary(boundary, 3);

for (int l = 0; l < k; ++l) {
if (s.contains(pts[l])) continue;

bcnt = 4;
boundary[3] = pts[l];
s = sphereFromBoundary(boundary, 4);
}
}
}
}
return s;
}

int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout << fixed << setprecision(4);
srand(time(nullptr));

int n;
while (cin >> n && n) {
vector<Point> pts(n);
for (int i = 0; i < n; ++i)
cin >> pts[i].x >> pts[i].y >> pts[i].z;

Sphere s = welzl(pts);
cout << s.r << " " << s.c.x << " " << s.c.y << " " << s.c.z << "\\n";
}
return 0;
}

复杂度分析

• 时间复杂度：期望 $O(n)$，最坏情况 $O(n^2)$（但随机化保证了最坏情况几乎不会出现）。
• 空间复杂度： $O(n)$，主要用于存储输入点集。

对于 $n=10000$ ，期望运算次数约为 $10^5$ 量级，可以在 $1$ 秒内完成。

总结

本题是计算几何中的经典问题——最小包围球。 $\text{Welzl}$ 算法的随机增量思想非常巧妙，通过随机化将平均时间复杂度降为线性。非递归迭代实现避免了函数调用开销，进一步提升了运行效率。

掌握 $\text{Welzl}$ 算法不仅对解决本题有帮助，对理解其他随机增量算法（如最小包围圆、最小凸包等）也具有重要参考价值。