最后更新于 2025-08-03 11:00:35

作者

zengyongxu

分类算法·理论

本文为作者原创，转载请注明出处

前置知识点

1 $\pi = \arccos (-1);$

2 数学知识

2.1 余弦定理 $a,b$ 夹角为 $\theta$ 度 $c^2 = a^2 + b^2 - 2ab\cos{\theta}$

2.2 三角函数

2.3 $\color{red}{足够的空间想象能力}$

3 对于一个点的结构体

code

struct point{
    double x, y;
    // 后面一般会重载运算符
};

4 对于一条线的结构体

code

struct line{
    point st, ed;
};

5 关于精度

加法，乘法精度较高
减法，除法，$\cos$，$\sin$，$\arccos$ 等精度较低

浮点数的比较

code

const double eps = 1e-8;

int sign(double x){
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}

int cmp(double x, double y){
    if (fabs(x - y) < eps)
        return 0;
    if (x < y)
        return -1;
    return 1;
}

向量

1 向量的加减法和数乘运算

1.1 $A + B = (A.x + B.x, A.y + B.y)$

code

point add(point a, point b)
{
    return {a.x + b.x, a.y + b.y};
}

1.2 $A - B = (A.x - B.x, A.y - B.y)$

code

point sub(point a, point b)
{
    return {a.x - b.x, a.y - b.y};
}

1.3 $A \times k = (A.x \times k, A.y \times k)$ 其中 $k$ 为常数

code

point mul(point a, double b)
{
    return {a.x * b, a.y * b};
}

2 内积（点积）设 $A,B$ 夹角为 $\theta$，则 $A \cdot B = \mid A\mid \mid B\mid \cos{\theta}$

2.1 几何意义：向量 $A$ 在向量 $B$ 上的投影与B的长度的乘积。

2.2 代码实现

code

double dot(point a, point b)
{
    return a.x * b.x + a.y * b.y;
}

3 外积（叉积） $A \times B = \mid A \mid \mid B \mid sin©$

3.1 几何意义：向量 $A$ 与 $B$ 张成的平行四边形的有向面积。$B$ 在 $A$ 的逆时针方向为正。

3.2 代码实现

code

double cross(point a, point b)
{
    return a.x * b.y - b.x * a.y;
}

4 常用函数

4.1 取模

code

double get_length(point a)
{
    return sqrt(dot(a, a));
}

4.2 计算向量夹角

code

double get_angle(point a, point b)
{
    return acos(dot(a, b) / get_length(a) / get_length(b));
}

4.3 计算两个向量构成的平行四边形有向面积

code

double area(point a, point b, Poinpointt c)
{
    return cross(b - a, c - a);
}

4.4 向量 $A$ 顺时针旋转 $C$ 的角度：

code

point rotate(point a, double angle)
{
    return point(a.x * cos(angle) + a.y * sin(angle), -a.x * sin(angle) + a.y * cos(angle));
}

点与线

1 直线定理

1.1 一般式 $ax + by + c = 0$

1.2 点向式 $p0 + vt$

1.3 斜截式 $y = kx + b$

2 常用操作

2.1 判断点在直线上 $A \times B = 0$

2.2 两直线相交

code

// cross(v, w) == 0则两直线平行或者重合
point get_line_intersection(point p, point v, point q, point w)
{
    point u = p - q;
    double t = cross(w, u) / cross(v, w);
    return p + v * t;
}

2.3 点到直线的距离

code

double distance_to_line(point p, point a, point b)
{
    point v1 = b - a, v2 = p - a;
    return fabs(cross(v1, v2) / get_length(v1));
}

2.4 点到线段的距离

code

double distance_to_segment(point p, point a, point b)
{
    if (a == b)
        return get_length(p - a);
    point v1 = b - a, v2 = p - a, v3 = p - b;
    if (sign(dot(v1, v2)) < 0)
        return get_length(v2);
    if (sign(dot(v1, v3)) > 0)
        return get_length(v3);
    return distance_to_line(p, a, b);
}

2.5 点在直线上的投影

code

point get_line_projection(point p, point a, point b)
{
    point v = b - a;
    return a + v * (dot(v, p - a) / dot(v, v));
}

2.6 点是否在线段上

code

bool on_segment(point p, point a, point b)
{
    return sign(cross(p - a, p - b)) == 0 && sign(dot(p - a, p - b)) <= 0;
}

2.7 判断两线段是否相交

code

bool segment_intersection(point a1, point a2, point b1, point b2)
{
    double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1);
    double c3 = cross(b2 - b1, a2 - b1), c4 = cross(b2 - b1, a1 - b1);
    return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0;
}

多边形

1 三角形

1.1 面积

1.1.1 叉积，略

1.1.2 海伦公式

code

p = (a + b + c) / 2;
S = sqrt(p * (p - a) * (p - b) * (p - c));

1.2 三角形四心

1.2.1 外心，外接圆圆心

三边中垂线交点。到三角形三个顶点的距离相等

1.2.2 内心，内切圆圆心

角平分线交点，到三边距离相等

1.2.3 垂心

三条垂线交点

1.2.4 重心

三条中线交点（到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点）

2 普通多边形

通常按逆时针存储所有点

2.1 定义

2.1.1 多边形

由在同一平面且不再同一直线上的多条线段首尾顺次连接且不相交所组成的图形叫多边形

2.1.2 简单多边形

简单多边形是除相邻边外其它边不相交的多边形

2.1.3 凸多边形

过多边形的任意一边做一条直线，如果其他各个顶点都在这条直线的同侧，则把这个多边形叫做凸多边形

任意凸多边形外角和均为 $360\degree$

任意凸多边形内角和为 $(n−2)180\degree$

2.2 常用函数

2.2.1 求多边形面积（不一定是凸多边形）

我们可以从第一个顶点除法把凸多边形分成n − 2个三角形，然后把面积加起来。

code

double polygon_area(point p[], int n)
{
    double s = 0;
    for (int i = 1; i + 1 < n; i ++ )
        s += cross(p[i] - p[0], p[i + 1] - p[i]);
    return s / 2;
}

2.2.2 判断点是否在多边形内（不一定是凸多边形）

a. 射线法，从该点任意做一条和所有边都不平行的射线。交点个数为偶数，则在多边形外，为奇数，则在多边形内。

b. 转角法

2.2.3 判断点是否在凸多边形内

只需判断点是否在所有边的左边（逆时针存储多边形）。

3 皮克定理

皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为:

$S = a + \frac{b}{2} - 1$

其中 $a$ 表示多边形内部的点数，$b$ 表示多边形边界上的点数，$S$ 表示多边形的面积。

6 圆

6.1 圆与直线交点

6.2 两圆交点

6.3 点到圆的切线

6.4 两圆公切线

6.5 两圆相交面积

凸包

Andrew 算法

1 将点排序：$x$ 为第一关键字从小到大排序，$y$ 为第二关键字从小到大排序

2 从左至右维护上半部分，再从右至左维护下半部分（此处顺序可变，代码中即反过来做的，但是需注意是顺时针方向弹栈）

2.1 使用栈，栈存点，每一个向量都是栈顶两个点组成的向量

2.2 判断栈顶边要不要弹栈。只需要判断当前点在栈顶向量的延长线的逆时针还是顺时针，在逆时针方向弹栈

2.3 注意下半部分算完之后要用第一个点（最左边的点）更新一下。

2.4 代码实现

code

#include <bits/stdc++.h>

using namespace std;

const int N = 1e4 + 5;

int n, st[N];
bool vis[N];

struct point
{
    double x, y;
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
    bool operator<(const point &t)
    {
        if (x == t.x)
            return y < t.y;
        return x < t.x;
    }
} a[N];

double get_dist(point a, point b)
{
    double xx = a.x - b.x;
    double yy = a.y - b.y;
    return sqrt(xx * xx + yy * yy);
}

double cross(point a, point b)
{
    return a.x * b.y - a.y * b.x;
}

double area(point a, point b, point c)
{
    return cross(b - a, c - a);
}

double andrew()
{
    sort(a + 1, a + n + 1);

    int tp = 0;
    for (int i = 1; i <= n; i++)
    {
        while (tp >= 2 && area(a[st[tp - 1]], a[st[tp]], a[i]) <= 0)
        {
            if (area(a[st[tp - 1]], a[st[tp]], a[i]) < 0)
                vis[st[tp--]] = false;
            else
                tp--;
        }
        st[++tp] = i;
        vis[i] = true;
    }

    vis[1] = false;
    for (int i = n; i >= 1; i--)
    {
        if (vis[i])
            continue;
        while (tp >= 2 && area(a[st[tp - 1]], a[st[tp]], a[i]) <= 0)
            tp--;
        st[++tp] = i;
    }

    double res = 0;
    for (int i = 2; i <= tp; i++)
        res += get_dist(a[st[i - 1]], a[st[i]]);
    return res;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> a[i].x >> a[i].y;

    double res = andrew();
    cout << fixed << setprecision(2) << res << "\n";
}

2.5 后记

当然，上面的代码是无法过模版题的，因为需要加特判，判断是否所有的都在一条线上。

此处给出 AC 代码：（注：下代码可能有缺陷，只特判了横着或竖着一条线，但是斜着没有 ~~但是我懒不想写了~~）

code

#include <bits/stdc++.h>

using namespace std;

const int N = 1e5 + 5;

int n, st[N];
bool vis[N];

struct point
{
    double x, y;
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
    bool operator<(const point &t)
    {
        if (x == t.x)
            return y < t.y;
        return x < t.x;
    }
} a[N];

double get_dist(point a, point b)
{
    double xx = a.x - b.x;
    double yy = a.y - b.y;
    return sqrt(xx * xx + yy * yy);
}

double cross(point a, point b)
{
    return a.x * b.y - a.y * b.x;
}

double area(point a, point b, point c)
{
    return cross(b - a, c - a);
}

double andrew()
{
    sort(a + 1, a + n + 1);

    if (a[1].y == a[n].y || a[1].x == a[n].x)
        return 2 * get_dist(a[1], a[n]);

    int tp = 0;
    for (int i = 1; i <= n; i++)
    {
        while (tp >= 2 && area(a[st[tp - 1]], a[st[tp]], a[i]) <= 0)
        {
            if (area(a[st[tp - 1]], a[st[tp]], a[i]) < 0)
                vis[st[tp--]] = false;
            else
                tp--;
        }
        st[++tp] = i;
        vis[i] = true;
    }

    vis[1] = false;
    for (int i = n; i >= 1; i--)
    {
        if (vis[i])
            continue;
        while (tp >= 2 && area(a[st[tp - 1]], a[st[tp]], a[i]) <= 0)
            tp--;
        st[++tp] = i;
    }

    double res = 0;
    for (int i = 2; i <= tp; i++)
        res += get_dist(a[st[i - 1]], a[st[i]]);

    return res;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> a[i].x >> a[i].y;

    double res = andrew();
    cout << fixed << setprecision(2) << res << "\n";
}

此外，还有一种信用卡凸包问题，经分析，其实外围一圈一定是一个圆，至于为什么，毕竟这不是题解。浅谈，多边形外角和 $180 \degree$

code

#include <bits/stdc++.h>

using namespace std;

const int N = 4e4 + 5;
const double pi = acos(-1);
const int dx[] = {1, 1, -1, -1};
const int dy[] = {1, -1, -1, 1};

int n, cnt, tp, st[N];
bool vis[N];

struct point
{
    double x, y;
    bool operator<(const point &t)
    {
        if (x == t.x)
            return y < t.y;
        return x < t.x;
    }
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
} q[N];

struct line
{
    point st, ed;
} ln[N];

point rotate(point a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

double cross(point a, point b)
{
    return a.x * b.y - a.y * b.x;
}

double area(point a, point b, point c)
{
    return cross(b - a, c - a);
}

double get_dist(point a, point b)
{
    double xx = a.x - b.x;
    double yy = a.y - b.y;
    return sqrt(xx * xx + yy * yy);
}

double andrew()
{
    sort(q + 1, q + cnt + 1);
    for (int i = 1; i <= cnt; i++)
    {
        while (tp >= 2 && area(q[st[tp - 1]], q[st[tp]], q[i]) <= 0)
        {
            if (area(q[st[tp - 1]], q[st[tp]], q[i]) < 0)
                vis[st[tp--]] = false;
            else
                tp--;
        }

        st[++tp] = i;
        vis[i] = true;
    }

    vis[1] = false;

    for (int i = cnt; i >= 1; i--)
    {
        if (vis[i])
            continue;
        while (tp >= 2 && area(q[st[tp - 1]], q[st[tp]], q[i]) <= 0)
            tp--;
        st[++tp] = i;
    }

    double res = 0;
    for (int i = 2; i <= tp; i++)
        res += get_dist(q[st[i - 1]], q[st[i]]);
    return res;
}

int main()
{
    double a, b, r;
    cin >> n >> a >> b >> r;

    a = a / 2 - r;
    b = b / 2 - r;

    for (int i = 1; i <= n; i++)
    {
        double x, y, z;
        cin >> x >> y >> z;
        for (int i = 0; i < 4; i++)
        {
            point t = rotate({dx[i] * b, dy[i] * a}, -z);
            q[++cnt] = {x + t.x, y + t.y};
        }
    }

    double res = andrew();
    cout << fixed << setprecision(2) << res + 2 * pi * r << "\n";
}

半平面交

1 将向量按角度排序，使用 atan2(y, x)

2 按照顺序扫描所有向量

2.1 用双端队列维护轮廓

2.2 如果在队头/队尾向量右侧则队头/队尾弹栈

2.3 注意最后用队尾更新队头，队头更新队尾

3.4 使用点向式存储点

3 代码实现

（代码含义为给出 $n$ 个多边形，求每个边向量的半平面交）

code

#include <bits/stdc++.h>

using namespace std;

const int N = 505;
const double eps = 1e-8;

int n, m, cnt;

struct point
{
    double x, y;
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
} pg[N], ans[N];

struct line
{
    point st, ed;
} ln[N];

int q[N];

int sign(double x)
{
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}

int cmp_double(double x, double y)
{
    if (fabs(x - y) < eps)
        return 0;
    if (x < y)
        return -1;
    return 1;
}

double cross(point a, point b)
{
    return a.x * b.y - a.y * b.x;
}

double area(point a, point b, point c)
{
    return cross(b - a, c - a);
}

double get_angle(line a)
{
    return atan2(a.ed.y - a.st.y, a.ed.x - a.st.x);
}

bool cmp(line a, line b)
{
    double A = get_angle(a), B = get_angle(b);
    if (!cmp_double(A, B))
        return area(a.st, a.ed, b.ed) < 0;
    return A < B;
}

point get_line_intersection(point p, point v, point q, point w)
{
    point u = p - q;
    double t = cross(w, u) / cross(v, w);
    return {p.x + v.x * t, p.y + v.y * t};
}

point get_line_intersection(line a, line b)
{
    return get_line_intersection(a.st, a.ed - a.st, b.st, b.ed - b.st);
}

bool on_right(line a, line b, line c)
{
    point o = get_line_intersection(b, c);
    return sign(area(a.st, a.ed, o)) <= 0;
}

double half_plane_intersection()
{
    sort(ln + 1, ln + cnt + 1, cmp);
    int qh = 1, qt = 0;
    for (int i = 1; i <= cnt; i++)
    {
        if (i > 1 && !cmp_double(get_angle(ln[i]), get_angle(ln[i - 1])))
            continue;
        while (qh < qt && on_right(ln[i], ln[q[qt - 1]], ln[q[qt]]))
            qt--;
        while (qh < qt && on_right(ln[i], ln[q[qh]], ln[q[qh + 1]]))
            qh++;
        q[++qt] = i;
    }

    while (qh < qt && on_right(ln[q[qh]], ln[q[qt - 1]], ln[q[qt]]))
        qt--;
    while (qh < qt && on_right(ln[q[qt]], ln[q[qh]], ln[q[qh + 1]]))
        qh++;

    q[++qt] = q[qh];
    int k = 0;
    for (int i = qh; i < qt; i++)
        ans[++k] = get_line_intersection(ln[q[i]], ln[q[i + 1]]);

    double res = 0;
    for (int i = 2; i < k; i++)
        res += area(ans[1], ans[i], ans[i + 1]);
    return res / 2;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> m;
        for (int i = 1; i <= m; i++)
            cin >> pg[i].x >> pg[i].y;
        for (int i = 1; i <= m; i++)
            ln[++cnt] = {pg[i], pg[(i == m ? 1 : i + 1)]};
    }

    double res = half_plane_intersection();
    cout << fixed << setprecision(3) << res << "\n";
}

最小圆覆盖

1 算法用途：

在一个二维平面上给定 $n$ 个点，找出最小的覆盖了所有点的圆。

2 性质（证明去通过这个、这个查询）

2.1 性质一：最小覆盖圆是唯一的。

2.2 性质二：若 $p$ 不在集合 $S$ 的最小覆盖圆内部，则 $p$ 一定在 ${ p } \cup S$ 的最小覆盖圆的边上。

2.3 $\color{red}{三个点确定一个圆（不共线）重要！}$

3 做法

3.1 将整个点集随机化

3.2 现将前一个点的最小覆盖圆设为那一个点

3.2.1 如果 $i$ 不在前 $i - 1$ 个点的最小覆盖圆的内部，则点 $i$ 一定在最小覆盖圆边上。

3.2.2 将圆设为 $i$

3.2.3 第二重枚举：枚举 $1 \leq j \leq i - 1$ 的 $j$

3.2.4 如果 $j$ 不在当前圆边上：则 $j$ 一定在能覆盖 $1 \sim j$ 和 $i$ 且 $i$ 在边上的最小圆的边上。（ $\color{red}{注：由于此处有对于\ i\ 在边上的限制，所以不是\ 1 \sim j\ 和\ i\ 的最小覆盖圆}$ ）

3.2.4.1 第三重枚举：枚举 $1 \leq k \leq j - 1$ 的 $k$

3.2.4.2 如果 $k$ 不在当前圆边上：则 $k$ 一定在能覆盖 $1 \sim k$ 和 $i,j$ 且 $i, j$ 在圆边上的最小圆的边上

3.2.4.3 $\color{red}{发现\ i,j,k\ 已经在圆边上，圆确定！}$ 此时 $break$ 返回上一层

3.2.5 否则就直接 $continue$

3.3 否则就直接 $continue$

3.4 实际上不是 $O(n^3)$ 而是 $O(n)$

4 那给定 3 个点我怎么求圆

4.1 令三个点为 $A,B,C$

4.2 不难发现圆心就是 $AB$ 的中垂线和 $AC$ 的中垂线交点

4.3 中垂线怎么求啊……

4.3.1 以 $AB$ 为例，中垂线的方向就是 $AB$ 向量顺时针旋转 $90 \degree$，起点就是 $\frac{A+B}{2}$ （此处加法是向量加法，除法是数乘，乘$\frac{1}{2}$即可）

5 代码实现

code

#include <bits/stdc++.h>

using namespace std;

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

int n;

struct point
{
    double x, y;
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
    point operator+(const point &t)
    {
        return {x + t.x, y + t.y};
    }
    point operator*(const double &t)
    {
        return {x * t, y * t};
    }
    double operator*(const point &t)
    {
        return x * t.y - y * t.x;
    }
    point operator/(const double &t)
    {
        return {x / t, y / t};
    }
} q[N];

struct circle
{
    point o;
    double r;
};

point get_line_intersection(point p, point v, point q, point w)
{
    auto u = p - q;
    double t = (w * u) / (v * w);
    return p + v * t;
}

point rotate(point a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

pair<point, point> get_line(point a, point b)
{
    return {(a + b) / 2, rotate(b - a, pi / 2)};
}

double get_dist(point a, point b)
{
    double xx = a.x - b.x;
    double yy = a.y - b.y;
    return sqrt(xx * xx + yy * yy);
}

circle get_circle(point a, point b, point c)
{
    pair<point, point> u = get_line(a, b), v = get_line(a, c);
    point p = get_line_intersection(u.first, u.second, v.first, v.second);
    return {p, get_dist(p, a)};
}

int sign(double x)
{
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps)
        return 0;
    if (x < y)
        return -1;
    return 1;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> q[i].x >> q[i].y;

    random_shuffle(q + 1, q + n + 1);
    circle c = {q[1], 0};

    for (int i = 2; i <= n; i++)
    {
        if (dcmp(c.r, get_dist(c.o, q[i])) < 0)
        {
            c = {q[i], 0};
            for (int j = 1; j < i; j++)
            {
                if (dcmp(c.r, get_dist(c.o, q[j])) < 0)
                {
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for (int k = 1; k < j; k++)
                        if (dcmp(c.r, get_dist(c.o, q[k])) < 0)
                            c = get_circle(q[i], q[j], q[k]);
                }
            }
        }
    }

    cout << fixed << setprecision(10) << c.r << "\n"
         << c.o.x << " " << c.o.y << "\n";
}

6 后记

模版题

此外，还有一种变形题是给定椭圆的，具体情况去题目里看。这种题的大致思路是：先旋转，再压缩，就又变成模版了。下面接上 AC 代码：

code

#include <bits/stdc++.h>

using namespace std;

const int N = 50005;
const double eps = 1e-12;
const double pi = acos(-1);

int n;

struct point
{
    double x, y;
    point operator+(const point &t)
    {
        return {x + t.x, y + t.y};
    }
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y};
    }
    point operator*(const double &t)
    {
        return {x * t, y * t};
    }
    point operator/(const double &t)
    {
        return {x / t, y / t};
    }
    double operator*(const point &t)
    {
        return x * t.y - y * t.x;
    }
} q[N];

struct circle
{
    point o;
    double r;
};

int sign(double x)
{
    if (fabs(x) < eps)
        return 0;
    if (x < 0)
        return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps)
        return 0;
    if (x < y)
        return -1;
    return 1;
}

double get_dist(point a, point b)
{
    double xx = a.x - b.x;
    double yy = a.y - b.y;
    return sqrt(xx * xx + yy * yy);
}

point rotate(point a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

point get_line_intersection(point p, point v, point q, point w)
{
    point u = p - q;
    double t = (w * u) / (v * w);
    return p + v * t;
}

pair<point, point> get_line(point a, point b)
{
    return {(a + b) / 2, rotate(b - a, pi / 2)};
}

circle get_circle(point a, point b, point c)
{
    pair<point, point> u = get_line(a, b), v = get_line(a, c);
    point p = get_line_intersection(u.first, u.second, v.first, v.second);
    return {p, get_dist(p, a)};
}

int main()
{
    double a, p;

    cin >> n;
    for (int i = 1; i <= n; i++)
        cin >> q[i].x >> q[i].y;
    cin >> a >> p;

    for (int i = 1; i <= n; i++)
    {
        q[i] = rotate(q[i], a / 180 * pi);
        q[i].x /= p;
    }

    random_shuffle(q + 1, q + n + 1);
    circle c = {q[1], 0};

    for (int i = 2; i <= n; i++)
    {
        if (dcmp(c.r, get_dist(c.o, q[i])) < 0)
        {
            c = {q[i], 0};
            for (int j = 1; j < i; j++)
            {
                if (dcmp(c.r, get_dist(c.o, q[j])) < 0)
                {
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for (int k = 1; k < j; k++)
                        if (dcmp(c.r, get_dist(c.o, q[k])) < 0)
                            c = get_circle(q[i], q[j], q[k]);
                }
            }
        }
    }

    cout << fixed << setprecision(3) << c.r << "\n";
}

$\color{blue}{P.S. 此题精度要求较高}$

三维计算几何

1 三维向量运用 $(x,y,z)$ 表示

code

struct point{
    double x, y, z;
};

2 加法，减法都一样

code

point add(point a, point b)
{
    return {a.x + b.x, a.y + b.y, a.z + b.z};
}

point sub(point a, point b)
{
    return {a.x - b.x, a.y + b.y, a.z + b.z};
}

3 点积（其实也差不多）

3.1 几何意义：令 $A,B$ 夹角为 $\theta$ 度，则 $A \cdot B = \mid A \mid \mid B \mid \cos{\theta}$

3.2 代数求解：令三维向量 $A$ 为 $(x_1, y_1, z_1)$，三维向量 $B$ 为 $(x_2, y_2, z_2)$，则 $A \cdot B = (x_1x_2, y_1y_2, z_1z_2)$

3.3 代码实现

code

double dot(point a, point b){
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

4 模长

求 $A$ 的模长，$\mid A \mid = \sqrt{x^2 + y^2 + z^2}$

code

double len(point a)
{
    return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
}

5 叉积

5.1 几何意义：令 $A,B$ 夹角为 $\theta$ 度，则 $A \cdot B = \mid A \mid \mid B \mid \sin{\theta}$

5.2 代数求解：令三维向量 $A$ 为 $(x_1, y_1, z_1)$，三维向量 $B$ 为 $(x_2, y_2, z_2)$，则 $A \times B = (y_1z_2 - z_1y_2, z_1x_2 - x_1z_2, x_1y_2 - x_2y_1)$

5.3 代码实现

code

point cross(point a, point b)
{
    return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}

5 法向量：任取平面上两个不共线的向量 $A,B$，法向量即 $A \times B$

6 判断一个点 $D$ 是否在平面里

任取平面上两个不共线的向量 $A,B$ ：先求法向量 $C$ = $A\times B$，然后求平面上任意一点到 $D$ 的向量 $E$ 与 $C$ 的点积，判断点积是否为 0

7 多面体欧拉定理：顶点数 - 棱长数 + 表面数 = 2

三维凸包

1 算法使用：

使用增量算法（即每次多考虑一个点）

2 具体细节

2.1 不妨设已有的凸包有 $m$ 个面，现在进行到第 $i$ 个点

2.2 判断当前点在哪些面的上面，则这些面是需要删除的。

2.3 找出一个轮廓线，将这个轮廓线上的每个向量两端向 $i$ 连线，组成一个新的面

2.4 实现过程中，注意用两个数组。

3 代码实现

code

#include <bits/stdc++.h>

using namespace std;

const int N = 205;
const double eps = 1e-10;

int n, m;
bool g[N][N];

double rand_eps()
{
    return ((double)rand() / RAND_MAX - 0.5) * eps;
}

struct point
{
    double x, y, z;
    void shake()
    {
        x += rand_eps(), y += rand_eps(), z += rand_eps();
    }
    point operator-(const point &t)
    {
        return {x - t.x, y - t.y, z - t.z};
    }
    double operator&(const point &t)
    {
        return x * t.x + y * t.y + z * t.z;
    }
    point operator*(const point &t)
    {
        return {y * t.z - z * t.y, z * t.x - x * t.z, x * t.y - y * t.x};
    }
    double len()
    {
        return sqrt(x * x + y * y + z * z);
    }
} q[N];

struct plane
{
    int v[3];
    point norm()
    {
        return (q[v[1]] - q[v[0]]) * (q[v[2]] - q[v[0]]);
    }
    double area()
    {
        return norm().len() / 2;
    }
    bool above(point a)
    {
        return ((a - q[v[0]]) & norm()) >= 0;
    }
} pln[N], np[N];

void get_convex_3d()
{
    pln[++m] = {1, 2, 3};
    pln[++m] = {3, 2, 1};
    for (int i = 4; i <= n; i++)
    {
        int cnt = 0;
        for (int j = 1; j <= m; j++)
        {
            bool t = pln[j].above(q[i]);
            if (!t)
                np[++cnt] = pln[j];
            for (int k = 0; k < 3; k++)
                g[pln[j].v[k]][pln[j].v[(k + 1) % 3]] = t;
        }
        for (int j = 1; j <= m; j++)
        {
            for (int k = 0; k < 3; k++)
            {
                int a = pln[j].v[k], b = pln[j].v[(k + 1) % 3];
                if (g[a][b] && !g[b][a])
                    np[++cnt] = {a, b, i};
            }
        }

        m = cnt;
        for (int j = 1; j <= m; j++)
            pln[j] = np[j];
    }
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; i++)
    {
        cin >> q[i].x >> q[i].y >> q[i].z;
        q[i].shake();
    }

    get_convex_3d();

    double res = 0;
    for (int i = 1; i <= m; i++)
        res += pln[i].area();
    cout << fixed << setprecision(6) << res << "\n";
}

4 小细节

注意到代码中有这样一段：

code

double rand_eps()
{
    return ((double)rand() / RAND_MAX - 0.5) * eps;
}

这样是为了让这些点“颤抖”一下，减小出现 4 点共平面的概率。

5 写在后面

这个题的模版题

warning

计算几何学习记录

前置知识点

浮点数的比较

向量

点与线

多边形

凸包

半平面交

最小圆覆盖

三维计算几何

三维凸包