计算几何学习笔记 [基础 V]

随机增量法

最小圆覆盖

题目：P1742 最小圆覆盖

在平面上给定 \(n\) 个点，求一个最小的圆，能覆盖所有点。要求输出半径和圆心坐标。

这里就需要“随机增量法”。类似数学归纳法。

假设我们已经解决了规模为 \(n-1\) 的问题，那么当问题规模变成 \(n\) 的时候，解决新加入的问题即可。

那么放到这个问题是什么呢？

先任取两个点，求出覆盖它们的最小圆。之后每次增加一个点，若点在圆里，那么不用管它；否则最小圆必过这个新点。为什么呢？

引理：新加入的点若不在圆内，那么最小圆必过新点。

证明就不必要了，可以感性理解。

但是这个办法不优！

新加入的点不在圆内的概率在 \(\dfrac{3}{i}\)，所以期望时间复杂度：

\[\sum_{i = 1}^n O(n^2) \times \dfrac{3}{i} = O(n^2) \]

那么怎么办？我们可以考虑再用增量法优化一次！

算法流程如下：

• 假设已经有一个圆 \(C\) 覆盖前 \(i-1\) 个点。枚举第 \(i\) 个点，那么：
• 若 \(P_i\) 不在圆 \(C\) 内，那么最小圆一定过 \(P_i\)。那么：
  • 更新圆 \(C\)。令 \(P_i\) 为圆心。半径为 \(0\)。
  • 再枚举第 \(j\) 个点。若 \(P_j\) 不在圆内，那么：
    • 更新圆 \(C\)。圆心设为 \(P_i\) 和 \(P_j\) 的中点，半径设为 \(\dfrac{1}{2}dis(P_i,P_j)\)。
    • 枚举第 \(k\) 个点：
      • 若 \(P_k\) 不在圆内，那么 \((P_i,P_j,P_k)\) 三点可以确定一个圆了。更新圆 \(C\) 即可。

伪代码是这样的：

void work()
{
    Circle C;
    for(i = 1 to n)
    {
        if(p[i] 不在圆 C 内)
        {
            C = {p[i], 0};
            for(j = 1 to i - 1)
            {
                if(p[j] 不在圆 C 内)
                {
                    C = {p[i] 和 p[j] 的中点, 0.5 * dis(p[i], p[j])};
                    for(k = 1 to j - 1)
                    {
                        if(k 不在圆 C 内)
                        {
                            C = 过 {p[i], p[j], p[k]} 的圆;
                        }
                    }
                }
            }
        }
    }
}

时间复杂度呢？每层循环都只有不超过 \(\dfrac{3}{i}\) 的概率进入下一层。

\[\begin{aligned} T_1(n) & = \sum_{i = 1}^n T_2(i) \times \dfrac{3}{i} \\ T_2(n) & = \sum_{i = 1}^n T_3(i) \times \dfrac{3}{i} \\ T_3(n) & = O(n) \end{aligned} \]

所以最后 \(T_1(n) = T_2(n) = T_3(n) = O(n)\)。

时间复杂度只有 \(O(n)\) 的优秀算法！

那么我们考虑一些算法细节：如何根据三点求一个圆？

首先，我们可以列个方程组：

\[\begin{cases} (x_1 - x)^2 + (y_1 - y)^2 = r^2 \\ (x_2 - x)^2 + (y_2 - y)^2 = r^2 \\ (x_3 - x)^2 + (y_3 - y)^2 = r^2 \\ \end{cases} \]

之后用 \((1) - (3)\) 和 \((2) - (3)\)。

\[\begin{cases} 2(x_1 - x_2)x + 2(y_1 - y_2)y = (x_1^2 - x_2^2) + (y_1^2 - y_2^2) \\ 2(x_1 - x_3)x + 2(y_1 - y_3)y = (x_1^2 - x_3^2) + (y_1^2 - y_3^2) \end{cases} \]

做个变量代换，写起来方便：

\[\begin{aligned} a_1 & = 2(x_1 - x_2) \\ b_1 & = 2(y_1 - y_2) \\ c_1 & = (x_1^2 - x_2^2) + (y_1^2 - y_2^2) \\ a_2 & = 2(x_1 - x_3) \\ b_2 & = 2(y_1 - y_3) \\ c_2 & = (x_1^2 - x_2^3) + (y_1^2 - y_3^2) \\ \end{aligned} \]

就有了：

\[\begin{cases} a_1 x + b_1 y = c_1 \\ a_2 x + b_2 y = c_2 \end{cases} \]

喜闻乐见的二元一次方程，解一下：

\[\begin{cases} x & = \dfrac{b_1 \times c_2 - b_2 \times c_1}{b_1 \times a_2 - b_2 \times a_1} \\ y & = \dfrac{a_2 \times c_1 - a_1 \times c_2}{a_2 \times b_1 - a_1 \times b_2} \end{cases} \]

解完了！

用计算几何解法也是可以的，因为我们其实很会“已知平面三点，作过三点的圆”的尺规作图做法，可以拿计算几何知识直接模拟。但是，在这里，解析几何也是可以的。

然后为什么是随机增量呢？

先随机打乱所有点，期望复杂度就是对的了。

代码：

#include 
#define DEBUG puts("QAQ")
#define openFile(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout)
#define NOSYNC ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
#define FOR(i,j,k)  for(int (i) = (j); (i) <= (k); ++ (i))
#define RFOR(i,j,k) for(int (i) = (j); (i) >= (k); -- (i))
#define For(i,j,k)  for(int (i) = (j); (i) < (k); ++ (i))
#define RFor(i,j,k) for(int (i) = (j); (i) > (k); -- (i))
#define SC(...) scanf(__VA_ARGS__)
#define PR(...) printf(__VA_ARGS__)
#define N 100005
using namespace std;

typedef double db;
const double eps = 1e-12;
const double PI = acos(-1.0);
int dcmp(db x)
{
    if(fabs(x) < eps) return 0;
    return x > 0 ? 1 : -1;    
}

struct Vector {
    db x, y;
    Vector(db x = 0, db y = 0) : x(x), y(y) {}
    Vector operator + (const Vector &A) const {
        return Vector(x + A.x, y + A.y);
    }
    Vector operator - (const Vector &A) const {
        return Vector(x - A.x, y - A.y);
    }
    Vector operator * (const db &A) const {
        return Vector(x * A, y * A);
    }
    Vector operator / (const db &A) const {
        return Vector(x / A, y / A);
    }
    bool operator == (const Vector &A) const {
        return dcmp(x - A.x) == 0 && dcmp(y - A.y) == 0;
    }
};

typedef Vector Point;

db dis(Point A, Point B)
{
    return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}

struct Circle {
    Point O;
    db r;
};

int n; Point p[N];

Point calc(Point A, Point B, Point C)
{
    db a1 = (A.x - B.x) * 2;
    db b1 = (A.y - B.y) * 2;
    db c1 = (A.x * A.x - B.x * B.x) + (A.y * A.y - B.y * B.y);
    db a2 = (A.x - C.x) * 2;
    db b2 = (A.y - C.y) * 2;
    db c2 = (A.x * A.x - C.x * C.x) + (A.y * A.y - C.y * C.y);
    db x = (b1 * c2 - b2 * c1) / (b1 * a2 - b2 * a1);
    db y = (a2 * c1 - a1 * c2) / (a2 * b1 - a1 * b2);
    return Point(x, y);
}

Circle C;

void work()
{
    C.O = p[1]; C.r = 0;
    FOR(i,2,n)
    {
        if(dcmp(dis(C.O, p[i]) - C.r) > 0)
        {
            C.O = p[i], C.r = 0;
            FOR(j,1,i - 1)
            {
                if(dcmp(dis(C.O, p[j]) - C.r) > 0)
                {
                    C.O = (p[i] + p[j]) / 2, C.r = dis(p[i], p[j]) / 2;
                    FOR(k,1,j - 1)
                    {
                        if(dcmp(dis(C.O, p[k]) - C.r) > 0)
                        {
                            C.O = calc(p[i], p[j], p[k]);
                            C.r = dis(C.O, p[i]);
                        }
                    }
                }
            }
        }
    }
}

int main()
{
	SC("%d", &n);
	FOR(i,1,n) SC("%lf%lf", &p[i].x, &p[i].y);
	random_shuffle(p + 1, p + n + 1);
	work();
	printf("%.10lf\n%.10lf %.10lf\n", C.r, C.O.x, C.O.y);
	return 0;
}