euclidean_mst.hpp¶
通过 Delaunay 候选边求整数平面点的欧氏最小生成树；坐标差须能用 64 位有符号整数表示，距离平方、叉积与四次圆内判定须能用 128 位有符号整数表示。
Complexity: Time: \(O(n \log n)\). Space: \(O(n)\).
Source: Library Checker
AC 记录：euclidean_mst。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log n).
/// Space: O(n).
#include "atcoder/dsu.hpp"
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
namespace euclidean_mst_internal {
template <class Int = long long, class Int2 = long long> struct VecI2 {
Int x, y;
VecI2() : x(0), y(0) {}
VecI2(std::pair<Int, Int> _p)
: x(std::move(_p.first)), y(std::move(_p.second)) {}
VecI2(Int _x, Int _y) : x(std::move(_x)), y(std::move(_y)) {}
VecI2 &operator+=(VecI2 r) {
x += r.x;
y += r.y;
return *this;
}
VecI2 &operator-=(VecI2 r) {
x -= r.x;
y -= r.y;
return *this;
}
VecI2 &operator*=(Int r) {
x *= r;
y *= r;
return *this;
}
VecI2 operator+(VecI2 r) const { return VecI2(x + r.x, y + r.y); }
VecI2 operator-(VecI2 r) const { return VecI2(x - r.x, y - r.y); }
VecI2 operator*(Int r) const { return VecI2(x * r, y * r); }
VecI2 operator-() const { return VecI2(-x, -y); }
Int2 operator*(VecI2 r) const {
return Int2(x) * Int2(r.x) + Int2(y) * Int2(r.y);
}
Int2 operator^(VecI2 r) const {
return Int2(x) * Int2(r.y) - Int2(y) * Int2(r.x);
}
bool operator<(VecI2 r) const { return x < r.x || (!(r.x < x) && y < r.y); }
Int2 norm() const { return Int2(x) * Int2(x) + Int2(y) * Int2(y); }
static bool compareYX(VecI2 a, VecI2 b) {
return a.y < b.y || (!(b.y < a.y) && a.x < b.x);
}
static bool compareXY(VecI2 a, VecI2 b) {
return a.x < b.x || (!(b.x < a.x) && a.y < b.y);
}
bool operator==(VecI2 r) const { return x == r.x && y == r.y; }
bool operator!=(VecI2 r) const { return x != r.x || y != r.y; }
};
// Int3 must be able to handle the value range :
// |x| <= | (any input - any input) ** 4 * 12 |
template <class Int = long long, class Int2 = long long, class Int3 = Int2>
class DelaunayTriangulation {
public:
using GPos2 = VecI2<Int, Int2>;
struct Edge {
int to;
int ccw;
int cw;
int rev;
bool on = false;
};
private:
static int isDinOABC(GPos2 a, GPos2 b, GPos2 c, GPos2 d) {
a = a - d;
b = b - d;
c = c - d;
auto val = Int3(b ^ c) * Int3(a.norm()) + Int3(c ^ a) * Int3(b.norm()) +
Int3(a ^ b) * Int3(c.norm());
return val > Int3(0) ? 1 : 0;
}
int getOpenAddress() {
if (oA.empty()) {
es.push_back({});
return (int)es.size() - 1;
}
int res = oA.back();
oA.pop_back();
return res;
}
std::pair<int, int> newEdge(int u, int v) {
int euv = getOpenAddress();
int evu = getOpenAddress();
es[euv].ccw = es[euv].cw = euv;
es[evu].ccw = es[evu].cw = evu;
es[euv].to = v;
es[evu].to = u;
es[euv].rev = evu;
es[evu].rev = euv;
es[euv].on = true;
es[evu].on = true;
return {euv, evu};
}
void eraseSingleEdge(int e) {
int ec = es[e].ccw;
int ecw = es[e].cw;
es[ec].cw = ecw;
es[ecw].ccw = ec;
es[e].on = false;
}
void eraseEdgeBidirectional(int e) {
int ex = es[e].rev;
eraseSingleEdge(e);
eraseSingleEdge(ex);
oA.push_back(e);
oA.push_back(ex);
}
void insertCcwAfter(int e, int x) {
int xc = es[x].ccw;
es[e].ccw = xc;
es[xc].cw = e;
es[e].cw = x;
es[x].ccw = e;
}
void insertCwAfter(int e, int x) {
int xcw = es[x].cw;
es[e].cw = xcw;
es[xcw].ccw = e;
es[e].ccw = x;
es[x].cw = e;
}
// move from ab to ac ... is this ccw?
int isCcw(int a, int b, int c) const {
auto ab = pos[b] - pos[a];
auto ac = pos[c] - pos[a];
auto cp = ab ^ ac;
if (0 < cp)
return 1;
if (cp < 0)
return -1;
return 0;
}
std::pair<int, int> goNext(int, int ea) {
int ap = es[ea].to;
int eap = es[es[ea].rev].ccw;
return {ap, eap};
}
std::pair<int, int> goPrev(int, int ea) {
int ap = es[es[ea].cw].to;
int eap = es[es[ea].cw].rev;
return {ap, eap};
}
std::tuple<int, int, int, int> goBottom(int a, int ea, int b, int eb) {
while (true) {
auto [ap, eap] = goPrev(a, ea);
if (isCcw(b, a, ap) > 0) {
std::tie(a, ea) = {ap, eap};
continue;
}
auto [bp, ebp] = goNext(b, eb);
if (isCcw(a, b, bp) < 0) {
std::tie(b, eb) = {bp, ebp};
continue;
}
break;
}
return {a, ea, b, eb};
}
std::pair<int, int> getMaximum(int a, int ea, bool mn) {
std::pair<int, int> ans = {a, ea};
int p = a, ep = ea;
do {
std::tie(p, ep) = goNext(p, ep);
if (mn)
ans = std::min(ans, std::make_pair(p, ep));
else
ans = std::max(ans, std::make_pair(p, ep));
} while (ep != ea);
return ans;
}
bool isDinOABC(int a, int b, int c, int d) {
return isDinOABC(pos[a], pos[b], pos[c], pos[d]);
}
std::pair<int, int> dfs(int a, int ea, int b, int eb) {
std::tie(a, ea) = getMaximum(a, ea, false);
std::tie(b, eb) = getMaximum(b, eb, true);
auto [al, eal, bl, ebl] = goBottom(a, ea, b, eb);
auto [bu, ebu, au, eau] = goBottom(b, eb, a, ea);
ebl = es[ebl].cw;
ebu = es[ebu].cw;
auto [abl, bal] = newEdge(al, bl);
insertCwAfter(abl, eal);
insertCcwAfter(bal, ebl);
if (al == au)
eau = abl;
if (bl == bu)
ebu = bal;
int ap = al, eap = eal;
int bp = bl, ebp = ebl;
while (ap != au || bp != bu) {
int a2 = es[eap].to;
int b2 = es[ebp].to;
int na = es[eap].ccw;
int nb = es[ebp].cw;
if (eap != eau && na != abl) {
int a1 = es[na].to;
if (isDinOABC(ap, bp, a2, a1)) {
eraseEdgeBidirectional(eap);
eap = na;
continue;
}
}
if (ebp != ebu && nb != bal) {
int b1 = es[nb].to;
if (isDinOABC(b2, ap, bp, b1)) {
eraseEdgeBidirectional(ebp);
ebp = nb;
continue;
}
}
bool ca = ebp == ebu;
if (eap != eau && ebp != ebu) {
if (isCcw(ap, bp, b2) < 0)
ca = true;
else if (isCcw(a2, ap, bp) < 0)
ca = false;
else
ca = isDinOABC(ap, bp, b2, a2);
}
if (ca) {
na = es[es[eap].rev].ccw;
auto [hab, hba] = newEdge(a2, bp);
insertCwAfter(hab, na);
insertCcwAfter(hba, ebp);
eap = na;
ap = a2;
} else {
nb = es[es[ebp].rev].cw;
auto [hba, hab] = newEdge(b2, ap);
insertCcwAfter(hba, nb);
insertCwAfter(hab, eap);
ebp = nb;
bp = b2;
}
}
return {al, abl};
}
std::pair<int, int> solveRange(int l, int r) {
if (r - l == 2) {
int u = l;
int v = l + 1;
auto [uv, vu] = newEdge(u, v);
return {u, uv};
}
if (r - l == 3) {
int u = l;
int v = l + 1;
int w = l + 2;
auto [uv, vu] = newEdge(u, v);
auto [vw, wv] = newEdge(v, w);
int ccw = isCcw(u, v, w);
if (ccw == 0) {
insertCcwAfter(vu, vw);
}
if (ccw > 0) {
auto [uw, wu] = newEdge(u, w);
insertCwAfter(uv, uw);
insertCwAfter(vw, vu);
insertCwAfter(wu, wv);
return {u, uv};
}
if (ccw < 0) {
auto [uw, wu] = newEdge(u, w);
insertCcwAfter(uv, uw);
insertCcwAfter(vw, vu);
insertCcwAfter(wu, wv);
return {v, vu};
}
return {u, uv};
}
int m = (l + r) / 2;
auto [a, ea] = solveRange(l, m);
auto [b, eb] = solveRange(m, r);
return dfs(a, ea, b, eb);
}
void solve() {
int sz = (int)pos.size();
if (sz <= 1)
return;
std::vector<int> pi(pos.size());
for (int i = 0; i < (int)pi.size(); i++)
pi[i] = i;
std::stable_sort(pi.begin(), pi.end(), [&](int l, int r) {
return pos[l].x != pos[r].x ? pos[l].x < pos[r].x : pos[l].y < pos[r].y;
});
auto bu1 = pos;
int ptr = 0;
mp.assign(sz, 0);
for (int i = 0; i < sz; i++) {
int v = pi[i];
if (i == 0 || !(bu1[pi[ptr - 1]] == bu1[v])) {
pi[ptr] = v;
pos[ptr++] = bu1[v];
mp[v] = v;
} else {
mp[v] = pi[ptr - 1];
}
}
if (ptr >= 2)
oOE = solveRange(0, ptr).second;
std::swap(pos, bu1);
for (auto &e : es)
e.to = pi[e.to];
}
std::vector<int> oA;
std::vector<GPos2> pos;
std::vector<Edge> es;
std::vector<int> mp;
int oOE = -1;
public:
DelaunayTriangulation() : pos() { solve(); }
DelaunayTriangulation(std::vector<GPos2> xp) : pos(std::move(xp)) { solve(); }
std::vector<std::pair<int, int>> getEdges() const {
std::vector<std::pair<int, int>> res;
for (int e = 0; e < (int)es.size(); e++)
if (es[e].on) {
int re = es[e].rev;
if (e < re)
continue;
res.push_back({es[e].to, es[re].to});
}
for (int v = 0; v < int(mp.size()); v++) {
if (mp[v] != v)
res.push_back({v, mp[v]});
}
return res;
}
};
} // namespace euclidean_mst_internal
/// @brief Return a Euclidean minimum spanning tree on integral planar points.
/// Divide-and-conquer Delaunay triangulation keeps only O(n) candidate edges:
/// the empty-circumcircle property guarantees that every Euclidean MST edge is
/// present. Kruskal on squared lengths then selects the tree; duplicate points
/// are connected by explicit zero-length edges. Coordinate differences must
/// fit in signed 64 bits; squared lengths, cross products and fourth-degree
/// incircle predicate intermediates must fit in signed 128 bits.
inline std::vector<std::pair<int, int>>
euclidean_mst(const std::vector<point<long long>> &pt) {
using internal_point = euclidean_mst_internal::VecI2<long long, __int128_t>;
using triangulation =
euclidean_mst_internal::DelaunayTriangulation<long long, __int128_t,
__int128_t>;
std::vector<internal_point> lhs;
lhs.reserve(pt.size());
for (const auto &va1 : pt) {
lhs.emplace_back(va1.x, va1.y);
}
struct weighted_edge {
std::pair<int, int> arc;
__int128_t len;
};
std::vector<weighted_edge> ord;
{
auto can = triangulation(std::move(lhs)).getEdges();
ord.reserve(can.size());
for (auto arc : can) {
__int128_t dx = __int128_t(pt[arc.first].x) - pt[arc.second].x;
__int128_t dy = __int128_t(pt[arc.first].y) - pt[arc.second].y;
ord.push_back({arc, dx * dx + dy * dy});
}
}
// Cache squared lengths so sorting never revisits point coordinates.
std::stable_sort(ord.begin(), ord.end(), [](const auto &a, const auto &b) {
return a.len < b.len;
});
atcoder::dsu ds1(int(pt.size()));
std::vector<std::pair<int, int>> ret;
ret.reserve(pt.empty() ? 0 : pt.size() - 1);
for (const auto &edge : ord) {
auto arc = edge.arc;
if (ds1.same(arc.first, arc.second)) {
continue;
}
ds1.merge(arc.first, arc.second);
if (arc.first > arc.second) {
std::swap(arc.first, arc.second);
}
ret.push_back(arc);
}
std::sort(ret.begin(), ret.end());
return ret;
}
} // namespace noya
// NOYA_EUCLIDEAN_MST_HPP
#ifndef NOYA_EUCLIDEAN_MST_HPP
#define NOYA_EUCLIDEAN_MST_HPP 1
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log n).
/// Space: O(n).
#include "atcoder/dsu.hpp"
#include "noya/geometry/geometry_base.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
namespace euclidean_mst_internal {
template <class Int = long long, class Int2 = long long> struct VecI2 {
Int x, y;
VecI2() : x(0), y(0) {}
VecI2(std::pair<Int, Int> _p)
: x(std::move(_p.first)), y(std::move(_p.second)) {}
VecI2(Int _x, Int _y) : x(std::move(_x)), y(std::move(_y)) {}
VecI2 &operator+=(VecI2 r) {
x += r.x;
y += r.y;
return *this;
}
VecI2 &operator-=(VecI2 r) {
x -= r.x;
y -= r.y;
return *this;
}
VecI2 &operator*=(Int r) {
x *= r;
y *= r;
return *this;
}
VecI2 operator+(VecI2 r) const { return VecI2(x + r.x, y + r.y); }
VecI2 operator-(VecI2 r) const { return VecI2(x - r.x, y - r.y); }
VecI2 operator*(Int r) const { return VecI2(x * r, y * r); }
VecI2 operator-() const { return VecI2(-x, -y); }
Int2 operator*(VecI2 r) const {
return Int2(x) * Int2(r.x) + Int2(y) * Int2(r.y);
}
Int2 operator^(VecI2 r) const {
return Int2(x) * Int2(r.y) - Int2(y) * Int2(r.x);
}
bool operator<(VecI2 r) const { return x < r.x || (!(r.x < x) && y < r.y); }
Int2 norm() const { return Int2(x) * Int2(x) + Int2(y) * Int2(y); }
static bool compareYX(VecI2 a, VecI2 b) {
return a.y < b.y || (!(b.y < a.y) && a.x < b.x);
}
static bool compareXY(VecI2 a, VecI2 b) {
return a.x < b.x || (!(b.x < a.x) && a.y < b.y);
}
bool operator==(VecI2 r) const { return x == r.x && y == r.y; }
bool operator!=(VecI2 r) const { return x != r.x || y != r.y; }
};
// Int3 must be able to handle the value range :
// |x| <= | (any input - any input) ** 4 * 12 |
template <class Int = long long, class Int2 = long long, class Int3 = Int2>
class DelaunayTriangulation {
public:
using GPos2 = VecI2<Int, Int2>;
struct Edge {
int to;
int ccw;
int cw;
int rev;
bool on = false;
};
private:
static int isDinOABC(GPos2 a, GPos2 b, GPos2 c, GPos2 d) {
a = a - d;
b = b - d;
c = c - d;
auto val = Int3(b ^ c) * Int3(a.norm()) + Int3(c ^ a) * Int3(b.norm()) +
Int3(a ^ b) * Int3(c.norm());
return val > Int3(0) ? 1 : 0;
}
int getOpenAddress() {
if (oA.empty()) {
es.push_back({});
return (int)es.size() - 1;
}
int res = oA.back();
oA.pop_back();
return res;
}
std::pair<int, int> newEdge(int u, int v) {
int euv = getOpenAddress();
int evu = getOpenAddress();
es[euv].ccw = es[euv].cw = euv;
es[evu].ccw = es[evu].cw = evu;
es[euv].to = v;
es[evu].to = u;
es[euv].rev = evu;
es[evu].rev = euv;
es[euv].on = true;
es[evu].on = true;
return {euv, evu};
}
void eraseSingleEdge(int e) {
int ec = es[e].ccw;
int ecw = es[e].cw;
es[ec].cw = ecw;
es[ecw].ccw = ec;
es[e].on = false;
}
void eraseEdgeBidirectional(int e) {
int ex = es[e].rev;
eraseSingleEdge(e);
eraseSingleEdge(ex);
oA.push_back(e);
oA.push_back(ex);
}
void insertCcwAfter(int e, int x) {
int xc = es[x].ccw;
es[e].ccw = xc;
es[xc].cw = e;
es[e].cw = x;
es[x].ccw = e;
}
void insertCwAfter(int e, int x) {
int xcw = es[x].cw;
es[e].cw = xcw;
es[xcw].ccw = e;
es[e].ccw = x;
es[x].cw = e;
}
// move from ab to ac ... is this ccw?
int isCcw(int a, int b, int c) const {
auto ab = pos[b] - pos[a];
auto ac = pos[c] - pos[a];
auto cp = ab ^ ac;
if (0 < cp)
return 1;
if (cp < 0)
return -1;
return 0;
}
std::pair<int, int> goNext(int, int ea) {
int ap = es[ea].to;
int eap = es[es[ea].rev].ccw;
return {ap, eap};
}
std::pair<int, int> goPrev(int, int ea) {
int ap = es[es[ea].cw].to;
int eap = es[es[ea].cw].rev;
return {ap, eap};
}
std::tuple<int, int, int, int> goBottom(int a, int ea, int b, int eb) {
while (true) {
auto [ap, eap] = goPrev(a, ea);
if (isCcw(b, a, ap) > 0) {
std::tie(a, ea) = {ap, eap};
continue;
}
auto [bp, ebp] = goNext(b, eb);
if (isCcw(a, b, bp) < 0) {
std::tie(b, eb) = {bp, ebp};
continue;
}
break;
}
return {a, ea, b, eb};
}
std::pair<int, int> getMaximum(int a, int ea, bool mn) {
std::pair<int, int> ans = {a, ea};
int p = a, ep = ea;
do {
std::tie(p, ep) = goNext(p, ep);
if (mn)
ans = std::min(ans, std::make_pair(p, ep));
else
ans = std::max(ans, std::make_pair(p, ep));
} while (ep != ea);
return ans;
}
bool isDinOABC(int a, int b, int c, int d) {
return isDinOABC(pos[a], pos[b], pos[c], pos[d]);
}
std::pair<int, int> dfs(int a, int ea, int b, int eb) {
std::tie(a, ea) = getMaximum(a, ea, false);
std::tie(b, eb) = getMaximum(b, eb, true);
auto [al, eal, bl, ebl] = goBottom(a, ea, b, eb);
auto [bu, ebu, au, eau] = goBottom(b, eb, a, ea);
ebl = es[ebl].cw;
ebu = es[ebu].cw;
auto [abl, bal] = newEdge(al, bl);
insertCwAfter(abl, eal);
insertCcwAfter(bal, ebl);
if (al == au)
eau = abl;
if (bl == bu)
ebu = bal;
int ap = al, eap = eal;
int bp = bl, ebp = ebl;
while (ap != au || bp != bu) {
int a2 = es[eap].to;
int b2 = es[ebp].to;
int na = es[eap].ccw;
int nb = es[ebp].cw;
if (eap != eau && na != abl) {
int a1 = es[na].to;
if (isDinOABC(ap, bp, a2, a1)) {
eraseEdgeBidirectional(eap);
eap = na;
continue;
}
}
if (ebp != ebu && nb != bal) {
int b1 = es[nb].to;
if (isDinOABC(b2, ap, bp, b1)) {
eraseEdgeBidirectional(ebp);
ebp = nb;
continue;
}
}
bool ca = ebp == ebu;
if (eap != eau && ebp != ebu) {
if (isCcw(ap, bp, b2) < 0)
ca = true;
else if (isCcw(a2, ap, bp) < 0)
ca = false;
else
ca = isDinOABC(ap, bp, b2, a2);
}
if (ca) {
na = es[es[eap].rev].ccw;
auto [hab, hba] = newEdge(a2, bp);
insertCwAfter(hab, na);
insertCcwAfter(hba, ebp);
eap = na;
ap = a2;
} else {
nb = es[es[ebp].rev].cw;
auto [hba, hab] = newEdge(b2, ap);
insertCcwAfter(hba, nb);
insertCwAfter(hab, eap);
ebp = nb;
bp = b2;
}
}
return {al, abl};
}
std::pair<int, int> solveRange(int l, int r) {
if (r - l == 2) {
int u = l;
int v = l + 1;
auto [uv, vu] = newEdge(u, v);
return {u, uv};
}
if (r - l == 3) {
int u = l;
int v = l + 1;
int w = l + 2;
auto [uv, vu] = newEdge(u, v);
auto [vw, wv] = newEdge(v, w);
int ccw = isCcw(u, v, w);
if (ccw == 0) {
insertCcwAfter(vu, vw);
}
if (ccw > 0) {
auto [uw, wu] = newEdge(u, w);
insertCwAfter(uv, uw);
insertCwAfter(vw, vu);
insertCwAfter(wu, wv);
return {u, uv};
}
if (ccw < 0) {
auto [uw, wu] = newEdge(u, w);
insertCcwAfter(uv, uw);
insertCcwAfter(vw, vu);
insertCcwAfter(wu, wv);
return {v, vu};
}
return {u, uv};
}
int m = (l + r) / 2;
auto [a, ea] = solveRange(l, m);
auto [b, eb] = solveRange(m, r);
return dfs(a, ea, b, eb);
}
void solve() {
int sz = (int)pos.size();
if (sz <= 1)
return;
std::vector<int> pi(pos.size());
for (int i = 0; i < (int)pi.size(); i++)
pi[i] = i;
std::stable_sort(pi.begin(), pi.end(), [&](int l, int r) {
return pos[l].x != pos[r].x ? pos[l].x < pos[r].x : pos[l].y < pos[r].y;
});
auto bu1 = pos;
int ptr = 0;
mp.assign(sz, 0);
for (int i = 0; i < sz; i++) {
int v = pi[i];
if (i == 0 || !(bu1[pi[ptr - 1]] == bu1[v])) {
pi[ptr] = v;
pos[ptr++] = bu1[v];
mp[v] = v;
} else {
mp[v] = pi[ptr - 1];
}
}
if (ptr >= 2)
oOE = solveRange(0, ptr).second;
std::swap(pos, bu1);
for (auto &e : es)
e.to = pi[e.to];
}
std::vector<int> oA;
std::vector<GPos2> pos;
std::vector<Edge> es;
std::vector<int> mp;
int oOE = -1;
public:
DelaunayTriangulation() : pos() { solve(); }
DelaunayTriangulation(std::vector<GPos2> xp) : pos(std::move(xp)) { solve(); }
std::vector<std::pair<int, int>> getEdges() const {
std::vector<std::pair<int, int>> res;
for (int e = 0; e < (int)es.size(); e++)
if (es[e].on) {
int re = es[e].rev;
if (e < re)
continue;
res.push_back({es[e].to, es[re].to});
}
for (int v = 0; v < int(mp.size()); v++) {
if (mp[v] != v)
res.push_back({v, mp[v]});
}
return res;
}
};
} // namespace euclidean_mst_internal
/// @brief Return a Euclidean minimum spanning tree on integral planar points.
/// Divide-and-conquer Delaunay triangulation keeps only O(n) candidate edges:
/// the empty-circumcircle property guarantees that every Euclidean MST edge is
/// present. Kruskal on squared lengths then selects the tree; duplicate points
/// are connected by explicit zero-length edges. Coordinate differences must
/// fit in signed 64 bits; squared lengths, cross products and fourth-degree
/// incircle predicate intermediates must fit in signed 128 bits.
inline std::vector<std::pair<int, int>>
euclidean_mst(const std::vector<point<long long>> &pt) {
using internal_point = euclidean_mst_internal::VecI2<long long, __int128_t>;
using triangulation =
euclidean_mst_internal::DelaunayTriangulation<long long, __int128_t,
__int128_t>;
std::vector<internal_point> lhs;
lhs.reserve(pt.size());
for (const auto &va1 : pt) {
lhs.emplace_back(va1.x, va1.y);
}
struct weighted_edge {
std::pair<int, int> arc;
__int128_t len;
};
std::vector<weighted_edge> ord;
{
auto can = triangulation(std::move(lhs)).getEdges();
ord.reserve(can.size());
for (auto arc : can) {
__int128_t dx = __int128_t(pt[arc.first].x) - pt[arc.second].x;
__int128_t dy = __int128_t(pt[arc.first].y) - pt[arc.second].y;
ord.push_back({arc, dx * dx + dy * dy});
}
}
// Cache squared lengths so sorting never revisits point coordinates.
std::stable_sort(ord.begin(), ord.end(), [](const auto &a, const auto &b) {
return a.len < b.len;
});
atcoder::dsu ds1(int(pt.size()));
std::vector<std::pair<int, int>> ret;
ret.reserve(pt.empty() ? 0 : pt.size() - 1);
for (const auto &edge : ord) {
auto arc = edge.arc;
if (ds1.same(arc.first, arc.second)) {
continue;
}
ds1.merge(arc.first, arc.second);
if (arc.first > arc.second) {
std::swap(arc.first, arc.second);
}
ret.push_back(arc);
}
std::sort(ret.begin(), ret.end());
return ret;
}
} // namespace noya
#endif // NOYA_EUCLIDEAN_MST_HPP
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O(n log n).
/// Space: O(n).
#include <algorithm>
#include <cassert>
#include <vector>
namespace atcoder {
// Implement (union by size) + (path compression)
// Reference:
// Zvi Galil and Giuseppe F. Italiano,
// Data structures and algorithms for disjoint set union problems
struct dsu {
public:
dsu() : _n(0) {}
explicit dsu(int n) : _n(n), parent_or_size(n, -1) {}
int merge(int a, int b) {
assert(0 <= a && a < _n);
assert(0 <= b && b < _n);
int x = leader(a), y = leader(b);
if (x == y) return x;
if (-parent_or_size[x] < -parent_or_size[y]) std::swap(x, y);
parent_or_size[x] += parent_or_size[y];
parent_or_size[y] = x;
return x;
}
bool same(int a, int b) {
assert(0 <= a && a < _n);
assert(0 <= b && b < _n);
return leader(a) == leader(b);
}
int leader(int a) {
assert(0 <= a && a < _n);
return _leader(a);
}
int size(int a) {
assert(0 <= a && a < _n);
return -parent_or_size[leader(a)];
}
std::vector<std::vector<int>> groups() {
std::vector<int> leader_buf(_n), group_size(_n);
for (int i = 0; i < _n; i++) {
leader_buf[i] = leader(i);
group_size[leader_buf[i]]++;
}
std::vector<std::vector<int>> result(_n);
for (int i = 0; i < _n; i++) {
result[i].reserve(group_size[i]);
}
for (int i = 0; i < _n; i++) {
result[leader_buf[i]].push_back(i);
}
result.erase(
std::remove_if(result.begin(), result.end(),
[&](const std::vector<int>& v) { return v.empty(); }),
result.end());
return result;
}
private:
int _n;
// root node: -1 * component size
// otherwise: parent
std::vector<int> parent_or_size;
int _leader(int a) {
if (parent_or_size[a] < 0) return a;
return parent_or_size[a] = _leader(parent_or_size[a]);
}
};
} // namespace atcoder
/// @complexity Time: O(1) per point/segment primitive, O(n) for polygon area or
/// point location, O(n log n) for convex hull. Space: O(1), except O(n) for hull.
#include <algorithm>
#include <cmath>
#include <optional>
#include <vector>
namespace noya {
/// @brief Two-dimensional point with vector arithmetic and lexicographic order.
template <class T> struct point {
T x{};
T y{};
point() = default;
point(T x_, T y_) : x(x_), y(y_) {}
point &operator+=(const point &rhs) {
x += rhs.x;
y += rhs.y;
return *this;
}
point &operator-=(const point &rhs) {
x -= rhs.x;
y -= rhs.y;
return *this;
}
point &operator*=(const T &scl) {
x *= scl;
y *= scl;
return *this;
}
point &operator/=(const T &scl) {
x /= scl;
y /= scl;
return *this;
}
friend point operator+(point l, const point &r) { return l += r; }
friend point operator-(point l, const point &r) { return l -= r; }
friend point operator*(point val, const T &scl) { return val *= scl; }
friend point operator*(const T &scl, point val) { return val *= scl; }
friend point operator/(point val, const T &scl) { return val /= scl; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &l, const point &r) {
return l.x < r.x || (l.x == r.x && l.y < r.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> o;
Real r1{};
};
/// @brief Return the dot product of two vectors.
template <class T> T dot(const point<T> &a, const point<T> &b) {
return a.x * b.x + a.y * b.y;
}
/// @brief Return the signed cross product of two vectors.
template <class T> T cross(const point<T> &a, const point<T> &b) {
return a.x * b.y - a.y * b.x;
}
/// @brief Return cross(a - o1, b - o1).
template <class T>
T cross(const point<T> &o1, const point<T> &a, const point<T> &b) {
return cross(a - o1, b - o1);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &val) { return dot(val, val); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &val, const T &eps = T{}) {
return (val > eps) - (val < -eps);
}
/// @brief Return -1, 0, or 1 for a clockwise, collinear, or counter-clockwise
/// turn.
template <class T>
int orientation(const point<T> &a, const point<T> &b, const point<T> &c,
const T &eps = T{}) {
return sign(cross(a, b, c), eps);
}
/// @brief Test whether p lies on the closed segment [a, b].
template <class T>
bool on_segment(const point<T> &p, const point<T> &a, const point<T> &b,
const T &eps = T{}) {
if (orientation(a, b, p, eps) != 0) {
return false;
}
return std::min(a.x, b.x) - eps <= p.x && p.x <= std::max(a.x, b.x) + eps &&
std::min(a.y, b.y) - eps <= p.y && p.y <= std::max(a.y, b.y) + eps;
}
/// @brief Test whether the closed segments [a, b] and [c, d] intersect.
template <class T>
bool segments_intersect(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, const T &eps = T{}) {
int s1 = orientation(a, b, c, eps);
int s2 = orientation(a, b, d, eps);
int s3 = orientation(c, d, a, eps);
int s4 = orientation(c, d, b, eps);
if (s1 == 0 && on_segment(c, a, b, eps)) {
return true;
}
if (s2 == 0 && on_segment(d, a, b, eps)) {
return true;
}
if (s3 == 0 && on_segment(a, c, d, eps)) {
return true;
}
if (s4 == 0 && on_segment(b, c, d, eps)) {
return true;
}
return s1 * s2 < 0 && s3 * s4 < 0;
}
/// @brief Intersect the infinite lines through (a, b) and (c, d), returning
/// nullopt when they are parallel or coincident.
template <class T>
std::optional<point<long double>>
line_intersection(const point<T> &a, const point<T> &b, const point<T> &c,
const point<T> &d, long double eps = 0) {
point<long double> lhs{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> b1{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> z{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> d1{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> da = b1 - lhs;
point<long double> db = d1 - z;
long double den = cross(da, db);
if (std::abs(den) <= eps) {
return std::nullopt;
}
long double rat = cross(z - lhs, db) / den;
return lhs + da * rat;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
/// Starts at the lexicographically smallest point; kc retains boundary points.
/// Coordinate differences and determinant arithmetic must fit in T.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> pt, bool kc = false) {
std::sort(pt.begin(), pt.end());
pt.erase(std::unique(pt.begin(), pt.end()), pt.end());
if (pt.size() <= 1) {
return pt;
}
if (kc) {
bool col = true;
for (int i = 2; i < int(pt.size()) && col; i++) {
col = orientation(pt[0], pt[1], pt[i]) == 0;
}
if (col) return pt;
}
std::vector<point<T>> res;
res.reserve(pt.size() + 1);
for (const point<T> &p : pt) {
while (res.size() >= 2) {
int rot = orientation(res[res.size() - 2], res.back(), p);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(p);
}
auto k = res.size(); // Keep the lower chain while building the upper chain.
for (int i = int(pt.size()) - 2; i >= 0; i--) {
while (res.size() > k) {
int rot = orientation(res[res.size() - 2], res.back(), pt[i]);
if (rot > 0 || (kc && rot == 0)) {
break;
}
res.pop_back();
}
res.push_back(pt[i]);
}
res.pop_back();
return res;
}
/// @brief Return twice the signed area of a polygon.
/// Coordinate differences, products and the accumulated area must fit in T.
template <class T> T polygon_area2(const std::vector<point<T>> &pg) {
T res{};
if (pg.size() < 3) return res;
auto prev = pg[1] - pg[0];
for (int i = 2; i < int(pg.size()); i++) {
auto cur = pg[i] - pg[0];
res += cross(prev, cur);
prev = cur;
}
return res;
}
/// @brief Classify a point relative to a polygon: -1 outside, 0 boundary, 1
/// inside.
template <class T>
int point_in_polygon(const point<T> &p, const std::vector<point<T>> &pg) {
bool in = false;
for (int i = 0; i < int(pg.size()); i++) {
point<T> a = pg[i];
point<T> b = pg[(i + 1) % pg.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
in = !in;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
in = !in;
}
}
return in ? 1 : -1;
}
/// @brief Return the Euclidean distance between two points.
template <class T> long double distance(const point<T> &a, const point<T> &b) {
return std::hypot(
static_cast<long double>(a.x) - static_cast<long double>(b.x),
static_cast<long double>(a.y) - static_cast<long double>(b.y));
}
} // namespace noya
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
namespace euclidean_mst_internal {
template <class Int = long long, class Int2 = long long> struct VecI2 {
Int x, y;
VecI2() : x(0), y(0) {}
VecI2(std::pair<Int, Int> _p)
: x(std::move(_p.first)), y(std::move(_p.second)) {}
VecI2(Int _x, Int _y) : x(std::move(_x)), y(std::move(_y)) {}
VecI2 &operator+=(VecI2 r) {
x += r.x;
y += r.y;
return *this;
}
VecI2 &operator-=(VecI2 r) {
x -= r.x;
y -= r.y;
return *this;
}
VecI2 &operator*=(Int r) {
x *= r;
y *= r;
return *this;
}
VecI2 operator+(VecI2 r) const { return VecI2(x + r.x, y + r.y); }
VecI2 operator-(VecI2 r) const { return VecI2(x - r.x, y - r.y); }
VecI2 operator*(Int r) const { return VecI2(x * r, y * r); }
VecI2 operator-() const { return VecI2(-x, -y); }
Int2 operator*(VecI2 r) const {
return Int2(x) * Int2(r.x) + Int2(y) * Int2(r.y);
}
Int2 operator^(VecI2 r) const {
return Int2(x) * Int2(r.y) - Int2(y) * Int2(r.x);
}
bool operator<(VecI2 r) const { return x < r.x || (!(r.x < x) && y < r.y); }
Int2 norm() const { return Int2(x) * Int2(x) + Int2(y) * Int2(y); }
static bool compareYX(VecI2 a, VecI2 b) {
return a.y < b.y || (!(b.y < a.y) && a.x < b.x);
}
static bool compareXY(VecI2 a, VecI2 b) {
return a.x < b.x || (!(b.x < a.x) && a.y < b.y);
}
bool operator==(VecI2 r) const { return x == r.x && y == r.y; }
bool operator!=(VecI2 r) const { return x != r.x || y != r.y; }
};
// Int3 must be able to handle the value range :
// |x| <= | (any input - any input) ** 4 * 12 |
template <class Int = long long, class Int2 = long long, class Int3 = Int2>
class DelaunayTriangulation {
public:
using GPos2 = VecI2<Int, Int2>;
struct Edge {
int to;
int ccw;
int cw;
int rev;
bool on = false;
};
private:
static int isDinOABC(GPos2 a, GPos2 b, GPos2 c, GPos2 d) {
a = a - d;
b = b - d;
c = c - d;
auto val = Int3(b ^ c) * Int3(a.norm()) + Int3(c ^ a) * Int3(b.norm()) +
Int3(a ^ b) * Int3(c.norm());
return val > Int3(0) ? 1 : 0;
}
int getOpenAddress() {
if (oA.empty()) {
es.push_back({});
return (int)es.size() - 1;
}
int res = oA.back();
oA.pop_back();
return res;
}
std::pair<int, int> newEdge(int u, int v) {
int euv = getOpenAddress();
int evu = getOpenAddress();
es[euv].ccw = es[euv].cw = euv;
es[evu].ccw = es[evu].cw = evu;
es[euv].to = v;
es[evu].to = u;
es[euv].rev = evu;
es[evu].rev = euv;
es[euv].on = true;
es[evu].on = true;
return {euv, evu};
}
void eraseSingleEdge(int e) {
int ec = es[e].ccw;
int ecw = es[e].cw;
es[ec].cw = ecw;
es[ecw].ccw = ec;
es[e].on = false;
}
void eraseEdgeBidirectional(int e) {
int ex = es[e].rev;
eraseSingleEdge(e);
eraseSingleEdge(ex);
oA.push_back(e);
oA.push_back(ex);
}
void insertCcwAfter(int e, int x) {
int xc = es[x].ccw;
es[e].ccw = xc;
es[xc].cw = e;
es[e].cw = x;
es[x].ccw = e;
}
void insertCwAfter(int e, int x) {
int xcw = es[x].cw;
es[e].cw = xcw;
es[xcw].ccw = e;
es[e].ccw = x;
es[x].cw = e;
}
// move from ab to ac ... is this ccw?
int isCcw(int a, int b, int c) const {
auto ab = pos[b] - pos[a];
auto ac = pos[c] - pos[a];
auto cp = ab ^ ac;
if (0 < cp)
return 1;
if (cp < 0)
return -1;
return 0;
}
std::pair<int, int> goNext(int, int ea) {
int ap = es[ea].to;
int eap = es[es[ea].rev].ccw;
return {ap, eap};
}
std::pair<int, int> goPrev(int, int ea) {
int ap = es[es[ea].cw].to;
int eap = es[es[ea].cw].rev;
return {ap, eap};
}
std::tuple<int, int, int, int> goBottom(int a, int ea, int b, int eb) {
while (true) {
auto [ap, eap] = goPrev(a, ea);
if (isCcw(b, a, ap) > 0) {
std::tie(a, ea) = {ap, eap};
continue;
}
auto [bp, ebp] = goNext(b, eb);
if (isCcw(a, b, bp) < 0) {
std::tie(b, eb) = {bp, ebp};
continue;
}
break;
}
return {a, ea, b, eb};
}
std::pair<int, int> getMaximum(int a, int ea, bool mn) {
std::pair<int, int> ans = {a, ea};
int p = a, ep = ea;
do {
std::tie(p, ep) = goNext(p, ep);
if (mn)
ans = std::min(ans, std::make_pair(p, ep));
else
ans = std::max(ans, std::make_pair(p, ep));
} while (ep != ea);
return ans;
}
bool isDinOABC(int a, int b, int c, int d) {
return isDinOABC(pos[a], pos[b], pos[c], pos[d]);
}
std::pair<int, int> dfs(int a, int ea, int b, int eb) {
std::tie(a, ea) = getMaximum(a, ea, false);
std::tie(b, eb) = getMaximum(b, eb, true);
auto [al, eal, bl, ebl] = goBottom(a, ea, b, eb);
auto [bu, ebu, au, eau] = goBottom(b, eb, a, ea);
ebl = es[ebl].cw;
ebu = es[ebu].cw;
auto [abl, bal] = newEdge(al, bl);
insertCwAfter(abl, eal);
insertCcwAfter(bal, ebl);
if (al == au)
eau = abl;
if (bl == bu)
ebu = bal;
int ap = al, eap = eal;
int bp = bl, ebp = ebl;
while (ap != au || bp != bu) {
int a2 = es[eap].to;
int b2 = es[ebp].to;
int na = es[eap].ccw;
int nb = es[ebp].cw;
if (eap != eau && na != abl) {
int a1 = es[na].to;
if (isDinOABC(ap, bp, a2, a1)) {
eraseEdgeBidirectional(eap);
eap = na;
continue;
}
}
if (ebp != ebu && nb != bal) {
int b1 = es[nb].to;
if (isDinOABC(b2, ap, bp, b1)) {
eraseEdgeBidirectional(ebp);
ebp = nb;
continue;
}
}
bool ca = ebp == ebu;
if (eap != eau && ebp != ebu) {
if (isCcw(ap, bp, b2) < 0)
ca = true;
else if (isCcw(a2, ap, bp) < 0)
ca = false;
else
ca = isDinOABC(ap, bp, b2, a2);
}
if (ca) {
na = es[es[eap].rev].ccw;
auto [hab, hba] = newEdge(a2, bp);
insertCwAfter(hab, na);
insertCcwAfter(hba, ebp);
eap = na;
ap = a2;
} else {
nb = es[es[ebp].rev].cw;
auto [hba, hab] = newEdge(b2, ap);
insertCcwAfter(hba, nb);
insertCwAfter(hab, eap);
ebp = nb;
bp = b2;
}
}
return {al, abl};
}
std::pair<int, int> solveRange(int l, int r) {
if (r - l == 2) {
int u = l;
int v = l + 1;
auto [uv, vu] = newEdge(u, v);
return {u, uv};
}
if (r - l == 3) {
int u = l;
int v = l + 1;
int w = l + 2;
auto [uv, vu] = newEdge(u, v);
auto [vw, wv] = newEdge(v, w);
int ccw = isCcw(u, v, w);
if (ccw == 0) {
insertCcwAfter(vu, vw);
}
if (ccw > 0) {
auto [uw, wu] = newEdge(u, w);
insertCwAfter(uv, uw);
insertCwAfter(vw, vu);
insertCwAfter(wu, wv);
return {u, uv};
}
if (ccw < 0) {
auto [uw, wu] = newEdge(u, w);
insertCcwAfter(uv, uw);
insertCcwAfter(vw, vu);
insertCcwAfter(wu, wv);
return {v, vu};
}
return {u, uv};
}
int m = (l + r) / 2;
auto [a, ea] = solveRange(l, m);
auto [b, eb] = solveRange(m, r);
return dfs(a, ea, b, eb);
}
void solve() {
int sz = (int)pos.size();
if (sz <= 1)
return;
std::vector<int> pi(pos.size());
for (int i = 0; i < (int)pi.size(); i++)
pi[i] = i;
std::stable_sort(pi.begin(), pi.end(), [&](int l, int r) {
return pos[l].x != pos[r].x ? pos[l].x < pos[r].x : pos[l].y < pos[r].y;
});
auto bu1 = pos;
int ptr = 0;
mp.assign(sz, 0);
for (int i = 0; i < sz; i++) {
int v = pi[i];
if (i == 0 || !(bu1[pi[ptr - 1]] == bu1[v])) {
pi[ptr] = v;
pos[ptr++] = bu1[v];
mp[v] = v;
} else {
mp[v] = pi[ptr - 1];
}
}
if (ptr >= 2)
oOE = solveRange(0, ptr).second;
std::swap(pos, bu1);
for (auto &e : es)
e.to = pi[e.to];
}
std::vector<int> oA;
std::vector<GPos2> pos;
std::vector<Edge> es;
std::vector<int> mp;
int oOE = -1;
public:
DelaunayTriangulation() : pos() { solve(); }
DelaunayTriangulation(std::vector<GPos2> xp) : pos(std::move(xp)) { solve(); }
std::vector<std::pair<int, int>> getEdges() const {
std::vector<std::pair<int, int>> res;
for (int e = 0; e < (int)es.size(); e++)
if (es[e].on) {
int re = es[e].rev;
if (e < re)
continue;
res.push_back({es[e].to, es[re].to});
}
for (int v = 0; v < int(mp.size()); v++) {
if (mp[v] != v)
res.push_back({v, mp[v]});
}
return res;
}
};
} // namespace euclidean_mst_internal
/// @brief Return a Euclidean minimum spanning tree on integral planar points.
/// Divide-and-conquer Delaunay triangulation keeps only O(n) candidate edges:
/// the empty-circumcircle property guarantees that every Euclidean MST edge is
/// present. Kruskal on squared lengths then selects the tree; duplicate points
/// are connected by explicit zero-length edges. Coordinate differences must
/// fit in signed 64 bits; squared lengths, cross products and fourth-degree
/// incircle predicate intermediates must fit in signed 128 bits.
inline std::vector<std::pair<int, int>>
euclidean_mst(const std::vector<point<long long>> &pt) {
using internal_point = euclidean_mst_internal::VecI2<long long, __int128_t>;
using triangulation =
euclidean_mst_internal::DelaunayTriangulation<long long, __int128_t,
__int128_t>;
std::vector<internal_point> lhs;
lhs.reserve(pt.size());
for (const auto &va1 : pt) {
lhs.emplace_back(va1.x, va1.y);
}
struct weighted_edge {
std::pair<int, int> arc;
__int128_t len;
};
std::vector<weighted_edge> ord;
{
auto can = triangulation(std::move(lhs)).getEdges();
ord.reserve(can.size());
for (auto arc : can) {
__int128_t dx = __int128_t(pt[arc.first].x) - pt[arc.second].x;
__int128_t dy = __int128_t(pt[arc.first].y) - pt[arc.second].y;
ord.push_back({arc, dx * dx + dy * dy});
}
}
// Cache squared lengths so sorting never revisits point coordinates.
std::stable_sort(ord.begin(), ord.end(), [](const auto &a, const auto &b) {
return a.len < b.len;
});
atcoder::dsu ds1(int(pt.size()));
std::vector<std::pair<int, int>> ret;
ret.reserve(pt.empty() ? 0 : pt.size() - 1);
for (const auto &edge : ord) {
auto arc = edge.arc;
if (ds1.same(arc.first, arc.second)) {
continue;
}
ds1.merge(arc.first, arc.second);
if (arc.first > arc.second) {
std::swap(arc.first, arc.second);
}
ret.push_back(arc);
}
std::sort(ret.begin(), ret.end());
return ret;
}
} // namespace noya
// NOYA_EUCLIDEAN_MST_HPP