平面最近点对 & 最小周长三角形 & 曼哈顿距离最近
Statement
求平面最近点对的距离，距离定义为欧几里德距离。
Solution
考虑按 \(x\) 排序，分治计算
先往左右递归，设得到的答案为 \(d\)，当前算跨过中点的答案，那么答案 \(\ge d\) 的点对可以不用枚举
设中点为 \(m\)
- 对 \(i\in[l..m]\)，\(x_i\le x_m-d\) 的不用枚举；对 \(j\in[m+1..r]\)，\(x_j\ge x_m+d\) 的不用枚举
- 对于所有 \(i\in[l..m]\)，不超过 6 个 \(j\in[m+1..r]\) 满足 \(x_j\in[x_m,x_m+d]\) 且 \(y_j\in[y_i-d,y_i+d]\)
于是对于需要枚举的点，按 \(y\) 排序，左边每个 \(i\) 枚举不超过 \(6\) 个候选点更新答案，\(O(n\log n)\)
Code
#include <bits/stdc++.h>
using namespace std;
#define rep(i, j, k) for (int i = (j); i <= (k); ++i)
#define reo(i, j, k) for (int i = (j); i >= (k); --i)
typedef long long ll;
const int N = 5e5 + 10;
struct Point {
ll x, y;
} a[N], b[N];
int n;
ll ans;
ll dis(int u, int v) {
return (a[u].x - a[v].x) * (a[u].x - a[v].x) + (a[u].y - a[v].y) * (a[u].y - a[v].y);
}
ll dis(Point u, Point v) {
return (u.x - v.x) * (u.x - v.x) + (u.y - v.y) * (u.y - v.y);
}
void Solve(int l, int r) {
if (r - l < 3) {
rep(i, l, r) rep(j, i + 1, r) ans = min(ans, dis(i, j));
sort(a + l, a + r + 1, [&](const Point& u, const Point& v) {
return u.y < v.y;
});
return;
}
int mid = (l + r) >> 1;
ll X = a[mid].x;
Solve(l, mid), Solve(mid + 1, r);
ll d = ceil(sqrt(ans));
vector<Point> L, R;
rep(i, l, mid)
if (a[i].x > X - d) L.push_back(a[i]);
rep(i, mid + 1, r)
if (a[i].x < X + d) R.push_back(a[i]);
int sz = R.size(), pl = 0, pr = 0;
if (sz) {
for (auto p : L) {
while (pl < sz && R[pl].y < p.y - d) ++pl;
while (pr < sz && R[pr].y <= p.y + d) ++pr;
pr = max(pr - 1, 0);
rep(i, pl, pr) ans = min(ans, dis(p, R[i]));
}
}
for (int i = l, j = mid + 1, tot = l; i <= mid || j <= r; )
if (i <= mid && (j > r || a[i].y <= a[j].y)) b[tot++] = a[i++];
else b[tot++] = a[j++];
rep(i, l, r) a[i] = b[i];
}
int main() {
ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> n;
rep(i, 1, n) {
cin >> a[i].x >> a[i].y;
}
sort(a + 1, a + n + 1, [&](const Point& u, const Point& v) {
return u.x < v.x;
});
ans = 1e18;
Solve(1, n);
cout << ans << '\n';
return 0;
}
Statement
找出最小周长三角形的周长。为减小难度，这里三角形也包含共线的三点。
Solution
这里是距离 \(\ge\frac d2\) 的不用枚举，注意到这里候选点数同样是常数。
Code
#include <bits/stdc++.h>
using namespace std;
#define rep(i, j, k) for (int i = (j); i <= (k); ++i)
#define reo(i, j, k) for (int i = (j); i >= (k); --i)
typedef long long ll;
const int N = 2e5 + 10;
struct Point {
double x, y;
} a[N];
int n;
double ans = 1e18;
double dis(int u, int v) {
return sqrt(1.0 * (a[u].x - a[v].x) * (a[u].x - a[v].x) + 1.0 * (a[u].y - a[v].y) * (a[u].y - a[v].y));
}
double dis(Point u, Point v) {
return sqrt(1.0 * (u.x - v.x) * (u.x - v.x) + 1.0 * (u.y - v.y) * (u.y - v.y));
}
void Solve(int l, int r) {
if (r - l < 3) {
if (r - l == 2)
ans = min(ans, dis(l, l + 1) + dis(l, l + 2) + dis(l + 1, l + 2));
return;
}
int mid = (l + r) >> 1;
Solve(l, mid), Solve(mid + 1, r);
double d = ans / 2.0;
vector<Point> L, R;
rep(i, l, mid)
if (a[i].x >= a[mid].x - d) L.push_back(a[i]);
rep(i, mid + 1, r)
if (a[i].x <= a[mid].x + d) R.push_back(a[i]);
sort(L.begin(), L.end(), [&](const Point& u, const Point& v) {
return u.y < v.y;
});
sort(R.begin(), R.end(), [&](const Point& u, const Point& v) {
return u.y < v.y;
});
int sz = R.size(), pl = 0, pr = 0;
for (auto i : L) {
while (pl < sz && R[pl].y < i.y - d) ++pl;
while (pr < sz && R[pr].y <= i.y + d) ++pr;
pr = max(0, pr - 1);
if (pr - pl + 1 >= 2)
rep(u, pl, pr)
rep(v, u + 1, pr) {
ans = min(ans, dis(i, R[u]) + dis(i, R[v]) + dis(R[u], R[v]));
}
}
sz = L.size(), pl = 0, pr = 0;
for (auto i : R) {
while (pl < sz && L[pl].y < i.y - d) ++pl;
while (pr < sz && L[pr].y <= i.y + d) ++pr;
pr = max(0, pr - 1);
if (pr - pl + 1 >= 2)
rep(u, pl, pr)
rep(v, u + 1, pr) {
ans = min(ans, dis(i, L[u]) + dis(i, L[v]) + dis(L[u], L[v]));
}
}
}
int main() {
ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> n;
rep(i, 1, n) {
cin >> a[i].x >> a[i].y;
}
sort(a + 1, a + n + 1, [&](const Point& u, const Point& v) {
return u.x < v.x;
});
Solve(1, n);
cout << fixed << setprecision(6) << ans << '\n';
return 0;
}
Statement
- 加点
- 求曼哈顿距离与 \(u\) 最近的距离
Solution
问题等价于三维偏序
因为若 \(x_i\ge x,y_i\ge y\)，\(\min_{i}\{|x_i-x|+|y_i-y|\}=\min\{x_i+y_i-(x+y)\}=\min\{x_i+y_i\}-(x+y)\)
Code
#include <bits/stdc++.h>
using namespace std;
#define rep(i, j, k) for (int i = (j); i <= (k); ++i)
#define reo(i, j, k) for (int i = (j); i >= (k); --i)
typedef long long ll;
const int N = 3e6 + 10;
struct Oper {
int op, x, y, id;
} a[N], b[N], c[N];
int n, m, tot, qtot, ans[N];
struct BIT {
int mn[N];
BIT() {
rep(i, 0, N - 1) mn[i] = 1e9;
}
void Upd(int x, int v) {
for (; x < N; x += x & -x) mn[x] = min(mn[x], v);
}
int Qry(int x) {
int res = 1e9;
for (; x; x -= x & -x) res = min(res, mn[x]);
return res;
}
void Clear(int x) {
for (; x < N; x += x & -x) mn[x] = 1e9;
}
} bit;
void Solve1(int l, int r) {
if (l == r) return;
int mid = (l + r) >> 1;
Solve1(l, mid), Solve1(mid + 1, r);
vector<Oper> Left, Right;
rep(i, l, mid)
if (b[i].op == 1) Left.push_back(b[i]);
rep(i, mid + 1, r)
if (b[i].op == 2) Right.push_back(b[i]);
int pos = Left.size();
reo(i, (int)Right.size() - 1, 0) {
while (pos && Left[pos - 1].x >= Right[i].x)
--pos, bit.Upd(N - Left[pos].y, Left[pos].x + Left[pos].y);
ans[Right[i].id] = min(ans[Right[i].id], bit.Qry(N - Right[i].y) - (Right[i].x + Right[i].y));
}
rep(i, pos, (int)Left.size() - 1) bit.Clear(N - Left[i].y);
int R = l - 1;
for (int i = l, j = mid + 1; i <= mid || j <= r; )
if (i <= mid && (j > r || b[i].x <= b[j].x)) c[++R] = b[i++];
else c[++R] = b[j++];
rep(i, l, r) b[i] = c[i];
}
void Solve2(int l, int r) {
if (l == r) return;
int mid = (l + r) >> 1;
Solve2(l, mid), Solve2(mid + 1, r);
vector<Oper> Left, Right;
rep(i, l, mid)
if (b[i].op == 1) Left.push_back(b[i]);
rep(i, mid + 1, r)
if (b[i].op == 2) Right.push_back(b[i]);
int pos = Left.size();
reo(i, (int)Right.size() - 1, 0) {
while (pos && Left[pos - 1].x >= Right[i].x)
--pos, bit.Upd(Left[pos].y, Left[pos].x - Left[pos].y);
ans[Right[i].id] = min(ans[Right[i].id], bit.Qry(Right[i].y) - (Right[i].x - Right[i].y));
}
rep(i, pos, (int)Left.size() - 1) bit.Clear(Left[i].y);
int R = l - 1;
for (int i = l, j = mid + 1; i <= mid || j <= r; )
if (i <= mid && (j > r || b[i].x <= b[j].x)) c[++R] = b[i++];
else c[++R] = b[j++];
rep(i, l, r) b[i] = c[i];
}
void Solve3(int l, int r) {
if (l == r) return;
int mid = (l + r) >> 1;
Solve3(l, mid), Solve3(mid + 1, r);
vector<Oper> Left, Right;
rep(i, l, mid)
if (b[i].op == 1) Left.push_back(b[i]);
rep(i, mid + 1, r)
if (b[i].op == 2) Right.push_back(b[i]);
int pos = -1;
rep(i, 0, (int)Right.size() - 1) {
while (pos < (int)Left.size() - 1 && Left[pos + 1].x <= Right[i].x)
++pos, bit.Upd(Left[pos].y, - Left[pos].x - Left[pos].y);
ans[Right[i].id] = min(ans[Right[i].id], Right[i].x + Right[i].y + bit.Qry(Right[i].y));
}
rep(i, 0, pos) bit.Clear(Left[i].y);
int R = l - 1;
for (int i = l, j = mid + 1; i <= mid || j <= r; )
if (i <= mid && (j > r || b[i].x <= b[j].x)) c[++R] = b[i++];
else c[++R] = b[j++];
rep(i, l, r) b[i] = c[i];
}
void Solve4(int l, int r) {
if (l == r) return;
int mid = (l + r) >> 1;
Solve4(l, mid), Solve4(mid + 1, r);
vector<Oper> Left, Right;
rep(i, l, mid)
if (b[i].op == 1) Left.push_back(b[i]);
rep(i, mid + 1, r)
if (b[i].op == 2) Right.push_back(b[i]);
int pos = -1;
rep(i, 0, (int)Right.size() - 1) {
while (pos < (int)Left.size() - 1 && Left[pos + 1].x <= Right[i].x)
++pos, bit.Upd(N - Left[pos].y, Left[pos].y - Left[pos].x);
ans[Right[i].id] = min(ans[Right[i].id], Right[i].x - Right[i].y + bit.Qry(N - Right[i].y));
}
rep(i, 0, pos) bit.Clear(N - Left[i].y);
int R = l - 1;
for (int i = l, j = mid + 1; i <= mid || j <= r; )
if (i <= mid && (j > r || b[i].x <= b[j].x)) c[++R] = b[i++];
else c[++R] = b[j++];
rep(i, l, r) b[i] = c[i];
}
int main() {
ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> n >> m;
rep(i, 1, n)
++tot, a[tot].op = 1, cin >> a[tot].x >> a[tot].y;
rep(i, 1, m) {
++tot, cin >> a[tot].op >> a[tot].x >> a[tot].y;
if (a[tot].op == 2)
a[tot].id = ++qtot;
}
rep(i, 1, qtot) ans[i] = 2e9;
rep(i, 1, tot) b[i] = a[i];
Solve1(1, tot);
rep(i, 1, tot) b[i] = a[i];
Solve2(1, tot);
rep(i, 1, tot) b[i] = a[i];
Solve3(1, tot);
rep(i, 1, tot) b[i] = a[i];
Solve4(1, tot);
rep(i, 1, qtot) cout << ans[i] << '\n';
return 0;
}
注意到 K-D Tree 直接秒了。
