计算几何——平面最近点对

问题描述

给定平面上 \(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;
}

这个算法的常数很大，主要在于哈希（代码里手写了哈希函数）。

posted @ 2024-05-23 12:23 RainPPR 阅读(15) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

RainPPR（ツ）

逄沛然

计算几何——平面最近点对

计算几何——平面最近点对

问题描述

朴素算法

简单剪枝

分治算法

期望线性

公告