题解 P1742 【最小圆覆盖】

· · 题解

最小圆覆盖 luoguP1742

给定n个点，求一个最小的圆包围所有的点。

随机增量法

定理1：如果点p不在集合S的最小覆盖圆内，则p一定在S\cup\{p\}的最小覆盖圆上。

根据这个定理，我们可以分三次确定前i个点的最小覆盖圆。

写成伪代码的形式非常简单：

圆 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 = {0.5*(P[i]+P[j]), 0.5*dist(P[i], P[j])};
                for(k=1 to j-1)
                {
                    if(P[k] 不在 C 内)
                    {
                        C = 外接圆(P[i], P[j], P[k]);
                    }
                }
            }
        }
    }
}

只需要三个模式完全相同的for循环就可以搞定。

但是对于这个算法，还有三个地方需要明确：如何求外接圆，以及复杂度是多少。

外接圆

可以使用初中的中垂线法。设三个不共线点为A,B,C，那么我们可以求AB,AC中垂线的交点。

以下我们需要四个工具去完成这个任务：

第一个任务，求两个向量的中点，将横纵坐标相加除以二即可。

第二个任务，将一个向量旋转90度，如果是顺时针旋转，即(x,y)变成(y,-x)

第三个任务，已经说白了，两个向量就可以表示一条直线。

第四个任务稍微有一点复杂，但是不是很复杂。众所周知，二维向量的叉积可以表示两向量所形成的平行四边形的有向面积。我们可以利用这一点解决直线交点的问题。

用p_1,v_1表示I，即：I=p_1+tv_1

又I在绿色直线上，所以又向量(I-p_2)和v_2构成的平行四边形面积为0，即(p_1+tv_1-p_2)\times v_2 = 0

解方程得，t=\frac{(p_2-p_1)\times v_2}{v_1\times v_2}，交点I就是p_1+tv_1

利用上述四个工具，我们可以求出三个点的外接圆。

复杂度证明

由于一堆点最多只有3个点确定了最小覆盖圆，因此n个点中每个点参与确定最小覆盖圆的概率不大于\frac{3}{n}

所以，每一层循环在第i个点处调用下一层的概率不大于\frac{3}{i}

那么设算法的三个循环的复杂度分别为T_1(n),T_2(n),T_3(n)，则有：

\begin{aligned}T_1(n) & = O(n) + \sum_{i=1}^{n}{\frac{3}{i}T_2(i)}\\T_2(n) & = O(n) + \sum_{i=1}^{n}{\frac{3}{i}T_3(i)}\\T_3(n) & = O(n)\end{aligned}

不难解得，T_1(n)=T_2(n)=T_3(n)=O(n)

代码

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

using namespace std;

struct vec
{
    double x, y;
    vec (const double& x0 = 0, const double& y0 = 0) : x(x0), y(y0) {}
    vec operator + (const vec& t) const {return vec(x+t.x, y+t.y);}
    vec operator - (const vec& t) const {return vec(x-t.x, y-t.y);}
    vec operator * (const double& t) const {return vec(x*t, y*t);}
    vec operator / (const double& t) const {return vec(x/t, y/t);}
    const double len2 () const {return x*x + y*y;}
    const double len () const {return sqrt(len2());}
    vec norm() const {return *this/len();}
    vec rotate_90_c () {return vec(y, -x);}
};

double dot(const vec& a, const vec& b) {return a.x*b.x + a.y*b.y;}
double crs(const vec& a, const vec& b) {return a.x*b.y - a.y*b.x;}

vec lin_lin_int(const vec& p0, const vec& v0, const vec& p1, const vec& v1)
{
    double t = crs(p1-p0, v1) / crs(v0, v1);
    return p0 + v0 * t;
}

vec circle(const vec& a, const vec& b, const vec& c)
{
    return lin_lin_int((a+b)/2, (b-a).rotate_90_c(), (a+c)/2, (c-a).rotate_90_c());
}

int n;
vec pot[100005];

int main()
{
    scanf("%d", &n);
    for(int i=1; i<=n; i++) scanf("%lf%lf", &pot[i].x, &pot[i].y);
    random_shuffle(pot+1, pot+n+1);
    vec o;
    double r2 = 0;
    for(int i=1; i<=n; i++)
    {
        if((pot[i]-o).len2() > r2)
        {
            o = pot[i], r2 = 0;
            for(int j=1; j<i; j++)
            {
                if((pot[j]-o).len2() > r2)
                {
                    o = (pot[i]+pot[j])/2, r2 = (pot[j]-o).len2();
                    for(int k=1; k<j; k++)
                    {
                        if((pot[k]-o).len2() > r2)
                        {
                            o = circle(pot[i], pot[j], pot[k]), r2 = (pot[k]-o).len2();
                        }
                    }
                }
            }
        }
    }
    printf("%.10lf\n%.10lf %.10lf\n", sqrt(r2), o.x, o.y);
    return 0;
}