如果你是一个学算法的人走到“计算几何”这一讲大概率已经有了一种隐约的心理准备前面那些排序、图论、动态规划好歹还能靠“背模板理解状态转移”硬啃下来。但计算几何不一样它把坐标、向量、浮点数全部堆到你面前很多代码在本地跑得好好的交上去就莫名其妙 WA。第九讲选了“凸包”和“最近点对”这两个经典问题其实就是在给你补两样最核心的能力一是用向量判断几何位置关系二是把“看起来只能暴力”的问题用分治或扫描拆到可接受的复杂度。这篇文章没有华丽的背景介绍直接讲原理、代码、边界条件和我在实际练习里踩过的坑适合正在学算法设计、准备机试或者刷竞赛题的人。1. 别急着写凸包先花半小时解决“精度”和“方向判断”很多人学计算几何第一个动作就是打开编辑器抄一份凸包模板结果怎么调都调不通。我个人的建议是先别管凸包是什么先把“判断一个点在一个向量的左侧还是右侧”这件事搞到滚瓜烂熟。因为凸包、最近点对、线段相交、多边形面积几乎所有计算几何问题最终都会落到这类基本谓词上。1.1 叉积是一个“带符号的面积”不是单纯的长度平面里有三个点 A、B、C我们想判断 C 在 AB 的顺时针方向还是逆时针方向。最稳的做法不是算角度而是算叉积struct Point { double x, y; }; // 计算向量 AB 与 AC 的叉积 // 返回 (B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x) double cross(const Point A, const Point B, const Point C) { return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x); }这个值的几何意义是“向量 AB 到向量 AC 的有向面积的两倍”。如果返回值大于 0说明 C 在 AB 的逆时针方向小于 0说明 C 在顺时针方向等于 0说明三点共线。你可能会问为什么不用 atan2 算出角度再比大小因为角度计算涉及反三角函数速度慢而且浮点误差更大。叉积只做几次加减乘在同样的 double 精度下误差累积要小得多而且不需要处理“360° 回绕”这种麻烦问题。这是计算几何所有模板的基石。1.2 EPS所有浮点比较都要留一条“容错带”直接比较两个 double 是否相等是计算几何最常见的翻车现场。比如两个点理论上应该重合但在连续运算后可能分别是 (0.999999999999, 1.000000000001) 和 (1.0, 1.0)。如果你用a.x b.x判断结果就是 false后面整个逻辑全乱。所以标准做法是引入一个极小量 eps通常取 1e-8 到 1e-10const double EPS 1e-8; bool same(const Point a, const Point b) { return fabs(a.x - b.x) EPS fabs(a.y - b.y) EPS; } int sgn(double x) { if (fabs(x) EPS) return 0; return x 0 ? 1 : -1; }sgn函数会把“接近 0”的数统一当作 0 处理。平时判断叉积时不要直接写if (cross 0)而是写if (sgn(cross) 0)。我在初学阶段经常省掉这一步结果构造凸包时遇到三点近乎共线的数据栈顶元素会被错误地弹出最后出来的多边形缺一个角。加了 eps 之后这类“玄学错误”立刻减少了一大半。1.3 排序比较器严格弱序里不能有“相等”无论是凸包算法还是最近点对第一步都是对点集排序。排序的比较器必须满足“严格弱序”也就是说对任意两个点比较器不能同时返回 a b 和 b a。一个典型的反例是直接按欧几里得距离排序或者用fabs(a.x - b.x) EPS当作“相等”后再比 y。这种做法在极少数数据下会导致sort行为未定义甚至死循环。最安全的写法是先按 x 排序如果 x 的差值大于 EPS直接返回 x 较小否则再按 y 排序。bool cmp(const Point a, const Point b) { if (fabs(a.x - b.x) EPS) return a.x b.x; return a.y b.y; }这种比较器看起来朴素但它保证了“唯一确定的全序关系”是后续所有算法稳定的前提。我推荐你在写完比较器之后随手构造一些点跑一遍sort再unique去重确保没有重复点。重复点不处理干净凸包会遇到除零或死循环最近点对会得到距离 0 但位置诡异的答案。2. 凸包用“弹栈”模拟橡皮筋收紧的过程凸包的定义很直观平面上有 n 个点用一个最小的凸多边形把所有点包住这个多边形的顶点就是凸包。可以想象成在钉了很多钉子的木板上拉一根橡皮筋松手后橡皮筋的形状就是凸包。凸包的用途比你想的更广求点集的最小包围范围、判断点是否在多边形内、计算凸多边形的直径旋转卡壳、做简单的碰撞检测甚至在地理信息系统里勾勒一片区域的轮廓。面试和竞赛里最常见的考题是“给定点集求凸包顶点”和“求凸包周长 / 面积”。2.1 三种主流凸包算法怎么选表格里对比的是三种最常见的实现算法时间复杂度编码难度稳定性适用场景Graham 扫描O(n log n)中依赖极角排序共线处理烦教科书教学Andrew 单调链O(n log n)低高按 x/y 排序即可最推荐日常使用Jarvis 步进O(nh)h 是凸包顶点数低高点少或者凸包顶点很少时我自己写题和比赛时基本只用 Andrew 单调链原因很实际Graham 扫描要先求极角再排序遇到多个点在同一极角方向上的时候排序规则要写得很小心否则顺序错了整个凸包就错。而 Andrew 只需要按 x-y 排序然后从左到右扫一遍构造下凸壳从右到左扫一遍构造上凸壳代码短逻辑直接。2.2 Andrew 单调链完整实现算法分三步按 x 升序x 相同按 y 升序排序并去除重复点。从左到右扫描所有点维护一个栈用“新点是否使得栈顶转弯方向不合法”来判断是否弹栈。从右到左再扫一遍构建上半部分最后把两部分拼起来。所谓“转弯方向不合法”就是检查当前栈顶是 p下一个点是 q新点是 r。如果向量 p-q 与 q-r 的叉积小于等于 0或小于 0取决于是否保留共线点说明 r 在 p-q 的顺时针方向或者同向这时候不应该继续留着 q 作为凸包顶点而要弹栈。这里有一个常见分歧保留共线点还是不保留。如果题目要求输出所有在凸包边上的点比较时用 0弹栈如果只要凸包顶点不关心共线点就用 0弹栈。下面的代码采用“不保留共线点”的常用版#include bits/stdc.h using namespace std; const double EPS 1e-8; struct Point { double x, y; bool operator(const Point other) const { if (fabs(x - other.x) EPS) return x other.x; return y other.y; } bool operator(const Point other) const { return fabs(x - other.x) EPS fabs(y - other.y) EPS; } }; double cross(const Point A, const Point B, const Point C) { return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x); } vectorPoint convexHull(vectorPoint pts) { int n pts.size(); if (n 1) return pts; sort(pts.begin(), pts.end()); pts.erase(unique(pts.begin(), pts.end()), pts.end()); n pts.size(); if (n 2) return pts; vectorPoint hull(2 * n); int k 0; // 下凸壳 for (int i 0; i n; i) { while (k 2 cross(hull[k-2], hull[k-1], pts[i]) EPS) { --k; } hull[k] pts[i]; } // 上凸壳 for (int i n - 2, t k 1; i 0; --i) { while (k t cross(hull[k-2], hull[k-1], pts[i]) EPS) { --k; } hull[k] pts[i]; } hull.resize(k - 1); return hull; }注意看第二个循环我从n-2开始因为最右边的点已经在第一个循环里放进去了。t k 1保留了下凸壳的右端点避免在上凸壳构建时把它弹掉。最后resize(k - 1)是因为起点被重复加了一次去掉重复的起点正好把首尾相连的多边形顶点序列凑齐。2.3 一个小例子手动跑一遍假设点集是(0, 0), (1, 1), (2, 2), (1, 0), (0, 2)。按 x-y 排序后是 (0,0), (0,2), (1,0), (1,1), (2,2)。第一个循环压入 (0,0)压入 (0,2)。遇到 (1,0)检查 (0,0)-(0,2)-(1,0) 的叉积这个向量是顺时针方向所以弹掉 (0,2)压入 (1,0)。压入 (1,1) 和 (2,2)。第二个循环从右往左会把中间凹陷的点弹掉最终得到 (0,0), (1,0), (2,2), (0,2)。这就是正确凸包。这个例子虽然小但足以说明“弹栈”的本质栈里保存的是目前已经确定可能在凸包上的点当新点让最近两个点的转弯方向变“错”那中间那个点就不可能是凸包顶点。整个过程和橡皮筋收拢的思路一致。2.4 边界情况不是所有点集都能生成“多边形”很多新手栽在边界上。点数只有 1 或 2凸包退化成点或线段这时候hull.size()返回 1 或 2代码里不能去访问hull[2]。所有点共线时第一个循环会把所有点顺序压进去第二个循环可能因为k t的条件而完全不弹栈最后返回的是一个包含两个端点、中间还夹着一些共线点的大集合。如果你只需要两端点建议在凸包构建后单独判断如果所有点都满足sgn(cross(pts[0], pts[1], pts[i])) 0直接返回{pts.front(), pts.back()}。另外题目如果只要求周长别忘了计算的周长是凸包上相邻点的欧几里得距离之和而hull最后一个点和第一个点也是相邻的。我记得第一次写凸包周长时候末尾点距丢掉了硬是差了一个边的长度最后对拍才查出来。3. 最近点对分治的真正威力在于“剪掉”没用的距离如果说凸包练的是“向量和栈”那最近点对练的就是“分治 合并时如何减少无效比较”。问题很简单给 n 个点找一对欧氏距离最小的点。暴力做法是 O(n²)n 到 1e5 就完全跑不动了。分治做法的复杂度是 O(n log n)代码量并不大但合并时的边界条件是很多人卡住的核心。3.1 分治流程三分治七分合并先把点按 x 坐标排序然后递归处理。假设现在处理的是区间 [l, r] 内的点如果区间里只剩一个点返回无穷大。如果区间里剩两个点直接返回它们距离。取中点 mid递归求左半部分最小距离 d1、右半部分最小距离 d2当前已知最小距离 d min(d1, d2)。合并时只有横坐标与中线距离小于 d 的点才可能刷新答案。把这些点按 y 坐标排序然后对每个点只和它后面“y 差距小于 d”的点比较。很多人会问合并时为什么只比较中线附近 d 距离以内的点因为如果某个点对的距离小于 d但它们的横坐标都距离中线超过 d那么两者横坐标差就超过 d距离必然超过 d不可能成为新的最优解。所以我们可以大胆地把距离中线超过 d 的点全部排除。3.2 为什么内层比较只需要看常数个点这是整个算法最精巧的地方。假设当前最小距离是 d我们选取中线附近的带状区域把区域内所有点按 y 排序。对于其中一个点 p只需要检查它后面的那些“y 坐标差小于 d”的点。为什么不会出现一个点后面有几十个点都要检查因为一旦某两个点的距离小于 d它们不会同时出现在一个 2d × d 的矩形区域内。更具体地说把 d × 2d 的区域划分成若干个 d/2 × d/2 的小方格每个小方格内最多只能有一个点在 2d × d 的范围内点的数量被限制在常数级别。实际编码里很多人直接写“每个点最多检查后面的 7 个点”或者“当 y 坐标差 d 就 break”两种写法都对。我推荐用 y 差 break 的写法不容易出错而且不依赖具体的点数量上界。这部分如果不理解也不影响你写出能 AC 的代码但如果面试官问你“为什么复杂度是 O(n log n) 而不是 O(n²)”你需要能讲清楚“常数点”这个核心。3.3 最近点对完整实现#include bits/stdc.h using namespace std; const double INF 1e20; const double EPS 1e-8; struct Point { double x, y; }; double dist(const Point a, const Point b) { double dx a.x - b.x; double dy a.y - b.y; return sqrt(dx * dx dy * dy); } bool cmpX(const Point a, const Point b) { if (fabs(a.x - b.x) EPS) return a.x b.x; return a.y b.y; } bool cmpY(const Point a, const Point b) { if (fabs(a.y - b.y) EPS) return a.y b.y; return a.x b.x; } double closestPair(vectorPoint pts, int l, int r) { if (r - l 1) { return INF; } int mid (l r) / 2; double midX pts[mid].x; double d min(closestPair(pts, l, mid), closestPair(pts, mid, r)); vectorPoint strip; for (int i l; i r; i) { if (fabs(pts[i].x - midX) d) { strip.push_back(pts[i]); } } sort(strip.begin(), strip.end(), cmpY); int m strip.size(); for (int i 0; i m; i) { for (int j i 1; j m; j) { if (strip[j].y - strip[i].y d) break; d min(d, dist(strip[i], strip[j])); } } return d; } double solve(vectorPoint pts) { sort(pts.begin(), pts.end(), cmpX); return closestPair(pts, 0, pts.size()); }这段代码的关键细节有两个。第一递归写的是closestPair(pts, l, mid)和closestPair(pts, mid, r)半个闭半个开避免mid点被重复计算。第二合并时strip用的是midX这是“当前区间中点的 x 坐标”不是整个点集的中线。如果你误用了全局中线会在递归子区间时把一些实际不在中线附近的点拖进来导致复杂度退化甚至答案错误。3.4 一个容易忽略的性能优化上面的代码每次合并都对新产生的strip做一次sort排序复杂度是 O(k log k)。在分治每一层所有 strip 加起来的长度是 O(n)所以总复杂度是 O(n log² n)。理论上仍然能过绝大多数题但如果是 1e5 以上的点可能有点悬。经典的优化方法是在递归过程中维护一个按 y 排序的辅助数组通过归并的方式让每层合并维持在 O(n)总复杂度严格降到 O(n log n)。不过说实话我大部分时候用上面这版已经能过题了。面试时可以先讲 O(n log² n) 版本再提一句“可以用归并优化到严格 O(n log n)”这反而显得你考虑周到。3.5 手推一个随机小数据点集是 (1, 2), (4, 7), (3, 1), (9, 8), (5, 3), (8, 2)按 x 排序后 (1,2), (3,1), (4,7), (5,3), (8,2), (9,8)mid 取 3左半是 (1,2),(3,1),(4,7)递归后左半最小距离应该是 (1,2) 和 (3,1)距离约 2.236右半是 (5,3),(8,2),(9,8)递归后最小距离约 3.162所以当前 d 2.236。合并时只有横坐标在 4±2.236 范围内的点进 strip。结果里 (4,7) 和 (5,3) 的距离约 4.123(1,2) 和 (3,1) 已经在左半算过不会重复计算。最终答案是 2.236。这个例子虽然简单但能看到合并的本质横向距离过大的点被早早淘汰不需要全部两两比较。4. 从模板到实战凸包和最近点对怎么写才不白练很多人刷题有一个毛病看完模板觉得懂了复制到编辑器里找一道模板题 AC然后就不管了。计算几何要是这样练下次遇到变体题必挂。4.1 一定要学会造随机数据对拍我强烈建议你写一个naive版本暴力 O(n²) 做最近点对或直接暴力枚举所有点对求最小距离然后用脚本生成随机点不断比较暴力结果和分治结果是否一致。对拍是让你“发现边界 bug”的最快方式。比如我曾经在最近点对代码里只对 x 排序、忘了重排 y导致少数随机数据下结果偏大这种错误靠肉眼根本看不出来只有对拍能逮住。下面是一个简单的 Python 对拍脚本思路import random for test in range(1000): n random.randint(1, 20) points [(random.randint(-100, 100), random.randint(-100, 100)) for _ in range(n)] # 用 C 程序跑一份输出和暴力结果对照别追求测试数据大20 个点以内就能暴露全部逻辑问题。等你对拍通过 1000 组再上 1e5 的大数据测性能这样心里才有底。4.2 共线点、重复点、浮点误差是三座大山重复点最近点对的最短距离会变成 0如果没有去重可能导致分治代码把结果算成“两个重复点之间”的距离虽然数值正确但在某些输出要求“输出点对”的问题里会输出一堆相同的点。共线点凸包如果不处理共线会出现凸包顶点上多出一些中间点。如果是求周长或面积问题不大如果是要求输出顶点序列就会 WA。浮点误差别忘了用sgn判断叉积符号。我在一次比赛里因为一个cross 0少写了sgn导致凸包在极大数据下偶尔会漏掉一个点那个 case 恰好卡在 0.0000001 的误差上调了整整一个晚上。4.3 两道题帮你串联这讲的两个算法一道是“Wall”POJ 1113给 n 个点求一个多边形多边形上的点到所有点的距离不小于 L求这个多边形的最小周长。答案其实就是凸包周长加上一个半径为 L 的圆的周长。你看凸包算完还要想到“外扩距离”用圆弧补上而不是简单地周长 L。另一道是“平面最近点对”洛谷 P1429 或 POJ 3714。前者是纯模板后者要求在一个点集和另一个点集之间求最近距离做合并时要注意区分两个集合里的点这是对分治的进一步理解递归时不能把同集合的点对算进去。建议先 AC 这两道再回头读一遍本文的边界情况部分很多当时没感觉的知识点就串起来了。5. 进一步的下一步旋转卡壳、扫描线以及我的学习心得如果第九讲的内容你已经完全拿下了下一步通常会接旋转卡壳Rotating Calipers用它求凸包直径、最大三角形面积、最小宽度等等。旋转卡壳的思想并不难在凸包上拿一条边作为基准让另一个点沿着凸包移动直到距离达到极值然后旋转到下一条边。它和凸包一样也是“维护极值、利用单调性避免重复扫描”的典型。还有一种常见扩展是扫描线类问题比如判断一堆线段是否相交、求矩形面积并。它们的共同点都是“先排序再用事件点/单调结构维护当前扫描线状态”和最近点对里“按 y 排序后只查有限范围”的思路一脉相承。所以我常说计算几何的第一课不是背模板而是理解“怎么减少没必要的比较”。我这里再分享一个个人体会计算几何是最容易靠“画图”来 debug 的算法领域。用 matplotlib 把点、凸包节点、中线和最近点对用不同颜色画出来再标上数字很多 bug 一眼就能看出来。我写的每个凸包模板都会先用随机点画一次图确认形状再上对拍。这比盯着代码空想要高效得多。如果你刚学到第九讲我的建议是先不用追求学完所有算法把叉积、eps、凸包、最近点对这四样彻底吃透配合对拍脚本刷 5 道左右的相关题目计算几何的地基就稳了。之后再去碰旋转卡壳、半平面交你会发现很多“新概念”实际上都是旧思维的变体。
