最小圆覆盖

前言

想不出骚话了。

算法讲解

定义

用半径最小的圆覆盖所有的点。

随机增量法

引理1

最小圆覆盖唯一存在。

引理2

给定点集中，最多根据三个点确定最小覆盖圆。

引理3

如果一个点 $P$ 不在集合 $S(P\notin S)$ 的最小圆覆盖内，那么 $P$ 一定在 $P\cup S$ 的最小圆覆盖上。

求解

那么根据以上引理。我们可以得出求解方法：

1. 当前已经求得前 $i-1$ 个点的最小覆盖圆为 $C$。
2. 如果第 $i$ 个点在圆 $C$ 内，那么最小覆盖圆仍然是 $C$。
3. 否则，第 $i$ 个点一定在最小覆盖圆上，接着确定前 $i-1$ 个点还有哪两个（或哪一个）在最小覆盖圆上。把当前的圆心设为 $P_i$，半径为 $0$。
4. 固定了一个点时：找到不在当前覆盖圆上的点，圆心为 $P_i$ 与 $P_j$ 中点，半径为两点距离的一半。
5. 固定了两个点时：找到不在当前覆盖圆上的点，当前覆盖圆为 $P_i$、$P_j$、$P_k$ 的外接圆。

其中求外接圆直接套用向量旋转、向量求交即可。

代码

洛谷最小圆覆盖模板。

Code

#include<iostream>
#include<cstring>
#include<vector>
#include<queue>
#include<string>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<set>
##define put putchar
##define ls(a) (a << 1)
##define rs(a) (a << 1 | 1)
##define pii pair<int, int>

using namespace std;

const int N = 1e5 + 10;
const double eps = 1e-10;

int n;
double r = 0;

struct point {
    double x, y;
    point (const double& x0 = 0, const double& y0 = 0) : x(x0), y(y0) { }
    point operator + (const point& t)  const {
        return point(x + t.x, y + t.y);
    }
    point operator - (const point& t)  const {
        return point(x - t.x, y - t.y);
    }
    point operator * (const double& t)  const {
        return point(x * t, y * t);
    }
    point operator / (const double& t) const {
        return point(x / t, y / t);
    }
    point rotate() {
        return point(y, -x);
    }
} p[N], o;

double dot(const point& a, const point& b) {
    return a.x * b.x + a.y * b.y;
}
double cross(const point& a, const point& b) {
    return a.x * b.y - a.y * b.x;
}
double len(const point& a) {
    return sqrt(dot(a, a));
}
double dis(const point &a, const point &b) {
    return sqrt(dot(a - b, a - b));
}

template<class T>
inline void read(T &x) {
    T res = 0, f = 1; char cc = getchar();
    while(! isdigit(cc)) {
        if(cc == '-') f = -1;
        cc = getchar();
    }
    while(isdigit(cc)) res = (res << 1) + (res << 3) + cc - 48, cc = getchar();
    x = res * f; return ;
}
template<class T, class ...T1>
inline void read(T &x, T1 &...x1) {
    read(x), read(x1...); return ;
}
template<class T>
inline void write(T x) {
    if(x < 0) x = -x, put('-');
    if(x > 9) write(x / 10);
    put(x % 10 + 48); return ;
}
template<>
inline void write(char x) {
    put(x); return ;
}
template<class T, class ...T1>
inline void write(T x, T1 ...x1) {
    write(x),write(x1...); return ;
}
template <class T>
inline bool chkmax(T &x, T y) {
    return x < y ? x = y , 1 : 0;
}

int sign(double a)
{
    if(fabs(a) < eps)
        return 0;
    return a > 0 ? 1 : -1;
}

point get_cross(point p0, point v0, point p1, point v1)
{
    double t = cross(p1 - p0, v1) / cross(v0, v1);
    return p0 + v0 * t;
}

point get_circle(point a, point b, point c)
{
    return get_cross((a + b) / 2, (b - a).rotate(), (a + c) / 2, (c - a).rotate());
}

signed main()
{
   // freopen(".in", "r", stdin);
   // freopen(".out", "w", stdout);
    read(n);
    for(int i = 1; i <= n; ++ i)
        scanf("%lf%lf", &p[i].x, &p[i].y);
    random_shuffle(p + 1, p + n + 1);
    for(int i = 1; i <= n; ++ i)
        if(sign(dis(p[i], o) - r) > 0)
        {
            o = p[i], r = 0;
            for(int j = 1; j < i; ++ j)
                if(sign(dis(p[j], o) - r) > 0)
                {
                    o = (p[i] + p[j]) / 2, r = dis(p[i], p[j]) / 2;
                    for(int k = 1; k < j; ++ k)
                        if(sign(dis(p[k], o) - r) > 0)
                            o = get_circle(p[i], p[j], p[k]), r = dis(o, p[i]);
                }
        }
    printf("%.10lf\n%.10lf %.10lf", r, o.x, o.y);
    return 0;
}

复杂度

这是一个比较 玄学 的复杂度。
点集中最多根据 3 个点即可确定最小覆盖圆，那么每个点参与确定最小点覆盖的概率不大于 $\frac{3}{n}$。
那么进入下一层循环的概率不大于 $\frac{3}{i}$。
则三层循环复杂度分别 $T_1(n)$，$T_2(n)$，$T_3(n)$，那么：
$T_1(n)=O(n)+\sum_{n}^{i=1}\frac{3}{i}T_2(i)$
$T_2(n)=O(n)+\sum_{n}^{i=1}\frac{3}{i}T_3(i)$
$T_3(n)=O(n)$

得 $T_1(n)=T_2(n)=T_3(n)=O(n)$，即复杂度 $O(n)$。

例题

P1742 最小圆覆盖

links，模板题。

P4288 [SHOI2014] 信号增幅仪

links，稍微转化一下即可。
把坐标轴顺时针旋转 $a$ 度，再把 x 轴按比例缩小至 $\frac{1}{p}$，做最小圆覆盖即可。

Code

#include<iostream>
#include<cstring>
#include<vector>
#include<queue>
#include<string>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<set>
##define put putchar
##define ls(a) (a << 1)
##define rs(a) (a << 1 | 1)
##define pii pair<int, int>

using namespace std;

const int N = 1e5 + 10;
const double eps = 1e-6;
const double pi = acos(-1);

int n;
double r = 0.0, the, mul;

struct point {
    double x, y;
    point (const double& x0 = 0.0, const double& y0 = 0.0) : x(x0), y(y0) { }
    point operator + (const point& t)  const {
        return point(x + t.x, y + t.y);
    }
    point operator - (const point& t)  const {
        return point(x - t.x, y - t.y);
    }
    point operator * (const double& t)  const {
        return point(x * t, y * t);
    }
    point operator / (const double& t) const {
        return point(x / t, y / t);
    }
    point rotate(double t) {
        return point(x * cos(t) - y * sin(t), x * sin(t) + y * cos(t));
    }
    point rotate90() {
        return point(-y, x);
    }
} p[N], o;

double dot(const point& a, const point& b) {
    return a.x * b.x + a.y * b.y;
}
double cross(const point& a, const point& b) {
    return a.x * b.y - a.y * b.x;
}
double len(const point& a) {
    return sqrt(dot(a, a));
}
double dis(const point &a, const point &b) {
    return sqrt(dot(a - b, a - b));
}

template<class T>
inline void read(T &x) {
    T res = 0, f = 1; char cc = getchar();
    while(! isdigit(cc)) {
        if(cc == '-') f = -1;
        cc = getchar();
    }
    while(isdigit(cc)) res = (res << 1) + (res << 3) + cc - 48, cc = getchar();
    x = res * f; return ;
}
template<class T, class ...T1>
inline void read(T &x, T1 &...x1) {
    read(x), read(x1...); return ;
}
template<class T>
inline void write(T x) {
    if(x < 0) x = -x, put('-');
    if(x > 9) write(x / 10);
    put(x % 10 + 48); return ;
}
template<>
inline void write(char x) {
    put(x); return ;
}
template<class T, class ...T1>
inline void write(T x, T1 ...x1) {
    write(x),write(x1...); return ;
}
template <class T>
inline bool chkmax(T &x, T y) {
    return x < y ? x = y , 1 : 0;
}

int sign(double a)
{
    if(fabs(a) < eps)
        return 0;
    return a > 0 ? 1 : -1;
}

point get_cross(point p0, point v0, point p1, point v1)
{
    double t = cross(p1 - p0, v1) / cross(v0, v1);
    return p0 + v0 * t;
}

point get_circle(point a, point b, point c)
{
    return get_cross((a + b) * 0.5, (b - a).rotate90(), (a + c) * 0.5, (c - a).rotate90());
}

signed main()
{
   // freopen(".in", "r", stdin);
   // freopen(".out", "w", stdout);
    read(n);
    for(int i = 1; i <= n; ++ i)
        scanf("%lf%lf", &p[i].x, &p[i].y);
    scanf("%lf%lf", &the, &mul), the = -the / 180.0 * pi, mul = 1.0 / mul;
    for(int i = 1; i <= n; ++ i)
        p[i] = p[i].rotate(the), p[i].x *= mul;
    random_shuffle(p + 1, p + n + 1);
    for(int i = 1; i <= n; ++ i)
        if(sign(dis(p[i], o) - r) > 0)
        {
            o = p[i], r = 0.0;
            for(int j = 1; j < i; ++ j)
                if(sign(dis(p[j], o) - r) > 0)
                {
                    o = (p[i] + p[j]) * 0.5, r = dis(p[i], p[j]) * 0.5;
                    for(int k = 1; k < j; ++ k)
                        if(sign(dis(p[k], o) - r) > 0)
                            o = get_circle(p[i], p[j], p[k]), r = dis(o, p[i]);
                }
        }
    printf("%.3lf", r);
    return 0;
}

尾声

Just this.