「Wallace 笔记」平面最近点对 解法汇总

征求更多的解法，应保证算法及其复杂度正确。

Description - 平面最近点对问题

给定平面上 (n) 个点 ({(x_1, y_1), (x_2, y_2), cdots , (x_n, y_n)})，找出其中的一对点的（欧几里得）距离，使得在这 (n) 个点的所有点对中，该距离为所有点对中最小的。

(1le nle 2 imes 10^5, 0le x_i, y_ile 10^9)。

Solution 1 - 平面分治

Analysis

考虑分治法。先按 (x) 坐标将当前点集 (P = {p_1, p_2, cdots , p_n}) 排序，然后以 (p_{lfloor{n/2} floor}) 的 (x) 坐标为基准将平面分割为两个部分。设这两个部分的点集分别为 (P_1, P_2)。

假如说我们使用递归将左右两部分的点集的最近点对的距离计算好了，分别为 (d_1, d_2)。那么我们设我们将要返回的答案为 (d = min(d_1, d_2))。

现在要做的就是分治算法中最重要的一步之一——合并答案：我们现在已经得到了点对中的两个点同时在 (P_1) 或 (P_2) 中的结果，接下来需要计算两个点分别在 (P_1, P_2) 中的答案。

然而直接二层循环枚举 (P) 中的所有点对的复杂度为平方级别的，需要优化。

首先我们 筛选出可能对答案有贡献的点：当前点集中所有 (x) 坐标与 (p_{lfloor{n/2} floor}) 的 (x) 坐标之差的绝对值 严格小于 (d) 的点，这里记为点集 (Q)。

显然我们的枚举范围已经大大缩小了，不过这还不可以直接枚举。

若我们可以通过某种手段使筛选出的点 (Q) 的 (y) 坐标有序，那么我们可以再次优化：

当枚举到 (Q) 中的某点 (q) 时，我们再次筛选—— 限定只有 (y) 坐标与 (q) 的 (y) 坐标的差值严格小于 (d) 的点才参与计算。

因为只要两点之间的 (y) 坐标之差大于等于 (d)，那么它们之间的距离也必然 (ge d)，显然不可能对答案有贡献。上一次对枚举范围的缩小的依据也是同理。

这里将最终筛选出的点的集合记为 (A)。可以发现点集 (A) 的范围实际上是一个矩形。

如果我们对 (A) 中的点再进行暴力枚举，那么本次递归的计算规模或被控制在 (O(n))。而其中 (|A|) 是 (O(1)) 级别，这个我们稍后再作证明。

对于上述的使 (y) 坐标有序的方法，只要在递归完归并排序即可。

那么总复杂度为 (T(n) = 2T(frac{n}{2}) + O(n) = O(nlog n))。

Proof

这里证明算法复杂度的正确性，即上文 (A) 的大小为常数。

我们根据将当前点 (q) 筛选出的点集 (A) 所对应的矩形分割为两个正方形，那么两个正方形的边长为 (d)。

由于 (d) 是在左右两边的递归的结果，那易知 对于两个正方形中的其中一个而言，其中的所有点对的距离均不超过 (d)。

如果我们再将这两个正方形分割为共 (8) 个小正方形，边长为 (frac{d}{2})。

我们尝试着在一个小正方形中放尽可能多的点，但显然我们连两个点都放不下。

为什么？因为两个点在同一个正方形区域中，要使两者距离最大，那就是 正方形相对的两个顶点，其最大距离就是其构成的对角线长。即 (frac{d}{sqrt{2}})。

而 (d > frac{d}{sqrt{2}}) ，因此 一个小正方形中至多包含一个点。

那么整个矩形就会至多有 (8) 个点，减去 (q)，可得 (max |A| = 7)。

Code

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Article : 平面最近点对 解法汇总
 */
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <iomanip>
#include <vector>

using namespace std;
const int N = 2e5 + 5;

struct Point {
	double x, y;
	void read() {
		scanf("%lf%lf", &x, &y);
	}
} dat[N];
int n;

struct cmp_x {
	bool operator () (const Point& x, const Point& y) {
		return x.x < y.x;
	}
};
struct cmp_y {
	bool operator () (const Point& x, const Point& y) {
		return x.y < y.y;
	}
};
double dist(const Point& x, const Point& y) {
	return sqrt((x.x - y.x) * (x.x - y.x) + (x.y - y.y) * (x.y - y.y));
}

double solve(int l, int r) {
	if (l == r) return 1e18;
	
	int mid = (l + r) >> 1;
	double ans = 1e18, x = dat[mid].x;
	
	ans = min(ans, solve(l, mid));
	ans = min(ans, solve(mid + 1, r));
	inplace_merge(dat + l, dat + mid + 1, dat + r + 1, cmp_y());
	
	vector<Point> rec;
	for (register int i = l; i <= r; i++)
		if (fabs(x - dat[i].x) < ans)
			rec.push_back(dat[i]);
	
	for (register int i = 0; i < rec.size(); i++)
		for (register int j = i + 1; j < rec.size(); j++) {
			if (rec[j].y - rec[i].y >= ans) break;
			ans = min(ans, dist(rec[i], rec[j]));
		}
	return ans;
}

signed main() {
	scanf("%d", &n);
	for (register int i = 1; i <= n; i++)
		dat[i].read();
	
	sort(dat + 1, dat + 1 + n, cmp_x());
	printf("%.4f
", solve(1, n));
	return 0;
}

Note

分治法至今任为求平面最近点对最主流的的算法。

另外，该算法所基于的平面分治算法还可以解决更多相似的平面问题。

Solution 2 - 直接扫描法

Analysis

其实还有一种不需要分治的做法。

这里我们先将输入的点按 (x) 坐标升序排序，然后从左往右扫。

当扫到某一个点 (p) 时，我们可以通过计算其右侧的点与之的最短距离来求解。

首先我们有一个二层循环的平方暴力，然后像上文一样考虑筛出可能对答案有贡献的点以减小枚举范围。

显然，与 (p) 的 (x) 坐标之差超过 (d)（当前答案）的可以直接省掉。

然后我们再次限定搜索范围为 与 (p) 的 (y) 坐标之差不超过 (d) 的点。

可以证明最终筛出的点的个数任为常数级别（证明见下文）。

在这里我们需要维护在当前点 (p) 之前且 (x) 坐标在范围内的点的集合。

不妨使用 multiset 来完成维护的任务——在刚扫到下一个点时，及时在 multiset 中删除 (x) 坐标过小的元素。

如何确定 (y) 坐标符合条件的点的集合？我们可以令 multiset 以 (y) 坐标为第一关键字保持有序，这样就可以通过 lower_bound 函数实现快速找到枚举的起始位置了。

由于每个点都只被（multiset）插入或删除一次，因此这一部分的复杂度为 (O(nlog n))。

若能证明枚举的点集大小为常数级别，那么这里的复杂度也为 (O(nlog n))。（lower_bound 需要一个 log）

加上排序，总复杂度为 (O(nlog n))。

Proof

仍然是复杂度证明。

根据我们限定的 (x, y) 坐标的范围，易知枚举的区域是一个边长为 (d) 和 (2d) 的矩形。

和上面的证明方式相似，我们将整个矩形分为左右两个大正方形，再分别分为 (4) 个边长为 (frac{d}{2}) 的小正方形。

由于每个大正方形中点对的距离不小于 (d)，于是每个小正方形中最大放 (1) 个点。

很显然，矩形区域中的点为常数规模。复杂度于是得证。

Code

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Article : 平面最近点对 解法汇总
 */
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <set>

using namespace std;
const int N = 2e5 + 5;

struct Point {
	double x, y;
	Point(int a = 0, int b = 0) {
		x = a, y = b;
	}
	void read() {
		scanf("%lf%lf", &x, &y);
	}
} dat[N];
int n;

struct cmp_x {
	bool operator () (const Point& x, const Point& y) {
		return x.x != y.x ? x.x < y.x : x.y < y.y;
	}
};
struct cmp_y {
	bool operator () (const Point& x, const Point& y) {
		return x.y != y.y ? x.y < y.y : x.x < y.x;
	}
};
double dist(const Point& x, const Point& y) {
	return sqrt((x.x - y.x) * (x.x - y.x) + (x.y - y.y) * (x.y - y.y));
}

multiset<Point, cmp_y> p;

signed main() {
	scanf("%d", &n);
	for (register int i = 1; i <= n; i++)
		dat[i].read();
	sort(dat + 1, dat + 1 + n, cmp_x());
	
	double ans = 1e18;
	for (register int r = 0, l = 1; r <= n; r++) {
		while (l <= r && dat[l].x - dat[r].x >= ans)
			p.erase(p.find(dat[l++]));
		for (multiset<Point, cmp_y>::iterator it = p.lower_bound(Point(dat[r].x, dat[r].y - ans));
			it != p.end() && it->y - ans < dat[r].y; it++) ans = min(ans, dist(dat[r], *it));
		p.insert(dat[r]);
	}
	
	printf("%.4f
", ans);
}

Note

相对好理解且好写的一种算法，但分治法仍然有学习的必要。

Postscript

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/-Wallace-/p/13268575.html