P2533 题解（随机增量法详解）

· · 题解

题目大意

给定 n 个点，求能覆盖这 n 个点的最小的圆。

题目分析

本题主要用到随机增量法，本蒟蒻将于本文详细介绍随机增量法。

随机增量算法是计算几何的一个重要算法，它对理论知识要求不高，算法时间复杂度低，应用范围广大。

增量法（Incremental Algorithm）

增量法，顾名思义就是先确定值，再将其增量的算法。

增量法的思想与第一数学归纳法类似，它的本质是将一个问题化为规模刚好小一层的子问题。解决子问题后加入当前的对象。

增量法形式简洁，可以应用于许多的几何题目中。

增量法往往结合随机化，可以避免最坏情况的出现。

希尔排序实际上就是一种缩小增量法。

最小圆问题

首先，我们不难发现，最后的圆一定是两点所连线段为直径的圆或三点所连成的三角形的外接圆。于是，我们考虑一个一个加入这些点后最小圆覆盖的变化。

一开始只有第一个点构成一个半径为 0 的圆，考虑加入第 i 个点，按下面步骤实施：

  1. 如果这个点在之前的最小圆中则跳过不管。
  2. 我们加入点 \forall j \in [1, i-1]，如果这个点已经在当前最小圆中我们不更新。否则以 j 与 i 所连直径为直径画圆。
  3. 但这时 \forall k \in [1, j-1] 有可能有出圆了，所以我们枚举这些点，然后以 \Delta ijk 的外接圆为新圆。

最后得到的圆即为答案。然后考虑时间复杂度。

乍一看复杂度是 \mathcal O(n ^ 3)，但是我们如果给这些点随机排列一下，再来仔细分析一下复杂度。

显然三个操作分别是 \mathcal O(n) 的时间复杂度，最主要的部分在于会不会进入下一层循环。

我们考虑 1 到 2 的过程，三点定一个圆，所以一个点是答案的概率为 \dfrac{3}{i}，也就是只有 \dfrac{3}{i} 的概率会进到下一层循环。而 2 到 3 的过程是类似，最后 3 的复杂度是 \mathcal O(j) 的，所以 2 的复杂度是 \sum\limits_{j = 1}^{j \le i} \mathcal O(j) \times \dfrac{3}{i} = \mathcal O(i)，那么可以得到总的时间复杂度是 \sum\limits_{i = 1}^{i \le n} \mathcal O(i) \times \dfrac{3}{i} = \mathcal O(n) 的。

双倍经验：P1742。

code

#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>

using namespace std;

const int N = 1e6 + 5, eps = 1e-12;
struct node{
    double x, y;
}p[N], o;
int n;
double r;

double dist(node u, node v)
{
    return sqrt((u.x-v.x) * (u.x-v.x) + (u.y-v.y) * (u.y-v.y));
}

node get1(node u, node v)
{
    node w;
    w.x = (u.x + v.x) / 2;
    w.y = (u.y + v.y) / 2;
    return w;
}

node get2(node u, node v, node w)
{
    double a, b, c, d, e, f;
    node res;
    a = v.x - u.x, b = v.y - u.y, c = w.x - v.x, d = w.y - v.y;
    e = v.x * v.x + v.y * v.y - u.x * u.x - u.y * u.y;
    f = w.x * w.x + w.y * w.y - v.x * v.x - v.y * v.y;
    res.x = (f * b - e * d) / (c * b - a * d) / 2; 
    res.y = (a * f - e * c) / (a * d - b * c) / 2;
    return res;
}

int main()
{
    scanf("%d", &n);
    for(int i = 1;i <= n;i++)
        scanf("%lf %lf", &p[i].x, &p[i].y);
    random_shuffle(p + 1, p + n + 1);
    o = p[1];
    for(int i = 1;i <= n;i++)
    {
        if(dist(o, p[i]) <= r - eps)
            continue;
        o = p[i], r = 0;
        for(int j = 1;j < i;j++)
        {
            if(dist(p[j], o) <= r - eps)
                continue;
            o = get1(p[i], p[j]);
            r = dist(p[i], p[j]) / 2;
            for(int k = 1;k < j;k++)
            {
                if(dist(p[k], o) <= r)
                    continue;
                o = get2(p[i], p[j], p[k]);
                r = dist(o, p[k]);
            }
        }
    }
    printf("%.2lf %.2lf %.2lf", o.x, o.y, r);
    return 0;
}