|
| 1 | +--- |
| 2 | +categories: |
| 3 | +- OI |
| 4 | +title: 期望线性的平面最近点对算法 |
| 5 | +date: 2024-09-28 23:19:00 |
| 6 | +tags: |
| 7 | +- geometry |
| 8 | +--- |
| 9 | + |
| 10 | +讲述一种**期望**线性复杂度的平面最近点对算法。 |
| 11 | + |
| 12 | +1. 将点打乱 |
| 13 | +2. 对于小常数 `$D$`,暴力计算前 `$D$` 个点的平面最近点对。 |
| 14 | +3. 考虑从前 `$i-1$` 个点推出前 `$i$` 个点的平面最近点对: |
| 15 | + + 设前 `$i-1$` 个点的平面最近点对距离为 `$s$`,将平面以 `$s$` 为边长划分成若干网格,用哈希表记录每个网格内的点。 |
| 16 | + + 检查第 `$i$` 个点所在的网格及其周围共 `$9$` 个网格内的点是否能更新答案,检查的点的数量是 `$O(1)$` 的,因为每个网格最多有 `$4$` 个点。 |
| 17 | + + 若答案更新,重构网格。前 `$i$` 个点的平面最近点对包含第 `$i$` 个点的概率为 `$O\left(\frac{1}{i}\right)$`,重构网格的代价为 `$O(i)$`,故每个点的期望代价为 `$O(1)$`。 |
| 18 | + |
| 19 | +若给定的点有重,似乎算法效率会减小,建议特判。 |
| 20 | + |
| 21 | +更奇怪的是,排序的效率似乎高于打乱的效率。 |
| 22 | + |
| 23 | +```cpp |
| 24 | +#include <algorithm> |
| 25 | +#include <cmath> |
| 26 | +#include <cstdint> |
| 27 | +#include <iomanip> |
| 28 | +#include <iostream> |
| 29 | +#include <random> |
| 30 | +#include <unordered_map> |
| 31 | +#include <unordered_set> |
| 32 | +#include <utility> |
| 33 | +#include <vector> |
| 34 | +#define x first |
| 35 | +#define y second |
| 36 | +using namespace std; |
| 37 | + |
| 38 | +const int M = 1000; |
| 39 | +int n; |
| 40 | +double ans = HUGE_VAL; |
| 41 | +vector<pair<int, int>> vec; |
| 42 | +unordered_map<uint64_t, vector<pair<int, int>>> grid; |
| 43 | + |
| 44 | +inline uint64_t mapping(unsigned a, unsigned b) { |
| 45 | + return (uint64_t)a << 32 | (uint64_t)b; |
| 46 | +} |
| 47 | + |
| 48 | +void prework() { |
| 49 | + int m = min(n, M); |
| 50 | + for (int i = 0; i < m; ++i) |
| 51 | + for (int j = i + 1; j < m; ++j) |
| 52 | + ans = min(ans, hypot(vec[i].x - vec[j].x, vec[j].y - vec[i].y)); |
| 53 | +} |
| 54 | + |
| 55 | +void remake(int m) { |
| 56 | + grid.clear(); |
| 57 | + for (int i = 0; i <= m; ++i) { |
| 58 | + int a = floor(vec[i].x / ceil(ans)), b = floor(vec[i].y / ceil(ans)); |
| 59 | + grid[mapping(a, b)].push_back(make_pair(vec[i].x, vec[i].y)); |
| 60 | + } |
| 61 | +} |
| 62 | + |
| 63 | +bool has_same() { |
| 64 | + unordered_set<uint64_t> ust; |
| 65 | + for (auto &&[u, v] : vec) { |
| 66 | + uint64_t hs = mapping(u, v); |
| 67 | + if (!ust.insert(hs).second) |
| 68 | + return true; |
| 69 | + } |
| 70 | + return false; |
| 71 | +} |
| 72 | + |
| 73 | +int main() { |
| 74 | + mt19937 rnd; rnd.seed(random_device()()); |
| 75 | + ios::sync_with_stdio(false); cin.tie(0); |
| 76 | + cin >> n; |
| 77 | + for (int i = 0; i < n; ++i) { |
| 78 | + int a, b; |
| 79 | + cin >> a >> b; |
| 80 | + vec.emplace_back(a, b); |
| 81 | + } |
| 82 | + if (has_same()) { |
| 83 | + ans = 0; |
| 84 | + goto jump; |
| 85 | + } |
| 86 | + shuffle(vec.begin(), vec.end(), rnd); |
| 87 | + prework(); remake(min(n, M) - 1); |
| 88 | + for (int i = M; i < n; ++i) { |
| 89 | + int a = floor(vec[i].x / ceil(ans)), b = floor(vec[i].y / ceil(ans)); |
| 90 | + double ret = ans; |
| 91 | + for (int dx = -1; dx <= 1; ++dx) |
| 92 | + for (int dy = -1; dy <= 1; ++dy) { |
| 93 | + int s = a + dx, t = b + dy; |
| 94 | + if (grid.count(mapping(s, t)) == 0) |
| 95 | + continue; |
| 96 | + for (auto &&[u, v] : grid[mapping(s, t)]) |
| 97 | + ret = min(ret, hypot(vec[i].x - u, vec[i].y - v)); |
| 98 | + } |
| 99 | + if (ret < ans) |
| 100 | + ans = ret, remake(i); |
| 101 | + else |
| 102 | + grid[mapping(a, b)].emplace_back(vec[i].x, vec[i].y); |
| 103 | + } |
| 104 | +jump: |
| 105 | + cout << fixed << setprecision(4) << ans << "\n"; |
| 106 | + return 0; |
| 107 | +} |
| 108 | +``` |
| 109 | +
|
| 110 | +也可以参考我在 [Qoj 上的提交记录](https://qoj.ac/submission/1017630)。 |
0 commit comments