题解：P2533 [AHOI2012] 信号塔

· · 题解

方法一：随机增量法：

分析：用最小的圆覆盖住所有的点：随机增量法复杂度是 \mathcal O(n)；

性质：在既定的给定点条件下，如果引入一张新的半平面，只要此前的最优解顶点（即唯一确定最小包围圆的几个关键顶点）能够包含于其中，则不必对此最优解进行修改，亦即此亦为新点集的最优解；否则，新的最优解顶点必然位于这个新的半空间的边界上。定理可以通过反证法证明。于是，基于此性质，我们便可得到一个类似于线性规划算法的随机增量式算法。

定义 D_i 为相对于 p_i 的最小包围圆。此算法实现的关键在于对于 p_i \notin D_{i-1} 时的处理。显然，如果 p_i \notin D_{i-1}，则 D_i= D_{i-1}；否则，需要对 D_i 另外更新。而且，D_i 的组成必然包含了 p_i；因此，此种情况下的最小包围圆是过 p_i 点且覆盖点集 p_1，p_2，p_3 \cdots p_{i-1} 的最小包围圆。

则仿照上述处理的思路，D_i=p_1， p_i，逐个判断点集 p_2，p_3 \cdots \cdots p_{i-1}，如果存在 p_j \notin D_i，则 D_i=p_j，p_i 。同时，再依次对点集 p_1 ，p_2，p_3 \cdots \cdots p_{j-1} 判断是否满足 p_k \notin D_i，若有不满足，则 D_i=p_k，p_j，p_i。

由于，三点唯一地确定一个圆，故而，只需在此基础上判断其他的点是否位于此包围圆内，不停地更新 p_k。当最内层循环完成时，退出循环，转而更新 p_j；当次内层循环结束时，退出循环，更新 p_i。当 i=n 时，表明对所有的顶点均已处理过 ，此时的 D_n 即表示覆盖了给定 n 个点的最小包围圆。

总结： 假设圆 O 是前 i-1 个点得最小覆盖圆，加入第 i 个点，如果在圆内或边上则什么也不做。否，新得到的最小覆盖圆肯定经过第 i 个点。然后以第 i 个点为基础（半径为 0），重复以上过程依次加入第 j 个点，若第 j 个点在圆外，则最小覆盖圆必经过第 j 个点。重复以上步骤三次。遍历完所有点之后，所得到的圆就是覆盖所有点得最小圆。证明可以考虑这么做：最小圆必定是可以通过不断放大半径，直到所有以任意点为圆心，半径为半径的圆存在交点，此时的半径就是最小圆。所以上述定理可以通过这个思想得到。这个做法复杂度是 \mathcal O(n) 的，当加入圆的顺序随机时，因为三点定一圆，所以不在圆内概率是 \frac{3}{i}，求出期望可得是 \mathcal O(n)。

方法二：模拟退火法：

分析：对于每个枚举的点找到改点到所给点的最远点的距离，然后保证这个距离最小，即为所求圆的半径。

AC 代码（方法 1）：

#include <bits/stdc++.h>
using namespace std;
const int N = 1000005;
const double eps = 1e-8;
struct Point {
    double x, y;
} a[N], O;
double R;
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));}
void work(Point p1, Point p2, Point p3) {
    double a, b, c, d, e, f;
    a = p2.y - p1.y; b = p3.y - p1.y; c = p2.x - p1.x; d = p3.x - p1.x;
    f = p3.x * p3.x + p3.y * p3.y - p1.x * p1.x - p1.y * p1.y;
    e = p2.x * p2.x + p2.y * p2.y - p1.x * p1.x - p1.y * p1.y;
    O.x = (a * f - b * e) / (2 * a * d - 2 * b * c); O.y = (d * e - c * f) / (2 * a * d - 2 * b * c);
    R = dis(O, p1);
}
int main() {
    int n;
    scanf("%d", &n);
    for (int i=1; i<=n; ++i) scanf("%lf%lf", &a[i].x, &a[i].y);
    random_shuffle(a+1, a+n+1);
    O=a[1]; R=0;
    for (int i=2; i<=n; ++i)
        if (dis(a[i], O) > R + eps) {
            O = a[i]; R = 0;
            for (int j=1; j<i; ++j)
                if (dis(a[j], O) > R + eps) {
                    O.x = (a[i].x + a[j].x) / 2; O.y = (a[i].y + a[j].y) / 2; R = dis(a[j], O);
                    for (int k=1; k<j; ++k) if (dis(a[k], O) > R + eps) work(a[i], a[j], a[k]);
                }
        }
    printf("%.2lf %.2lf %.2lf\n", O.x, O.y, R);
    return 0;
}