实战详解)
算法设计与分析这门课走到分治那一章实验二基本都会被安排成最近点对问题。题目描述通常只有两三行平面上给定 n 个点求距离最近的一对点以及它们的距离。我第一次拿到这个题的时候心里想的是这玩意儿有什么好做的两重循环十几行就写完了直到看见数据规模写着 n ≤ 10⁵运行时间限制 1 到 2 秒才反应过来自己把问题想简单了。暴力枚举在 n 1000 的时候还能秒出结果到了 10 万这个量级运算次数是 5×10⁹ 级别等它跑完饭都凉了。分治法求解最近点对是分治思想里少数几个原理一看就懂、代码一写就错的经典案例。它不像归并排序那样把递归写对就万事大吉中间那一步跨中线的合并藏着好几个能把人折磨半天的细节。下面我按自己当年做这个实验的顺序把思路、代码、踩过的坑和实验报告怎么写一次说清楚。适合正在做这个实验的同学也适合想把分治真正用起来、而不是只会背递推式的人。1. 暴力枚举撞墙之后分治的切入点在哪1.1 先算一笔账看看 O(n²) 究竟能撑到多大很多人对复杂度的直觉是模糊的觉得平方慢一点但也能跑。我们把账算具体一点一次距离计算包含两次减法、两次乘法、一次加法和一次开方在普通 PC 上大约几十纳秒。n 10⁴ 时点对数量约 5×10⁷勉强能在半秒到一秒之间跑完n 10⁵ 时点对数量约 5×10⁹按每次 20 纳秒算就是 100 秒而实际上因为内存访问模式和开方指令的延迟真实耗时往往还要更久。数据规模 n暴力枚举点对数实测大致耗时1,000约 5×10⁵10 毫秒以内10,000约 5×10⁷0.5 到 1 秒100,000约 5×10⁹100 秒以上1,000,000约 5×10¹¹彻底别想这张表里的时间随机器差异会有波动但量级关系是确定的暴力法的时间随规模是平方增长规模翻十倍时间翻一百倍。分治法把复杂度压到 O(n log n)n 从 10⁴ 涨到 10⁵时间只涨大概十倍多一点点。10 万点的分治实现实测通常在几十毫秒以内和暴力法不是一个世界的东西。1.2 一分为二之后最近点对只有三种可能分治的套路就是分、治、合。把点集按 x 坐标排序从中间切开左边一半右边一半各自递归求出左边内部的最近距离 dₗ 和右边内部的最近距离 dᵣ取 d min(dₗ, dᵣ)。到这里为止都很顺。关键问题在合最近的那一对点有没有可能一个在左边、一个在右边当然有可能。所以整张图上的最近点对只可能是三种情况之一完全落在左半边、完全落在右半边、横跨中线。前两种已经被递归处理掉了我们唯一需要额外操心的就是第三种。这时候一个非常自然但错误的简化冒出来了有人说既然点已经按 x 排好序了那最近点对肯定出现在 x 相邻的两个点之间扫一遍相邻点不就行了这是这个实验里最经典的思维陷阱。反例随手就能构造三个点是 (0, 0)、(0, 100)、(1, 0)。按 x 排序后是 (0,0)、(0,100)、(1,0)x 相邻的两对距离分别是 100 和约 100.005但真正的最近点对是 (0,0) 和 (1,0)距离 1。x 相邻完全不意味着距离近因为 y 方向可能差得很远。按 y 排序后比较相邻点同样错。反例是把上面的例子转 90 度。所以任何排一次序扫一遍的简化都是不成立的必须老老实实做合并。1.3 跨中线的那一步才是整个算法的灵魂真正让 O(n log n) 成立的是合并阶段的一个剪枝观察如果一个点对横跨中线并且它的距离小于当前的 d min(dₗ, dᵣ)那么这两个点各自到中线的水平距离一定都小于 d。原因很简单点对距离小于 d那么它在 x 方向和 y 方向的投影分量也都小于 d。于是我们只需要考察一条竖直的带子以中线为中心左右各宽 d总宽 2d。只有落在这条带子里的点才有可能构成优于 d 的跨中线点对。带子外面的点它到中线的距离已经不小于 d 了和对面任何一个点的距离必然不小于 d直接排除。这一步把候选点从 n 个压缩到了期望 O(1) 个对随机分布的数据而言。剩下的问题就是在这条带子里怎么高效地找最近点对这就引出了下一节的核心技巧。2. 递归骨架里的三处设计每一处都能写错2.1 排序只做一次之后靠归并维护有序性最直觉的实现是每一层递归都把带子里的点拿出来按 y 排一次序然后暴力比较。这样做是对的复杂度是 O(n log² n)10 万点跑起来大概几百毫秒交作业其实够用了很多人就是这么写的。但如果你想拿到满分的效率正确做法是把按 y 排序这件事也交给归并来做。思路是递归函数返回时保证它管辖的这段区间内的点已经按 y 有序。那么父层拿到左右两个已经按 y 有序的子区间后只需要做一次归并O(n)就能得到整个区间按 y 有序的序列。这样一来每层的额外开销从排序的 O(n log n)降到了归并的 O(n)总复杂度就是标准的 O(n log n)。这个技巧的本质是排序和分治的递归结构长得一模一样那就没必要排两次。归并排序本身就是分治最近点对也是分治把两者的递归树合并成一棵就能省掉一个 log。这是分治优化里非常通用的思路后面处理别的题目也能反复用到。2.2 中线的位置要在递归之前取出来这是我见过最多人翻车的地方而且翻车之后现象特别诡异小数据全对大数据偶尔错加个打印一看结果又对了。原因就出在中线取在哪。如果你用同一个数组同时承载按 x 排序和按 y 排序两种状态那么在递归返回之后数组的物理顺序已经被按 y 重排过了。这时候你再去取p[mid].x作为中线坐标取到的根本不是原来那个 x 中位数而是一个被 y 序打乱之后的随机点的坐标。中线偏了带子就歪了被剪掉的点里可能就藏着最优解。正确的写法是进入函数后、调用递归之前先把midx p[mid].x存下来。这个值在整个函数体里都不会变递归怎么折腾数组都不影响它。我在注释里一般会特意写一句必须在递归前取因为过两个月回头看代码自己都会怀疑这里为什么要这么写。2.3 带子的宽度用严格小于边界上的点可以不管收集带子里的点时判断条件是fabs(p[i].x - midx) d。用严格小于而不是小于等于是安全的。因为距离恰好等于 d 的点对对我们求最小值没有贡献剪掉它不会让答案变差。同理内层循环里的剪枝条件写dy d就 break因为 y 方向的距离已经不小于 d 了而完整距离一定不小于任何单个方向上的分量所以这对点以及后面 y 更大的点都不可能刷新答案直接跳出内层循环。这两个不等号的方向如果写反会带来两类不同的后果。带子条件写成 d只是稍微多算了几个点答案不会错属于性能小损但内层剪枝如果写成dy d并且用了非严格的排序严格相等的边界情况会多算几个点同样不影响正确性。真正致命的是把带子写成了fabs(...) d/2这种想当然的缩小那就直接漏解了。2.4 为什么每个点只需要和常数个邻居比较带子里按 y 升序排列之后对每个点我们只需要往后看常数个点就够了不需要看完整条带子。这个常数从哪来给一个能在实验报告里写清楚的论证。假设当前最近距离是 d。拿出以中线为中心、宽 2d左右各 d、高 d 的一个矩形。对这个矩形里的任意一个点 p把它放在矩形的下边界上我们考察这个矩形内、y 坐标不小于 p 的所有点。把这 2d × d 的矩形沿中线劈成左右两个 d × d 的正方形每个正方形再切成四个边长 d/2 的小正方形。关键观察来了每个 d/2 边长的小正方形里最多只能有一个点。因为同一个小正方形内任意两点的距离最大是它的对角线也就是 (√2/2)d ≈ 0.707d小于 d。而 d 是左半边或右半边内部的最近距离同一侧出现两个点距离小于 d矛盾。所以每个 d × d 正方形最多 4 个点整个 2d × d 矩形最多 8 个点。除去 p 自己每个点最多只需要和 7 个点比较。这就是为什么内层循环可以放心地在有限的常数次之后 break也正是 O(n log n) 里那个 O(n) 由来的依据。在实际代码里有人直接硬编码往后看 7 个点有人用dy d的剪枝动态判断。我更推荐后者因为它不依赖那个常数的具体数值逻辑自洽改数据结构时也不容易出错。3. 一份能直接跑通的实现逐段拆开讲3.1 点的结构体和比较函数#include bits/stdc.h using namespace std; struct Point { double x, y; int id; // 原始编号便于输出是第几个点 }; vectorPoint p; // 全局数组递归过程中反复重排 double bestDist 1e300; int bestI -1, bestJ -1; inline 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 (a.x ! b.x) return a.x b.x; return a.y b.y; // 次级关键字别省 } bool cmpy(const Point a, const Point b) { if (a.y ! b.y) return a.y b.y; return a.x b.x; }两个比较函数都写了次级关键字cmpx在 x 相等时比 ycmpy在 y 相等时比 x。这不是可有可无的美化而是为了保证排序结果在任何平台上都完全一致。只写主关键字的话标准库的排序是不保证相同元素相对顺序的同一份数据在不同编译器或不同运行次数下可能排出不同的顺序对于带重复点的测试用例会让哪一对点被输出变得不可预测。如果题目对相等距离时的输出有编号要求这个不确定性会直接让你 WA。3.2 递归函数的主体// 处理 p[l..r]返回时保证 p[l..r] 按 y 升序 void solve(int l, int r) { int n r - l 1; if (n 3) { // 递归基小规模直接暴力同时按 y 排好 for (int i l; i r; i) for (int j i 1; j r; j) { double d dist(p[i], p[j]); if (d bestDist) { bestDist d; bestI min(p[i].id, p[j].id); bestJ max(p[i].id, p[j].id); } } sort(p.begin() l, p.begin() r 1, cmpy); return; } int mid (l r) 1; double midx p[mid].x; // 必须在这里取递归之后就不是它了 solve(l, mid); solve(mid 1, r); // 左右两段各自已按 y 有序一次归并搞定本层 inplace_merge(p.begin() l, p.begin() mid 1, p.begin() r 1, cmpy); double d bestDist; // 取当前全局最优用于带子剪枝 vectorint strip; strip.reserve(n); for (int i l; i r; i) if (fabs(p[i].x - midx) d) strip.push_back(i); for (size_t i 0; i strip.size(); i) { for (size_t j i 1; j strip.size(); j) { double dy p[strip[j]].y - p[strip[i]].y; if (dy d) break; // y 方向已拉开后面只会更远 double cur dist(p[strip[i]], p[strip[j]]); if (cur bestDist) { bestDist cur; bestI min(p[strip[i]].id, p[strip[j]].id); bestJ max(p[strip[i]].id, p[strip[j]].id); } } } }递归基选 n ≤ 3 有两个考虑一是小规模暴力比较的常数极小比继续递归划算二是必须在这个出口把这一段按 y 排好否则父层的inplace_merge前提不成立。我见过有人忘了这个 sort结果归并出来的顺序是乱的带子里的 y 序剪枝失效答案就飘了。double d bestDist;这一行放在两次递归之后是必须的。递归过程中 bestDist 可能已经被更新得更小用最新的值做带子宽度能让带子更窄、剪枝更狠。如果习惯用返回值传递子问题的最近距离那就写double d min(dl, dr);效果一样。3.3 主函数与数据读入int main() { int n; scanf(%d, n); p.resize(n); for (int i 0; i n; i) { scanf(%lf %lf, p[i].x, p[i].y); p[i].id i; // 0-based输出时记得 1 } sort(p.begin(), p.end(), cmpx); // 进入递归前必须按 x 有序 solve(0, n - 1); printf(最近距离: %.6f\n, bestDist); printf(最近点对: 第 %d 个点 和 第 %d 个点\n, bestI 1, bestJ 1); return 0; }sort(p.begin(), p.end(), cmpx)这一行是整个算法能够成立的前提绝对不能漏。递归里的mid (l r) 1只有在数组按 x 有序时才真的把点集沿 x 方向切成了两半。如果忘了预排序那段代码切出来的就是随机两半结果错得毫无规律而且小数据量下因为递归基是暴力还可能歪打正着给出正确答案非常具有迷惑性。编号的输出要留意题目要求。有些题目从 1 开始编号有些从 0 开始还有些要求输出坐标而不是编号。这个小细节每年都能送走一批人。3.4 不想动脑时的保底写法如果时间紧或者对inplace_merge的边界没把握可以用一个更笨但更不容易错的版本递归函数不再维护 y 序而是在每次合并时把带子里的点单独取出来复制到一个临时数组用sort按 y 排一次然后暴力扫相邻的若干个。复杂度退化到 O(n log² n)但 10 万点的量级下实测通常也就在几百毫秒这个区间大部分评测机是能过的。这种写法的好处是把数组重排和归并维护有序这两件事解耦了中线取值的时机问题自然消失因为原数组的顺序从头到尾没被动过。代价是常数大一些逻辑上多了一次临时数组的拷贝。4. 从样例通过到全面正确中间隔着这些坑4.1 精度问题的处理位置距离计算里有一次开方理论上会引入浮点误差。但要注意这个误差并不会影响正确性判断因为开方是单调运算如果真实距离 a b那么算出来的 sqrt(a²) 和 sqrt(b²) 之间的大小关系仍然一致除非两者差距小于浮点精度。坐标是普通范围内的整数时这种情况不会出现。真正需要留神的是两件事。一是不要自作聪明地加 eps。有人写成if (cur bestDist - 1e-9)觉得这样可以避免浮点抖动实际上这会让一些本来应该被记录的点对被忽略如果题目要求输出具体的点对编号就会输出一个不是最优的答案。二是输出格式距离通常要求保留 6 位小数用printf(%.6f)就对别用cout的默认精度默认只保留 6 位有效数字坐标大一点比如 1e5 量级时小数部分会被全砍掉。如果坐标范围很大比如到 10⁹x 和 y 都是整数那么 dx * dx 会到 10¹⁸ 量级用 int 或者 float 都会溢出或丢精度必须用 double 或者 long long 来存中间结果。这一点在三维扩展或者曼哈顿距离的变体里尤其容易出问题。4.2 递归基写成 n ≤ 2 会怎样把递归基设成 n ≤ 2 也能跑对但会让递归层数多一层而且当 n 3 时继续切会切出 1 和 2 两段白白增加函数调用开销。设成 n ≤ 3 是常见的折中。设得更大比如 n ≤ 8 也可以小规模暴力的常数比递归小能省一点时间但要注意递归基越小暴力那部分的比较次数就越少收益递减一般 3 到 5 之间就够了。需要强调的是无论递归基设成几那个返回前按 y 排好序的动作都不能少。如果用了保底的每层单独排序写法那就无所谓因为不存在跨层的顺序约定。4.3 重复点和共线点输入里有重复点的时候最近距离是 0这个结果本身没问题。但如果题目允许重复点那么带子剪枝里dy d这个条件在 d 已经等于 0 的情况下会立刻 break内层循环几乎不执行。所以对重复点的情况算法给出的 0 是靠递归基里的暴力比较或者更早的某一层拿到并更新到 bestDist 的后面的剪枝只是不再做事而已。如果 d 一开始就是 0还有一种可能出问题fabs(p[i].x - midx) d这个条件在 d 0 时会筛不出任何点带子为空。这没问题因为答案已经找到了。但如果你在代码里除了求最小值还要求输出距离最小的点对中编号最小的那一对那么在 d 0 且存在多对重复点的时候需要在相等时比较编号而不是简单地用。题目不要求的话就别加加了容易引入新的 bug。共线的点所有点 x 坐标相同会触发另一个现象。此时按 x 排序后递归的切分点完全靠索引 midpoint 决定左右两边的点 x 坐标全都相等带子宽度判断fabs(x - midx) d会把所有点都收进带子里剪枝退化成暴力。这种情况的复杂度会从 O(n log n) 退化到 O(n²)但一般不会有人拿这种极端数据来卡你知道有这回事就行。4.4 对拍让随机数据替你找 bug手写五六个用例是发现不了逻辑漏洞的必须上对拍。做法是写一个绝对正确的暴力程序再写一个数据生成器和一个循环脚本跑上几百组随机数据比对输出。import random, subprocess def brute(pts): best float(inf) bi bj -1 for i in range(len(pts)): for j in range(i 1, len(pts)): d ((pts[i][0] - pts[j][0]) ** 2 (pts[i][1] - pts[j][1]) ** 2) ** 0.5 if d best: best, bi, bj d, i, j return best, bi, bj for t in range(500): n random.randint(2, 40) R random.choice([5, 50, 1000]) pts [(random.randint(-R, R), random.randint(-R, R)) for _ in range(n)] data f{n}\n \n.join(f{x} {y} for x, y in pts) \n out subprocess.run([./solve], inputdata, capture_outputTrue, textTrue).stdout got float(out.strip().split()[0]) if out.strip() else None exp, _, _ brute(pts) if got is None or abs(got - exp) 1e-6: print(不一致的错误用例:) print(data) print(期望, exp, 实际, got) break else: print(500 组随机数据全部通过)这里有两个设计上的小心思。一是坐标系范围要变换着来用小范围比如 ±5是为了制造大量重复点和共线情况用大范围是为了覆盖一般情形。二是随机数据的规模要控制在几百以内这样才能用暴力程序算期望值如果一上来就是 10 万点对拍就无从谈起了。跑通对拍之后再把规模放大到 10⁵ 做性能测试这样正确性和性能都有了保障。5. 计时实验和复杂度论证报告里怎么写才站得住5.1 把递推式展开一遍分治的复杂度分析在报告里必须写清楚而且不能只写一句由主定理得 T(n) O(n log n)那样太敷衍也体现不出你真的理解了。推导过程是这样的每层递归把问题分成两个规模为 n/2 的子问题加上合并阶段的 O(n) 开销得到 T(n) 2T(n/2) O(n)递归基是常数。把递归树一层层展开第一层代价是 cn第二层是 2 × c(n/2) cn第三层是 4 × c(n/4) cn每一层的总代价都是 cn。递归树一共有多少层每次规模减半从 n 到 1 需要 log₂n 层。所以总代价是 cn × log₂n即 O(n log n)。这个展开过程比直接套主定理有说服力得多因为它顺带说明了每一层的合并为什么必须是 O(n)。如果合并阶段退化成了 O(n log n)比如每层重新排序那么总复杂度就变成 O(n log² n)正好对应前面提到的那种保底写法。5.2 计时实验怎么设计计时部分的设计比很多人想的讲究。几个要点数据规模至少要取 6 到 8 个点比如 10³、5×10³、10⁴、5×10⁴、10⁵、5×10⁵这样才能看出增长趋势。只测两三个点画出来就是一条直线段说明不了问题。每个规模至少跑 5 次去掉最大值后取平均。单次计时受系统调度、缓存状态影响很大尤其是毫秒级的程序波动可能有几倍。取多次平均是最基本的严谨性。生成数据时坐标范围要合理。如果坐标范围取 0 到 n那么点数越多点越密最近距离越小带子里的点数也越少测出来的时间可能比理论预期还漂亮。如果想让实验更有代表性可以让坐标范围固定比如 0 到 10⁶这样点密度随规模变化更接近真实场景。计时函数用chrono::high_resolution_clock单位换算成毫秒输出时保留三位小数。别用clock()它统计的是 CPU 时间在多核或者有系统调用的情况下容易给出奇怪的结果。把实测数据和理论曲线画在同一张图上纵轴可以用对数坐标如果两条线基本平行就说明实测符合 O(n log n)这是报告里最有说服力的一页。5.3 报告里容易被扣分的写法几个我见过的高频问题。一是只贴代码不解释关键变量的作用尤其是midx提前取、递归基返回前排序这两处不写清楚就等于把最难的部分藏起来了。二是复杂度分析里把带子里最多 8 个点当成理所当然不做论证评阅人一看就知道是抄的。三是性能对比只说分治比暴力快很多不给数据、不给规模、不给机器环境这种结论没有意义。还有一个隐蔽的问题不要在报告里写该算法稳定这类含糊的表述容易引起歧义。想说排序结果可复现就直接说比较函数里带了次级关键字保证相同主关键字的元素顺序确定。6. 这套分治框架还能往哪搬6.1 三维最近点对为什么不能照搬把问题提到三维点的距离变成三个方向分量的平方和开方分治的骨架还是按 x 切、递归、考察跨中线但合并阶段会变得麻烦很多。二维时我们靠一个 2d × d 的矩形论证出里面最多 8 个点三维时要考察的是一个 2d × d × d 的立方体里面的点数上限会涨到 48 个左右常数大了不少代码里的剪枝逻辑也不再是简单的y 方向拉开就 break还要考虑 z 方向。实际工程里遇到高维最近点对一般不会硬上分治而是用 KD 树配合剪枝或者直接用随机化的近似算法。分治法在高维上的常数代价太大了理论复杂度好看但实际跑不过。6.2 换了距离度量会怎样如果把欧氏距离换成曼哈顿距离分治框架可以直接用带子剪枝的逻辑也不变因为曼哈顿距离同样满足分量的和小于 d 蕴含每个分量都小于 d这个性质。换成切比雪夫距离各分量差的最大值也一样。但如果你要做的是最大距离点对也就是求距离最远的一对点分治就完全不适用了。那个问题有更漂亮的解法最大距离点对一定在凸包的顶点上先求凸包再旋转卡壳复杂度 O(n log n)。同一个点集求最小距离用分治求最大距离用凸包这个对比本身就很能说明算法的选择取决于问题结构而不是套模板。6.3 分治里跨中线这个套路的迁移跨中线这一步的本质是把全局最优解按照是否跨越切分线分成两类递归负责不跨越的那一类合并阶段负责跨越的那一类并且用递归得到的当前最优值来剪枝。这个套路在逆序对计数、最大子段和、平面凸包的分治版本里都能看到影子。以最大子段和为例分治的做法就是分别求左半、右半的最大子段和再求必然跨越中点的最大子段和三者取最大。跨中点的那一项从中点往两边各扫一遍就能算出来。你看结构和最近点对几乎一模一样递归负责单侧合并负责跨中线跨中线部分要设计一个能在 O(n) 内算完的方法。真正掌握分治的标志不是能默写最近点对的代码而是遇到一个新问题时能立刻反应过来这题的最优解是不是可以按跨越分界线来分类然后设计出合并阶段的做法。这才是这门课想让你学会的东西。最后说一个我自己的习惯。写完分治代码之后我会先拿五个点、十个点的随机数据跑一遍暴力对拍对拍过了再上规模。因为分治的错误几乎全都是小数据看不出来、大数据必现的类型越早用随机数据把 bug 逼出来后面省的时间越多。当年我就是因为跳过了对拍这一步一个midx取值时机的问题被卡了整整两天最后是靠着一组只有 17 个点、坐标范围在 ±3 之间的随机数据才复现出来。坐标范围取小一点重复点和共线情况多反而是最容易暴露边界问题的数据。