计算几何——平面最近点对
问题描述
给定平面上 \(n\)（\(n\ge2\)）个点，找出一堆点，使得其间距离最短。
下午将介绍分治做法、非分治做法，以及期望线性做法。
其中运行速度（P1429 平面最近点对 加强版）大致为，非分治做法最快，期望线性做法最慢。
朴素算法
非常显然了，\(\mathcal O(n^2)\) 的遍历每一对点。
非常简略的一个代码示例，
function solve():
min_dist = inf
for i in point_set:
for j in point_set:
min_dist = min(min_dist, dist(i,j))
return min_dist
简单剪枝
考虑一种常见的统计序列的思想：
依次加入每一个元素，统计它和其左边所有元素的贡献。
具体地，
-
我们把所有点按照 \(x_i\) 为第一关键字、\(y_i\) 为第二关键字排序。
-
同时，建立一个以 \(y_i\) 为第一关键字、 \(x_i\) 为第二关键字排序的
multiset。 -
对于每一个位置 \(i\)，我们执行以下操作：
-
假设我们已经算出来的最小距离是 \(d\)，
-
将所有满足 \(|x_i-x_j|\ge d\) 的点从集合中删除。它们不会再对答案有贡献。
-
对于集合内满足 \(|y_i-y_j|< d\) 的所有点，统计它们和 \((x_i,y_i)\) 的距离。
-
将 \((x_i,y_i)\) 插入到集合中。
这个算法的复杂度为 \(\mathcal O(n\log n)\)，比分治做法常数略小，证明略。
代码：
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
struct emm {
ll x, y;
emm() = default;
emm(ll x, ll y): x(x), y(y) {}
friend ll operator *(const emm &a, const emm &b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
};
struct cmp_x {
bool operator ()(const emm &a, const emm &b) const {
return a.x == b.x ? a.y < b.y : a.x < b.x;
}
};
struct cmp_y {
bool operator ()(const emm &a, const emm &b) const {
return a.y < b.y;
}
};
double ans = 1e10;
void upd_ans(const emm &a, const emm &b) {
double dist = sqrt(a * b);
if (dist < ans) ans = dist;
}
vector<emm> a;
signed main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
int n; cin >> n; a.resize(n);
for (int i = 0; i < n; ++i) cin >> a[i].x >> a[i].y;
sort(a.begin(), a.end(), cmp_x());
multiset<emm, cmp_y> s;
for (int i = 0, l = 0; i < n; ++i) {
while (l < i && a[i].x - a[l].x >= ans) s.erase(s.find(a[l++]));
auto it = s.lower_bound(emm(a[i].x, a[i].y - ans));
for (; it != s.end() && it->y - a[i].y < ans; ++it) upd_ans(a[i], *it);
s.insert(a[i]);
}
printf("%.4lf\n", ans);
return 0;
}
这个做法的问题在于 multiset 的大常数，但是好写。
分治算法
考虑如果我们把所有点按照 \(x_i\) 排序，分治解决。
-
假设我们已经算出来的最小距离是 \(d\)。
-
考虑如何合并，显然只有两个集合分界线处各 \(d\) 距离内的点需要考虑。
-
我们枚举这个小集合内的点，计算每个点向下最多 \(d\) 个单位的点的贡献。
因为当前最小距离 \(d\)、向下枚举的是 \(d\times2d\) 的矩阵，其内部的点的个数是 \(\mathcal O(1)\) 的。
因此，整体复杂度即考虑分治的复杂度，即 \(\mathcal O(n\log n)\)，但是常数比非分治略大。
代码：
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using db = double;
struct emm {
ll x, y;
emm() = default;
emm(ll x, ll y): x(x), y(y) {}
friend ll operator *(const emm &a, const emm &b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
};
struct cmp_x {
bool operator ()(const emm &a, const emm &b) const {
return a.x == b.x ? a.y < b.y : a.x < b.x;
}
};
struct cmp_y {
bool operator ()(const emm &a, const emm &b) const {
return a.y < b.y;
}
};
double ans = 1e10;
void upd_ans(const emm &a, const emm &b) {
double dist = sqrt(a * b);
if (dist < ans) ans = dist;
}
vector<emm> a;
void merge(int l, int r) {
if (l == r) return;
int m = l + r >> 1; ll mx = a[m].x;
merge(l, m), merge(m + 1, r);
inplace_merge(a.begin() + l, a.begin() + m + 1, a.begin() + r + 1, cmp_y());
vector<emm> t;
for (int i = l; i <= r; ++i) {
if (abs(a[i].x - mx) >= ans) continue;
for (auto j = t.rbegin(); j != t.rend(); ++j) {
if (a[i].y - j->y >= ans) break;
upd_ans(a[i], *j);
}
t.push_back(a[i]);
}
}
signed main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
int n; cin >> n; a.resize(n);
for (int i = 0; i < n; ++i) cin >> a[i].x >> a[i].y;
sort(a.begin(), a.end(), cmp_x());
merge(0, n - 1);
printf("%.4lf\n", ans);
return 0;
}
使用了 std::inplace_merge 作为归并，详见 cppreference。
期望线性
注意是期望线性做法，复杂度理论期望值是 \(\mathcal O(n)\) 的。
但是实际上常数巨大，而且容易被卡，实测速度反而最慢。
-
同样我们考虑加入一个点的贡献，但是这里需要先随机打乱。
-
记前 \(i-1\) 个点的最近点对距离为 \(d\)，将平面以 \(d\) 为边长划分为若干个网格。
-
检查第 \(i\) 个点所在网格的周围九个网格中的所有点，并更新答案。
-
使用哈希表存下每个网格内的点，如果答案被更新，就重构网格图，否则不重构。
因为前 \(i-1\) 个点的最近点对距离为 \(d\)，从而每个网格不超过 \(4\) 个点。
注意到需检查的点的个数是 \(\mathcal O(1)\) 的，在前 \(i\) 个点中，最近点对包含 \(i\) 的概率为
\(\mathcal O(1/i)\)。
而重构网格的代价为 \(\mathcal O(i)\)，从而第 \(i\) 个点的期望代价为 \(\mathcal O(1)\)。
于是对于 \(n\) 个点，该算法期望为 \(\mathcal O(n)\)。
代码：
#include <bits/stdc++.h>
using namespace std;
struct my_hash {
static uint64_t splitmix64(uint64_t x) {
x += 0x9e3779b97f4a7c15;
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
return x ^ (x >> 31);
}
size_t operator()(uint64_t x) const {
static const uint64_t FIXED_RANDOM =
chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
size_t operator()(pair<uint64_t, uint64_t> x) const {
static const uint64_t FIXED_RANDOM =
chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x.first + FIXED_RANDOM) ^
(splitmix64(x.second + FIXED_RANDOM) >> 1);
}
};
mt19937 rng(time(0));
using ll = long long;
using grid = pair<int, int>;
struct emm {
ll x, y;
emm() = default;
emm(ll x, ll y): x(x), y(y) {}
friend ll operator *(const emm &a, const emm &b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
};
double ans = 1e10;
void upd_ans(const emm &a, const emm &b) {
double dist = sqrt(a * b);
if (dist < ans) ans = dist;
}
vector<emm> a;
unordered_map<grid, vector<emm>, my_hash> ump;
#define group(e, t) make_pair(e.x / (int)t, e.y / (int)t)
signed main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
int n; cin >> n; a.resize(n);
for (int i = 0; i < n; ++i) cin >> a[i].x >> a[i].y;
shuffle(a.begin(), a.end(), rng);
for (int i = 0; i < n; ++i) {
double lt = ans;
int tx, ty; tie(tx, ty) = group(a[i], lt);
for (int kx = tx - 1; kx <= tx + 1; ++kx) {
for (int ky = ty - 1; ky <= ty + 1; ++ky) {
auto eq = make_pair(kx, ky);
if (!ump.count(eq)) continue;
for (emm j : ump[eq]) upd_ans(a[i], j);
}
}
if (ans == 0) break;
if (ans != lt) {
ump = decltype(ump)();
for (int j = 0; j < i; ++j) ump[group(a[j], ans)].push_back(a[j]);
}
ump[group(a[i], ans)].push_back(a[i]);
}
printf("%.4lf\n", ans);
return 0;
}
这个算法的常数很大，主要在于哈希（代码里手写了哈希函数）。
标签：const,emm,ll,几何,int,最近,ans,平面,dist From： https://www.cnblogs.com/RainPPR/p/18208157