description: 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.¶
euclidean_mst.hpp¶
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.
Verified by euclidean_mst.
求整数平面点的欧氏最小生成树；通过 Delaunay 候选边避免完全图。
Implementation¶
#ifndef NOYA_EUCLIDEAN_MST_HPP
#define NOYA_EUCLIDEAN_MST_HPP 1
/// @complexity Time: O(n log n).
/// Space: O(n).
#include "atcoder/dsu.hpp"
#include "noya/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; }
};
template<class Elem>
class CsrArray{
public:
struct ListRange{
using iterator = typename std::vector<Elem>::iterator;
iterator begi, endi;
iterator begin() const { return begi; }
iterator end() const { return endi; }
int size() const { return (int)std::distance(begi, endi); }
Elem& operator[](int i) const { return begi[i]; }
};
struct ConstListRange{
using iterator = typename std::vector<Elem>::const_iterator;
iterator begi, endi;
iterator begin() const { return begi; }
iterator end() const { return endi; }
int size() const { return (int)std::distance(begi, endi); }
const Elem& operator[](int i) const { return begi[i]; }
};
private:
int m_n;
std::vector<Elem> m_list;
std::vector<int> m_pos;
public:
CsrArray() : m_n(0), m_list(), m_pos() {}
static CsrArray Construct(int n, std::vector<std::pair<int, Elem>> items){
CsrArray res;
res.m_n = n;
std::vector<int> buf(n+1, 0);
for(auto& [u,v] : items){ ++buf[u]; }
for(int i=1; i<=n; i++) buf[i] += buf[i-1];
res.m_list.resize(buf[n]);
for(int i=(int)items.size()-1; i>=0; i--){
res.m_list[--buf[items[i].first]] = std::move(items[i].second);
}
res.m_pos = std::move(buf);
return res;
}
static CsrArray FromRaw(std::vector<Elem> list, std::vector<int> pos){
CsrArray res;
res.m_n = pos.size() - 1;
res.m_list = std::move(list);
res.m_pos = std::move(pos);
return res;
}
ListRange operator[](int u) { return ListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
ConstListRange operator[](int u) const { return ConstListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
int size() const { return m_n; }
int fullSize() const { return (int)m_list.size(); }
};
// 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 enabled = 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(openAddress.empty()){
edges.push_back({});
return (int)edges.size() - 1;
}
int res = openAddress.back();
openAddress.pop_back();
return res;
}
std::pair<int, int> newEdge(int u, int v){
int euv = getOpenAddress();
int evu = getOpenAddress();
edges[euv].ccw = edges[euv].cw = euv;
edges[evu].ccw = edges[evu].cw = evu;
edges[euv].to = v;
edges[evu].to = u;
edges[euv].rev = evu;
edges[evu].rev = euv;
edges[euv].enabled = true;
edges[evu].enabled = true;
return { euv, evu };
}
void eraseSingleEdge(int e){
int eccw = edges[e].ccw;
int ecw = edges[e].cw;
edges[eccw].cw = ecw;
edges[ecw].ccw = eccw;
edges[e].enabled = false;
}
void eraseEdgeBidirectional(int e){
int ex = edges[e].rev;
eraseSingleEdge(e);
eraseSingleEdge(ex);
openAddress.push_back(e);
openAddress.push_back(ex);
}
void insertCcwAfter(int e, int x){
int xccw = edges[x].ccw;
edges[e].ccw = xccw;
edges[xccw].cw = e;
edges[e].cw = x;
edges[x].ccw = e;
}
void insertCwAfter(int e, int x){
int xcw = edges[x].cw;
edges[e].cw = xcw;
edges[xcw].ccw = e;
edges[e].ccw = x;
edges[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 = edges[ea].to;
int eap = edges[edges[ea].rev].ccw;
return { ap, eap };
}
std::pair<int, int> goPrev(int , int ea){
int ap = edges[edges[ea].cw].to;
int eap = edges[edges[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 toMin){
std::pair<int, int> ans = { a, ea };
int p = a, ep = ea;
do {
std::tie(p, ep) = goNext(p, ep);
if(toMin) 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 = edges[ebl].cw;
ebu = edges[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 = edges[eap].to;
int b2 = edges[ebp].to;
int nxeap = edges[eap].ccw;
int nxebp = edges[ebp].cw;
if(eap != eau && nxeap != abl){
int a1 = edges[nxeap].to;
if(isDinOABC(ap, bp, a2, a1)){
eraseEdgeBidirectional(eap);
eap = nxeap;
continue;
}
}
if(ebp != ebu && nxebp != bal){
int b1 = edges[nxebp].to;
if(isDinOABC(b2, ap, bp, b1)){
eraseEdgeBidirectional(ebp);
ebp = nxebp;
continue;
}
}
bool chooseA = ebp == ebu;
if(eap != eau && ebp != ebu){
if(isCcw(ap, bp, b2) < 0) chooseA = true;
else if(isCcw(a2, ap, bp) < 0) chooseA = false;
else chooseA = isDinOABC(ap, bp, b2, a2);
}
if(chooseA){
nxeap = edges[edges[eap].rev].ccw;
auto [hab, hba] = newEdge(a2, bp);
insertCwAfter(hab, nxeap);
insertCcwAfter(hba, ebp);
eap = nxeap; ap = a2;
}
else {
nxebp = edges[edges[ebp].rev].cw;
auto [hba, hab] = newEdge(b2, ap);
insertCcwAfter(hba, nxebp);
insertCwAfter(hab, eap);
ebp = nxebp; 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 posbuf = pos;
int posptr = 0;
mappings.assign(sz, 0);
for(int i=0; i<sz; i++){
int v = pi[i];
if(i == 0 || !(posbuf[pi[posptr-1]] == posbuf[v])){
pi[posptr] = v;
pos[posptr++] = posbuf[v];
mappings[v] = v;
} else {
mappings[v] = pi[posptr-1];
}
}
if(posptr >= 2) outerOneEdge = solveRange(0, posptr).second;
std::swap(pos, posbuf);
for(auto& e : edges) e.to = pi[e.to];
}
std::vector<int> openAddress;
std::vector<GPos2> pos;
std::vector<Edge> edges;
std::vector<int> mappings;
int outerOneEdge = -1;
public:
DelaunayTriangulation()
: pos()
{ solve(); }
DelaunayTriangulation(std::vector<GPos2> x_points)
: pos(std::move(x_points))
{
solve();
}
std::vector<std::pair<int, int>> getEdges() const {
std::vector<std::pair<int, int>> res;
for(int e=0; e<(int)edges.size(); e++) if(edges[e].enabled){
int re = edges[e].rev;
if(e < re) continue;
res.push_back({ edges[e].to, edges[re].to });
}
for(int v=0; v<int(mappings.size()); v++){
if(mappings[v] != v) res.push_back({ v, mappings[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.
inline std::vector<std::pair<int, int>>
euclidean_mst(const std::vector<point<long long>> &points) {
using internal_point =
euclidean_mst_internal::VecI2<long long, long long>;
using triangulation = euclidean_mst_internal::DelaunayTriangulation<
long long, long long, __int128_t>;
std::vector<internal_point> converted;
converted.reserve(points.size());
for (const auto &value : points) {
converted.emplace_back(value.x, value.y);
}
std::vector<std::pair<int, int>> candidates =
triangulation(std::move(converted)).getEdges();
auto squared_distance = [&](const std::pair<int, int> &edge) {
__int128_t dx = __int128_t(points[edge.first].x) - points[edge.second].x;
__int128_t dy = __int128_t(points[edge.first].y) - points[edge.second].y;
return dx * dx + dy * dy;
};
std::stable_sort(candidates.begin(), candidates.end(),
[&](const auto &first, const auto &second) {
return squared_distance(first) < squared_distance(second);
});
atcoder::dsu components(int(points.size()));
std::vector<std::pair<int, int>> result;
result.reserve(points.empty() ? 0 : points.size() - 1);
for (auto edge : candidates) {
if (components.same(edge.first, edge.second)) {
continue;
}
components.merge(edge.first, edge.second);
if (edge.first > edge.second) {
std::swap(edge.first, edge.second);
}
result.push_back(edge);
}
std::sort(result.begin(), result.end());
return result;
}
} // namespace noya
#endif // NOYA_EUCLIDEAN_MST_HPP
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <numeric>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>
/// @complexity Time: O(n log n).
/// Space: O(n).
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 primitive; O(n log n) for convex hull.
/// Space: O(1) per primitive and O(n) for hull construction.
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 &other) {
x += other.x;
y += other.y;
return *this;
}
point &operator-=(const point &other) {
x -= other.x;
y -= other.y;
return *this;
}
point &operator*=(const T &scale) {
x *= scale;
y *= scale;
return *this;
}
point &operator/=(const T &scale) {
x /= scale;
y /= scale;
return *this;
}
friend point operator+(point left, const point &right) {
return left += right;
}
friend point operator-(point left, const point &right) {
return left -= right;
}
friend point operator*(point value, const T &scale) { return value *= scale; }
friend point operator*(const T &scale, point value) { return value *= scale; }
friend point operator/(point value, const T &scale) { return value /= scale; }
friend bool operator==(const point &, const point &) = default;
friend bool operator<(const point &left, const point &right) {
return left.x < right.x || (left.x == right.x && left.y < right.y);
}
};
/// @brief Circle represented by a center and a nonnegative radius.
template <class Real> struct circle {
point<Real> center;
Real radius{};
};
/// @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 - origin, b - origin).
template <class T>
T cross(const point<T> &origin, const point<T> &a, const point<T> &b) {
return cross(a - origin, b - origin);
}
/// @brief Return the squared Euclidean norm.
template <class T> T norm2(const point<T> &value) { return dot(value, value); }
/// @brief Compare a value with zero using an optional absolute tolerance.
template <class T> int sign(const T &value, const T &epsilon = T{}) {
return (value > epsilon) - (value < -epsilon);
}
/// @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 &epsilon = T{}) {
return sign(cross(a, b, c), epsilon);
}
/// @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 &epsilon = T{}) {
if (orientation(a, b, p, epsilon) != 0) {
return false;
}
return std::min(a.x, b.x) - epsilon <= p.x &&
p.x <= std::max(a.x, b.x) + epsilon &&
std::min(a.y, b.y) - epsilon <= p.y &&
p.y <= std::max(a.y, b.y) + epsilon;
}
/// @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 &epsilon = T{}) {
int ab_c = orientation(a, b, c, epsilon);
int ab_d = orientation(a, b, d, epsilon);
int cd_a = orientation(c, d, a, epsilon);
int cd_b = orientation(c, d, b, epsilon);
if (ab_c == 0 && on_segment(c, a, b, epsilon)) {
return true;
}
if (ab_d == 0 && on_segment(d, a, b, epsilon)) {
return true;
}
if (cd_a == 0 && on_segment(a, c, d, epsilon)) {
return true;
}
if (cd_b == 0 && on_segment(b, c, d, epsilon)) {
return true;
}
return ab_c * ab_d < 0 && cd_a * cd_b < 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 epsilon = 0) {
point<long double> first{static_cast<long double>(a.x),
static_cast<long double>(a.y)};
point<long double> second{static_cast<long double>(b.x),
static_cast<long double>(b.y)};
point<long double> third{static_cast<long double>(c.x),
static_cast<long double>(c.y)};
point<long double> fourth{static_cast<long double>(d.x),
static_cast<long double>(d.y)};
point<long double> direction_a = second - first;
point<long double> direction_b = fourth - third;
long double denominator = cross(direction_a, direction_b);
if (std::abs(denominator) <= epsilon) {
return std::nullopt;
}
long double ratio = cross(third - first, direction_b) / denominator;
return first + direction_a * ratio;
}
/// @brief Return the convex hull in counter-clockwise order without repetition.
template <class T>
std::vector<point<T>> convex_hull(std::vector<point<T>> points,
bool keep_collinear = false) {
std::sort(points.begin(), points.end());
points.erase(std::unique(points.begin(), points.end()), points.end());
if (points.size() <= 1) {
return points;
}
bool all_collinear = true;
for (int i = 2; i < int(points.size()); i++) {
all_collinear &= orientation(points[0], points[1], points[i]) == 0;
}
if (keep_collinear && all_collinear) {
return points;
}
std::vector<point<T>> lower, upper;
for (const point<T> &p : points) {
while (lower.size() >= 2) {
int turn = orientation(lower[lower.size() - 2], lower.back(), p);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
lower.pop_back();
}
lower.push_back(p);
}
for (auto it = points.rbegin(); it != points.rend(); ++it) {
while (upper.size() >= 2) {
int turn = orientation(upper[upper.size() - 2], upper.back(), *it);
if (turn > 0 || (keep_collinear && turn == 0)) {
break;
}
upper.pop_back();
}
upper.push_back(*it);
}
lower.pop_back();
upper.pop_back();
lower.insert(lower.end(), upper.begin(), upper.end());
return lower;
}
/// @brief Return twice the signed area of a polygon.
template <class T> T polygon_area2(const std::vector<point<T>> &polygon) {
T result{};
for (int i = 0; i < int(polygon.size()); i++) {
result += cross(polygon[i], polygon[(i + 1) % polygon.size()]);
}
return result;
}
/// @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>> &polygon) {
bool inside = false;
for (int i = 0; i < int(polygon.size()); i++) {
point<T> a = polygon[i];
point<T> b = polygon[(i + 1) % polygon.size()];
if (on_segment(p, a, b)) {
return 0;
}
if (a.y <= p.y && p.y < b.y && orientation(a, b, p) > 0) {
inside = !inside;
}
if (b.y <= p.y && p.y < a.y && orientation(a, b, p) < 0) {
inside = !inside;
}
}
return inside ? 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
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; }
};
template<class Elem>
class CsrArray{
public:
struct ListRange{
using iterator = typename std::vector<Elem>::iterator;
iterator begi, endi;
iterator begin() const { return begi; }
iterator end() const { return endi; }
int size() const { return (int)std::distance(begi, endi); }
Elem& operator[](int i) const { return begi[i]; }
};
struct ConstListRange{
using iterator = typename std::vector<Elem>::const_iterator;
iterator begi, endi;
iterator begin() const { return begi; }
iterator end() const { return endi; }
int size() const { return (int)std::distance(begi, endi); }
const Elem& operator[](int i) const { return begi[i]; }
};
private:
int m_n;
std::vector<Elem> m_list;
std::vector<int> m_pos;
public:
CsrArray() : m_n(0), m_list(), m_pos() {}
static CsrArray Construct(int n, std::vector<std::pair<int, Elem>> items){
CsrArray res;
res.m_n = n;
std::vector<int> buf(n+1, 0);
for(auto& [u,v] : items){ ++buf[u]; }
for(int i=1; i<=n; i++) buf[i] += buf[i-1];
res.m_list.resize(buf[n]);
for(int i=(int)items.size()-1; i>=0; i--){
res.m_list[--buf[items[i].first]] = std::move(items[i].second);
}
res.m_pos = std::move(buf);
return res;
}
static CsrArray FromRaw(std::vector<Elem> list, std::vector<int> pos){
CsrArray res;
res.m_n = pos.size() - 1;
res.m_list = std::move(list);
res.m_pos = std::move(pos);
return res;
}
ListRange operator[](int u) { return ListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
ConstListRange operator[](int u) const { return ConstListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
int size() const { return m_n; }
int fullSize() const { return (int)m_list.size(); }
};
// 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 enabled = 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(openAddress.empty()){
edges.push_back({});
return (int)edges.size() - 1;
}
int res = openAddress.back();
openAddress.pop_back();
return res;
}
std::pair<int, int> newEdge(int u, int v){
int euv = getOpenAddress();
int evu = getOpenAddress();
edges[euv].ccw = edges[euv].cw = euv;
edges[evu].ccw = edges[evu].cw = evu;
edges[euv].to = v;
edges[evu].to = u;
edges[euv].rev = evu;
edges[evu].rev = euv;
edges[euv].enabled = true;
edges[evu].enabled = true;
return { euv, evu };
}
void eraseSingleEdge(int e){
int eccw = edges[e].ccw;
int ecw = edges[e].cw;
edges[eccw].cw = ecw;
edges[ecw].ccw = eccw;
edges[e].enabled = false;
}
void eraseEdgeBidirectional(int e){
int ex = edges[e].rev;
eraseSingleEdge(e);
eraseSingleEdge(ex);
openAddress.push_back(e);
openAddress.push_back(ex);
}
void insertCcwAfter(int e, int x){
int xccw = edges[x].ccw;
edges[e].ccw = xccw;
edges[xccw].cw = e;
edges[e].cw = x;
edges[x].ccw = e;
}
void insertCwAfter(int e, int x){
int xcw = edges[x].cw;
edges[e].cw = xcw;
edges[xcw].ccw = e;
edges[e].ccw = x;
edges[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 = edges[ea].to;
int eap = edges[edges[ea].rev].ccw;
return { ap, eap };
}
std::pair<int, int> goPrev(int , int ea){
int ap = edges[edges[ea].cw].to;
int eap = edges[edges[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 toMin){
std::pair<int, int> ans = { a, ea };
int p = a, ep = ea;
do {
std::tie(p, ep) = goNext(p, ep);
if(toMin) 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 = edges[ebl].cw;
ebu = edges[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 = edges[eap].to;
int b2 = edges[ebp].to;
int nxeap = edges[eap].ccw;
int nxebp = edges[ebp].cw;
if(eap != eau && nxeap != abl){
int a1 = edges[nxeap].to;
if(isDinOABC(ap, bp, a2, a1)){
eraseEdgeBidirectional(eap);
eap = nxeap;
continue;
}
}
if(ebp != ebu && nxebp != bal){
int b1 = edges[nxebp].to;
if(isDinOABC(b2, ap, bp, b1)){
eraseEdgeBidirectional(ebp);
ebp = nxebp;
continue;
}
}
bool chooseA = ebp == ebu;
if(eap != eau && ebp != ebu){
if(isCcw(ap, bp, b2) < 0) chooseA = true;
else if(isCcw(a2, ap, bp) < 0) chooseA = false;
else chooseA = isDinOABC(ap, bp, b2, a2);
}
if(chooseA){
nxeap = edges[edges[eap].rev].ccw;
auto [hab, hba] = newEdge(a2, bp);
insertCwAfter(hab, nxeap);
insertCcwAfter(hba, ebp);
eap = nxeap; ap = a2;
}
else {
nxebp = edges[edges[ebp].rev].cw;
auto [hba, hab] = newEdge(b2, ap);
insertCcwAfter(hba, nxebp);
insertCwAfter(hab, eap);
ebp = nxebp; 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 posbuf = pos;
int posptr = 0;
mappings.assign(sz, 0);
for(int i=0; i<sz; i++){
int v = pi[i];
if(i == 0 || !(posbuf[pi[posptr-1]] == posbuf[v])){
pi[posptr] = v;
pos[posptr++] = posbuf[v];
mappings[v] = v;
} else {
mappings[v] = pi[posptr-1];
}
}
if(posptr >= 2) outerOneEdge = solveRange(0, posptr).second;
std::swap(pos, posbuf);
for(auto& e : edges) e.to = pi[e.to];
}
std::vector<int> openAddress;
std::vector<GPos2> pos;
std::vector<Edge> edges;
std::vector<int> mappings;
int outerOneEdge = -1;
public:
DelaunayTriangulation()
: pos()
{ solve(); }
DelaunayTriangulation(std::vector<GPos2> x_points)
: pos(std::move(x_points))
{
solve();
}
std::vector<std::pair<int, int>> getEdges() const {
std::vector<std::pair<int, int>> res;
for(int e=0; e<(int)edges.size(); e++) if(edges[e].enabled){
int re = edges[e].rev;
if(e < re) continue;
res.push_back({ edges[e].to, edges[re].to });
}
for(int v=0; v<int(mappings.size()); v++){
if(mappings[v] != v) res.push_back({ v, mappings[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.
inline std::vector<std::pair<int, int>>
euclidean_mst(const std::vector<point<long long>> &points) {
using internal_point =
euclidean_mst_internal::VecI2<long long, long long>;
using triangulation = euclidean_mst_internal::DelaunayTriangulation<
long long, long long, __int128_t>;
std::vector<internal_point> converted;
converted.reserve(points.size());
for (const auto &value : points) {
converted.emplace_back(value.x, value.y);
}
std::vector<std::pair<int, int>> candidates =
triangulation(std::move(converted)).getEdges();
auto squared_distance = [&](const std::pair<int, int> &edge) {
__int128_t dx = __int128_t(points[edge.first].x) - points[edge.second].x;
__int128_t dy = __int128_t(points[edge.first].y) - points[edge.second].y;
return dx * dx + dy * dy;
};
std::stable_sort(candidates.begin(), candidates.end(),
[&](const auto &first, const auto &second) {
return squared_distance(first) < squared_distance(second);
});
atcoder::dsu components(int(points.size()));
std::vector<std::pair<int, int>> result;
result.reserve(points.empty() ? 0 : points.size() - 1);
for (auto edge : candidates) {
if (components.same(edge.first, edge.second)) {
continue;
}
components.merge(edge.first, edge.second);
if (edge.first > edge.second) {
std::swap(edge.first, edge.second);
}
result.push_back(edge);
}
std::sort(result.begin(), result.end());
return result;
}
} // namespace noya