Euclidean Minimum Spanning Tree
(geometry/euclidean_mst.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/euclidean_mst.hpp"
Overview
T may be a built-in integer or an exact coordinate class such as
math::Rational<long long> or math::Rational<utilities::BigInt>.
Rational intermediate arithmetic uses T without floating-point conversion.
All intermediate fractions must be representable. Listed complexities count
scalar operations; rational arithmetic adds its gcd and integer arithmetic costs.
Returned lengths or constructed coordinates still use long double.
This header constructs a minimum spanning tree of two-dimensional exact-coordinate points under Euclidean distance. It first builds a Delaunay triangulation by divide and conquer, then applies Kruskal’s algorithm to its $O(N)$ edges.
Duplicate points are connected to one representative by zero-length edges. Collinear and cocircular point sets are supported.
Types
template <class T>
struct EuclideanMstEdge {
int from;
int to;
T squared_distance;
};
template <class T>
struct EuclideanMst {
long double cost;
std::vector<EuclideanMstEdge<T>> edges;
};
For integral input coordinates, T in the returned types is wide_type<T>,
which is __int128_t. Edge selection compares squared_distance exactly.
EuclideanMst::cost is the sum of the selected Euclidean lengths and is stored
as long double.
All orientation, squared-distance, and degree-four incircle expressions must fit in signed 128-bit arithmetic. This condition holds for the Library Checker constraints.
Functions
| Function | Description | Complexity |
|---|---|---|
template <ExactCoordinate T> std::vector<EuclideanMstEdge<wide_type<T>>> euclidean_mst_edges(const std::vector<Point<T>>& points) |
Returns the $O(N)$ Delaunay and duplicate edges containing a Euclidean MST. | $O(N\log N)$ time and $O(N)$ memory |
template <ExactCoordinate T> EuclideanMst<wide_type<T>> euclidean_mst(const std::vector<Point<T>>& points) |
Returns the MST length and its selected edges. | $O(N\log N)$ time and $O(N)$ memory |
Vertices are identified by their zero-based indices in points. For zero or
one point, the result has cost zero and no edges. For $N$ points with $N>0$,
euclidean_mst(points).edges contains exactly $N-1$ edges.
Example
#include "geometry/euclidean_mst.hpp"
#include <iostream>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> points;
points.emplace_back(0, 0);
points.emplace_back(3, 0);
points.emplace_back(0, 4);
auto mst = m1une::geometry::euclidean_mst(points);
std::cout << static_cast<double>(mst.cost) << "\n";
for (const auto& edge : mst.edges) {
std::cout << edge.from << ' ' << edge.to << "\n";
}
}
Depends on
- DSU (Disjoint Set Union)
(ds/dsu/dsu.hpp)
- geometry/detail/floating_predicate.hpp
- 2D Point and Predicates
(geometry/point.hpp)
Required by
- Geometry Bundle
(geometry/all.hpp)
- Delaunay Triangulation
(geometry/delaunay_triangulation.hpp)
- Voronoi Diagram
(geometry/voronoi_diagram.hpp)
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/delaunay_triangulation.test.cpp
- verify/geometry/euclidean_mst.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/rational.test.cpp
- verify/geometry/voronoi_diagram.test.cpp
Code
#ifndef M1UNE_GEOMETRY_EUCLIDEAN_MST_HPP
#define M1UNE_GEOMETRY_EUCLIDEAN_MST_HPP 1
#include <algorithm>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <limits>
#include <tuple>
#include <utility>
#include <vector>
#include "../ds/dsu/dsu.hpp"
#include "point.hpp"
namespace m1une {
namespace geometry {
template <class T>
struct EuclideanMstEdge {
int from;
int to;
T squared_distance;
};
template <class T>
struct EuclideanMst {
long double cost;
std::vector<EuclideanMstEdge<T>> edges;
};
namespace detail {
template <ExactCoordinate T>
class EuclideanDelaunay {
private:
using W = wide_type<T>;
struct InternalPoint {
W x;
W y;
friend bool operator==(const InternalPoint&, const InternalPoint&) = default;
};
struct Edge {
int to;
int ccw;
int cw;
int reverse;
bool enabled = false;
};
std::vector<int> open_addresses;
std::vector<InternalPoint> points;
std::vector<Edge> edges;
std::vector<int> duplicate_representative;
static InternalPoint subtract(const InternalPoint& a, const InternalPoint& b) {
return InternalPoint{a.x - b.x, a.y - b.y};
}
static W cross_product(const InternalPoint& a, const InternalPoint& b) {
return a.x * b.y - a.y * b.x;
}
static W squared_norm(const InternalPoint& point) {
return point.x * point.x + point.y * point.y;
}
static bool inside_circumcircle(
InternalPoint a,
InternalPoint b,
InternalPoint c,
const InternalPoint& d
) {
a = subtract(a, d);
b = subtract(b, d);
c = subtract(c, d);
W determinant = cross_product(b, c) * squared_norm(a)
+ cross_product(c, a) * squared_norm(b)
+ cross_product(a, b) * squared_norm(c);
return determinant > 0;
}
int get_open_address() {
if (open_addresses.empty()) {
edges.push_back(Edge());
return int(edges.size()) - 1;
}
int result = open_addresses.back();
open_addresses.pop_back();
return result;
}
std::pair<int, int> add_edge(int from, int to) {
int forward = get_open_address();
int backward = get_open_address();
edges[forward].to = to;
edges[forward].ccw = forward;
edges[forward].cw = forward;
edges[forward].reverse = backward;
edges[forward].enabled = true;
edges[backward].to = from;
edges[backward].ccw = backward;
edges[backward].cw = backward;
edges[backward].reverse = forward;
edges[backward].enabled = true;
return {forward, backward};
}
void erase_directed_edge(int edge) {
int ccw = edges[edge].ccw;
int cw = edges[edge].cw;
edges[ccw].cw = cw;
edges[cw].ccw = ccw;
edges[edge].enabled = false;
}
void erase_edge(int edge) {
int reverse = edges[edge].reverse;
erase_directed_edge(edge);
erase_directed_edge(reverse);
open_addresses.push_back(edge);
open_addresses.push_back(reverse);
}
void insert_ccw_after(int edge, int position) {
int next = edges[position].ccw;
edges[edge].ccw = next;
edges[next].cw = edge;
edges[edge].cw = position;
edges[position].ccw = edge;
}
void insert_cw_after(int edge, int position) {
int next = edges[position].cw;
edges[edge].cw = next;
edges[next].ccw = edge;
edges[edge].ccw = position;
edges[position].cw = edge;
}
int orientation(int a, int b, int c) const {
InternalPoint ab = subtract(points[b], points[a]);
InternalPoint ac = subtract(points[c], points[a]);
W value = cross_product(ab, ac);
return (value > 0) - (value < 0);
}
std::pair<int, int> go_next(int edge) const {
int vertex = edges[edge].to;
int next_edge = edges[edges[edge].reverse].ccw;
return {vertex, next_edge};
}
std::pair<int, int> go_previous(int edge) const {
int vertex = edges[edges[edge].cw].to;
int next_edge = edges[edges[edge].cw].reverse;
return {vertex, next_edge};
}
std::tuple<int, int, int, int> lower_tangent(
int left_vertex,
int left_edge,
int right_vertex,
int right_edge
) const {
while (true) {
auto [next_left_vertex, next_left_edge] = go_previous(left_edge);
if (orientation(right_vertex, left_vertex, next_left_vertex) > 0) {
left_vertex = next_left_vertex;
left_edge = next_left_edge;
continue;
}
auto [next_right_vertex, next_right_edge] = go_next(right_edge);
if (orientation(left_vertex, right_vertex, next_right_vertex) < 0) {
right_vertex = next_right_vertex;
right_edge = next_right_edge;
continue;
}
break;
}
return {left_vertex, left_edge, right_vertex, right_edge};
}
std::pair<int, int> extreme_vertex(int vertex, int edge, bool minimum) const {
std::pair<int, int> result = {vertex, edge};
int current_vertex = vertex;
int current_edge = edge;
do {
std::tie(current_vertex, current_edge) = go_next(current_edge);
std::pair<int, int> candidate = {current_vertex, current_edge};
if ((minimum && candidate < result) || (!minimum && result < candidate)) {
result = candidate;
}
} while (current_edge != edge);
return result;
}
bool inside_circumcircle(int a, int b, int c, int d) const {
return inside_circumcircle(points[a], points[b], points[c], points[d]);
}
std::pair<int, int> merge_triangulations(
int left_vertex,
int left_edge,
int right_vertex,
int right_edge
) {
std::tie(left_vertex, left_edge) = extreme_vertex(left_vertex, left_edge, false);
std::tie(right_vertex, right_edge) = extreme_vertex(right_vertex, right_edge, true);
auto [lower_left, lower_left_edge, lower_right, lower_right_edge]
= lower_tangent(left_vertex, left_edge, right_vertex, right_edge);
auto [upper_right, upper_right_edge, upper_left, upper_left_edge]
= lower_tangent(right_vertex, right_edge, left_vertex, left_edge);
lower_right_edge = edges[lower_right_edge].cw;
upper_right_edge = edges[upper_right_edge].cw;
auto [base, reverse_base] = add_edge(lower_left, lower_right);
insert_cw_after(base, lower_left_edge);
insert_ccw_after(reverse_base, lower_right_edge);
if (lower_left == upper_left) upper_left_edge = base;
if (lower_right == upper_right) upper_right_edge = reverse_base;
int left = lower_left;
int left_candidate = lower_left_edge;
int right = lower_right;
int right_candidate = lower_right_edge;
while (left != upper_left || right != upper_right) {
int next_left = edges[left_candidate].to;
int next_right = edges[right_candidate].to;
int next_left_candidate = edges[left_candidate].ccw;
int next_right_candidate = edges[right_candidate].cw;
if (left_candidate != upper_left_edge && next_left_candidate != base) {
int second_left = edges[next_left_candidate].to;
if (inside_circumcircle(left, right, next_left, second_left)) {
erase_edge(left_candidate);
left_candidate = next_left_candidate;
continue;
}
}
if (right_candidate != upper_right_edge && next_right_candidate != reverse_base) {
int second_right = edges[next_right_candidate].to;
if (inside_circumcircle(next_right, left, right, second_right)) {
erase_edge(right_candidate);
right_candidate = next_right_candidate;
continue;
}
}
bool choose_left = right_candidate == upper_right_edge;
if (left_candidate != upper_left_edge && right_candidate != upper_right_edge) {
if (orientation(left, right, next_right) < 0) {
choose_left = true;
} else if (orientation(next_left, left, right) < 0) {
choose_left = false;
} else {
choose_left = inside_circumcircle(left, right, next_right, next_left);
}
}
if (choose_left) {
next_left_candidate = edges[edges[left_candidate].reverse].ccw;
auto [new_base, new_reverse_base] = add_edge(next_left, right);
insert_cw_after(new_base, next_left_candidate);
insert_ccw_after(new_reverse_base, right_candidate);
left_candidate = next_left_candidate;
left = next_left;
} else {
next_right_candidate = edges[edges[right_candidate].reverse].cw;
auto [new_reverse_base, new_base] = add_edge(next_right, left);
insert_ccw_after(new_reverse_base, next_right_candidate);
insert_cw_after(new_base, left_candidate);
right_candidate = next_right_candidate;
right = next_right;
}
}
return {lower_left, base};
}
std::pair<int, int> solve_range(int left, int right) {
if (right - left == 2) {
auto [forward, backward] = add_edge(left, left + 1);
(void)backward;
return {left, forward};
}
if (right - left == 3) {
int middle = left + 1;
int last = left + 2;
auto [first_middle, middle_first] = add_edge(left, middle);
auto [middle_last, last_middle] = add_edge(middle, last);
int direction = orientation(left, middle, last);
if (direction == 0) {
insert_ccw_after(middle_first, middle_last);
return {left, first_middle};
}
auto [first_last, last_first] = add_edge(left, last);
if (direction > 0) {
insert_cw_after(first_middle, first_last);
insert_cw_after(middle_last, middle_first);
insert_cw_after(last_first, last_middle);
return {left, first_middle};
}
insert_ccw_after(first_middle, first_last);
insert_ccw_after(middle_last, middle_first);
insert_ccw_after(last_first, last_middle);
return {middle, middle_first};
}
int middle = (left + right) / 2;
auto [left_vertex, left_edge] = solve_range(left, middle);
auto [right_vertex, right_edge] = solve_range(middle, right);
return merge_triangulations(left_vertex, left_edge, right_vertex, right_edge);
}
void solve() {
int size = int(points.size());
if (size <= 1) return;
std::vector<int> order(size);
for (int i = 0; i < size; i++) order[i] = i;
std::stable_sort(order.begin(), order.end(), [&](int left, int right) {
if (points[left].x != points[right].x) {
return points[left].x < points[right].x;
}
return points[left].y < points[right].y;
});
std::vector<InternalPoint> original_points = points;
duplicate_representative.assign(size, 0);
int unique_size = 0;
for (int i = 0; i < size; i++) {
int vertex = order[i];
if (i == 0 || !(original_points[order[unique_size - 1]] == original_points[vertex])) {
order[unique_size] = vertex;
points[unique_size] = original_points[vertex];
unique_size++;
duplicate_representative[vertex] = vertex;
} else {
duplicate_representative[vertex] = order[unique_size - 1];
}
}
if (unique_size >= 2) solve_range(0, unique_size);
points.swap(original_points);
for (auto& edge : edges) edge.to = order[edge.to];
}
public:
explicit EuclideanDelaunay(const std::vector<Point<T>>& input_points) {
assert(input_points.size() <= std::size_t(std::numeric_limits<int>::max()));
points.reserve(input_points.size());
edges.reserve(std::size_t(6) * input_points.size());
for (const auto& point : input_points) {
points.push_back(InternalPoint{W(point.x), W(point.y)});
}
solve();
}
bool has_duplicates() const {
for (
int vertex = 0;
vertex < int(duplicate_representative.size());
++vertex
) {
if (duplicate_representative[vertex] != vertex) return true;
}
return false;
}
std::vector<std::pair<int, int>> get_edges() const {
std::vector<std::pair<int, int>> result;
result.reserve(edges.size() / 2 + duplicate_representative.size());
for (int edge = 0; edge < int(edges.size()); edge++) {
if (!edges[edge].enabled) continue;
int reverse = edges[edge].reverse;
if (edge < reverse) continue;
result.emplace_back(edges[edge].to, edges[reverse].to);
}
for (int vertex = 0; vertex < int(duplicate_representative.size()); vertex++) {
if (duplicate_representative[vertex] != vertex) {
result.emplace_back(vertex, duplicate_representative[vertex]);
}
}
return result;
}
};
} // namespace detail
// Returns O(n) Delaunay edges containing a Euclidean minimum spanning tree.
template <ExactCoordinate T>
std::vector<EuclideanMstEdge<wide_type<T>>> euclidean_mst_edges(
const std::vector<Point<T>>& points
) {
using W = wide_type<T>;
auto delaunay_edges = detail::EuclideanDelaunay<T>(points).get_edges();
std::vector<EuclideanMstEdge<W>> result;
result.reserve(delaunay_edges.size());
for (auto [from, to] : delaunay_edges) {
result.push_back(EuclideanMstEdge<W>{from, to, distance2(points[from], points[to])});
}
return result;
}
// Returns a Euclidean minimum spanning tree.
template <ExactCoordinate T>
EuclideanMst<wide_type<T>> euclidean_mst(const std::vector<Point<T>>& points) {
using W = wide_type<T>;
auto candidates = euclidean_mst_edges(points);
std::sort(candidates.begin(), candidates.end(), [](const auto& left, const auto& right) {
if (left.squared_distance != right.squared_distance) {
return left.squared_distance < right.squared_distance;
}
if (left.from != right.from) return left.from < right.from;
return left.to < right.to;
});
m1une::ds::Dsu dsu(int(points.size()));
EuclideanMst<W> result;
result.cost = 0;
result.edges.reserve(points.empty() ? 0 : points.size() - 1);
for (const auto& edge : candidates) {
if (dsu.same(edge.from, edge.to)) continue;
dsu.merge(edge.from, edge.to);
result.cost += std::sqrt(static_cast<long double>(edge.squared_distance));
result.edges.push_back(edge);
if (result.edges.size() + 1 == points.size()) break;
}
assert(points.empty() || result.edges.size() + 1 == points.size());
return result;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_EUCLIDEAN_MST_HPP#line 1 "geometry/euclidean_mst.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <limits>
#include <tuple>
#include <utility>
#include <vector>
#line 1 "ds/dsu/dsu.hpp"
#line 5 "ds/dsu/dsu.hpp"
#include <numeric>
#line 8 "ds/dsu/dsu.hpp"
namespace m1une {
namespace ds {
struct Dsu {
private:
int _n;
// parent_or_size[i] is the parent of i if it's >= 0.
// If it's < 0, then i is a root and -parent_or_size[i] is the size of the group.
std::vector<int> parent_or_size;
// Returns {new leader, absorbed leader}. The absorbed leader is -1 when
// both vertices already belong to the same component.
std::pair<int, int> merge_leaders(int a, int b) {
int x = leader(a), y = leader(b);
if (x == y) return {x, -1};
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, y};
}
public:
Dsu() : _n(0) {}
explicit Dsu(int n) : _n(n), parent_or_size(n, -1) {}
// Merges the group containing 'a' with the group containing 'b'.
// Returns the leader of the merged group.
int merge(int a, int b) {
return merge_leaders(a, b).first;
}
// Invokes callback(new_leader, absorbed_leader) after an actual merge.
// Returns the leader of the merged group.
template <class Callback>
int merge(int a, int b, Callback&& callback) {
std::pair<int, int> merged = merge_leaders(a, b);
if (merged.second != -1) callback(merged.first, merged.second);
return merged.first;
}
// Returns true if 'a' and 'b' belong to the same group.
bool same(int a, int b) {
return leader(a) == leader(b);
}
// Returns the leader (representative) of the group containing 'a'.
int leader(int a) {
if (parent_or_size[a] < 0) return a;
// Path compression
return parent_or_size[a] = leader(parent_or_size[a]);
}
// Returns the size of the group containing 'a'.
int size(int a) {
return -parent_or_size[leader(a)];
}
// Returns a list of all groups, where each group is a vector of its elements.
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;
}
};
} // namespace ds
} // namespace m1une
#line 1 "geometry/point.hpp"
#line 7 "geometry/point.hpp"
#include <type_traits>
#line 1 "geometry/detail/floating_predicate.hpp"
namespace m1une {
namespace geometry {
namespace predicate_detail {
template <typename T>
constexpr T absolute(T value) {
return value < T(0) ? -value : value;
}
template <typename T>
constexpr T max_value(T first, T second) {
return first < second ? second : first;
}
template <typename T>
constexpr T vector_scale(T x, T y) {
return max_value(absolute(x), absolute(y));
}
template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
if constexpr (Exact) {
return (value > T(0)) - (value < T(0));
} else {
const T tolerance = T(eps) * scale;
return (value > tolerance) - (value < -tolerance);
}
}
template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
if constexpr (Exact) {
return T(0);
} else {
return vector_scale(ax, ay) * vector_scale(bx, by);
}
}
template <bool Exact, typename T>
constexpr int determinant_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T determinant = ax * by - ay * bx;
return scaled_sign<Exact>(
determinant,
determinant_scale<Exact>(ax, ay, bx, by),
eps
);
}
template <bool Exact, typename T>
constexpr int orientation_sign(
T direction_x,
T direction_y,
T offset_x,
T offset_y,
long double eps
) {
const T determinant =
direction_x * offset_y - direction_y * offset_x;
T scale = T(0);
if constexpr (!Exact) {
const T direction_scale =
vector_scale(direction_x, direction_y);
scale = direction_scale * max_value(
direction_scale,
vector_scale(offset_x, offset_y)
);
}
return scaled_sign<Exact>(determinant, scale, eps);
}
template <bool Exact, typename T>
constexpr int dot_sign(
T ax,
T ay,
T bx,
T by,
long double eps
) {
const T value = ax * bx + ay * by;
T scale = T(0);
if constexpr (!Exact) {
scale = vector_scale(ax, ay) * vector_scale(bx, by);
}
return scaled_sign<Exact>(value, scale, eps);
}
} // namespace predicate_detail
} // namespace geometry
} // namespace m1une
#line 10 "geometry/point.hpp"
namespace m1une {
namespace geometry {
template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
(std::is_arithmetic_v<T> ||
(std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
T(0);
T(1);
static_cast<long double>(a);
{ +a } -> std::same_as<T>;
{ -a } -> std::same_as<T>;
{ a + b } -> std::same_as<T>;
{ a - b } -> std::same_as<T>;
{ a * b } -> std::same_as<T>;
{ a / b } -> std::same_as<T>;
{ a += b } -> std::same_as<T&>;
{ a -= b } -> std::same_as<T&>;
}));
// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;
template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
std::conditional_t<std::floating_point<T>, long double, T>>;
template <Coordinate T>
struct Point {
T x;
T y;
constexpr Point() : x(0), y(0) {}
constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}
template <Coordinate U>
explicit constexpr Point(const Point<U>& other)
: x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}
constexpr Point& operator+=(const Point& other) {
x += other.x;
y += other.y;
return *this;
}
constexpr Point& operator-=(const Point& other) {
x -= other.x;
y -= other.y;
return *this;
}
constexpr Point operator+() const {
return *this;
}
constexpr Point operator-() const {
return Point(-x, -y);
}
friend constexpr Point operator+(Point left, const Point& right) {
return left += right;
}
friend constexpr Point operator-(Point left, const Point& right) {
return left -= right;
}
friend constexpr bool operator==(const Point&, const Point&) = default;
friend constexpr bool operator<(const Point& left, const Point& right) {
if (left.x != right.x) return left.x < right.x;
return left.y < right.y;
}
};
template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
return Point<long double>(point);
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) * Result(scalar),
Result(point.y) * Result(scalar)
);
}
template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
return point * scalar;
}
template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
using Result = std::common_type_t<T, Scalar>;
return Point<Result>(
Result(point.x) / Result(scalar),
Result(point.y) / Result(scalar)
);
}
template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}
template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}
template <Coordinate T>
constexpr wide_type<T> cross(
const Point<T>& origin,
const Point<T>& a,
const Point<T>& b
) {
using W = wide_type<T>;
W ax = W(a.x) - W(origin.x);
W ay = W(a.y) - W(origin.y);
W bx = W(b.x) - W(origin.x);
W by = W(b.y) - W(origin.y);
return ax * by - ay * bx;
}
template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
return dot(point, point);
}
template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
using W = wide_type<T>;
W dx = W(a.x) - W(b.x);
W dy = W(a.y) - W(b.y);
return dx * dx + dy * dy;
}
template <Coordinate T>
long double norm(const Point<T>& point) {
return std::hypot(
static_cast<long double>(point.x),
static_cast<long double>(point.y)
);
}
template <Coordinate 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)
);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio + second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
(std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
const Point<T>& a,
const Point<T>& b,
M m,
N n
) {
long double first_ratio = static_cast<long double>(m);
long double second_ratio = static_cast<long double>(n);
long double denominator = first_ratio - second_ratio;
assert(denominator != 0);
Point<long double> first(a);
Point<long double> direction = Point<long double>(b) - first;
return first + direction * (first_ratio / denominator);
}
template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
return predicate_detail::scaled_sign<ExactCoordinate<T>>(
value,
wide_type<T>(1),
eps
);
}
template <Coordinate T>
constexpr int orientation(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
using W = wide_type<T>;
const W first_x = W(b.x) - W(a.x);
const W first_y = W(b.y) - W(a.y);
const W second_x = W(c.x) - W(a.x);
const W second_y = W(c.y) - W(a.y);
return predicate_detail::orientation_sign<ExactCoordinate<T>>(
first_x,
first_y,
second_x,
second_y,
eps
);
}
template <Coordinate T>
constexpr bool collinear(
const Point<T>& a,
const Point<T>& b,
const Point<T>& c,
long double eps = 1e-12L
) {
return orientation(a, b, c, eps) == 0;
}
template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
long double cosine = std::cos(angle);
long double sine = std::sin(angle);
return Point<long double>(
static_cast<long double>(point.x) * cosine -
static_cast<long double>(point.y) * sine,
static_cast<long double>(point.x) * sine +
static_cast<long double>(point.y) * cosine
);
}
template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
long double length = norm(point);
assert(length != 0);
return Point<long double>(
static_cast<long double>(point.x) / length,
static_cast<long double>(point.y) / length
);
}
} // namespace geometry
} // namespace m1une
#line 16 "geometry/euclidean_mst.hpp"
namespace m1une {
namespace geometry {
template <class T>
struct EuclideanMstEdge {
int from;
int to;
T squared_distance;
};
template <class T>
struct EuclideanMst {
long double cost;
std::vector<EuclideanMstEdge<T>> edges;
};
namespace detail {
template <ExactCoordinate T>
class EuclideanDelaunay {
private:
using W = wide_type<T>;
struct InternalPoint {
W x;
W y;
friend bool operator==(const InternalPoint&, const InternalPoint&) = default;
};
struct Edge {
int to;
int ccw;
int cw;
int reverse;
bool enabled = false;
};
std::vector<int> open_addresses;
std::vector<InternalPoint> points;
std::vector<Edge> edges;
std::vector<int> duplicate_representative;
static InternalPoint subtract(const InternalPoint& a, const InternalPoint& b) {
return InternalPoint{a.x - b.x, a.y - b.y};
}
static W cross_product(const InternalPoint& a, const InternalPoint& b) {
return a.x * b.y - a.y * b.x;
}
static W squared_norm(const InternalPoint& point) {
return point.x * point.x + point.y * point.y;
}
static bool inside_circumcircle(
InternalPoint a,
InternalPoint b,
InternalPoint c,
const InternalPoint& d
) {
a = subtract(a, d);
b = subtract(b, d);
c = subtract(c, d);
W determinant = cross_product(b, c) * squared_norm(a)
+ cross_product(c, a) * squared_norm(b)
+ cross_product(a, b) * squared_norm(c);
return determinant > 0;
}
int get_open_address() {
if (open_addresses.empty()) {
edges.push_back(Edge());
return int(edges.size()) - 1;
}
int result = open_addresses.back();
open_addresses.pop_back();
return result;
}
std::pair<int, int> add_edge(int from, int to) {
int forward = get_open_address();
int backward = get_open_address();
edges[forward].to = to;
edges[forward].ccw = forward;
edges[forward].cw = forward;
edges[forward].reverse = backward;
edges[forward].enabled = true;
edges[backward].to = from;
edges[backward].ccw = backward;
edges[backward].cw = backward;
edges[backward].reverse = forward;
edges[backward].enabled = true;
return {forward, backward};
}
void erase_directed_edge(int edge) {
int ccw = edges[edge].ccw;
int cw = edges[edge].cw;
edges[ccw].cw = cw;
edges[cw].ccw = ccw;
edges[edge].enabled = false;
}
void erase_edge(int edge) {
int reverse = edges[edge].reverse;
erase_directed_edge(edge);
erase_directed_edge(reverse);
open_addresses.push_back(edge);
open_addresses.push_back(reverse);
}
void insert_ccw_after(int edge, int position) {
int next = edges[position].ccw;
edges[edge].ccw = next;
edges[next].cw = edge;
edges[edge].cw = position;
edges[position].ccw = edge;
}
void insert_cw_after(int edge, int position) {
int next = edges[position].cw;
edges[edge].cw = next;
edges[next].ccw = edge;
edges[edge].ccw = position;
edges[position].cw = edge;
}
int orientation(int a, int b, int c) const {
InternalPoint ab = subtract(points[b], points[a]);
InternalPoint ac = subtract(points[c], points[a]);
W value = cross_product(ab, ac);
return (value > 0) - (value < 0);
}
std::pair<int, int> go_next(int edge) const {
int vertex = edges[edge].to;
int next_edge = edges[edges[edge].reverse].ccw;
return {vertex, next_edge};
}
std::pair<int, int> go_previous(int edge) const {
int vertex = edges[edges[edge].cw].to;
int next_edge = edges[edges[edge].cw].reverse;
return {vertex, next_edge};
}
std::tuple<int, int, int, int> lower_tangent(
int left_vertex,
int left_edge,
int right_vertex,
int right_edge
) const {
while (true) {
auto [next_left_vertex, next_left_edge] = go_previous(left_edge);
if (orientation(right_vertex, left_vertex, next_left_vertex) > 0) {
left_vertex = next_left_vertex;
left_edge = next_left_edge;
continue;
}
auto [next_right_vertex, next_right_edge] = go_next(right_edge);
if (orientation(left_vertex, right_vertex, next_right_vertex) < 0) {
right_vertex = next_right_vertex;
right_edge = next_right_edge;
continue;
}
break;
}
return {left_vertex, left_edge, right_vertex, right_edge};
}
std::pair<int, int> extreme_vertex(int vertex, int edge, bool minimum) const {
std::pair<int, int> result = {vertex, edge};
int current_vertex = vertex;
int current_edge = edge;
do {
std::tie(current_vertex, current_edge) = go_next(current_edge);
std::pair<int, int> candidate = {current_vertex, current_edge};
if ((minimum && candidate < result) || (!minimum && result < candidate)) {
result = candidate;
}
} while (current_edge != edge);
return result;
}
bool inside_circumcircle(int a, int b, int c, int d) const {
return inside_circumcircle(points[a], points[b], points[c], points[d]);
}
std::pair<int, int> merge_triangulations(
int left_vertex,
int left_edge,
int right_vertex,
int right_edge
) {
std::tie(left_vertex, left_edge) = extreme_vertex(left_vertex, left_edge, false);
std::tie(right_vertex, right_edge) = extreme_vertex(right_vertex, right_edge, true);
auto [lower_left, lower_left_edge, lower_right, lower_right_edge]
= lower_tangent(left_vertex, left_edge, right_vertex, right_edge);
auto [upper_right, upper_right_edge, upper_left, upper_left_edge]
= lower_tangent(right_vertex, right_edge, left_vertex, left_edge);
lower_right_edge = edges[lower_right_edge].cw;
upper_right_edge = edges[upper_right_edge].cw;
auto [base, reverse_base] = add_edge(lower_left, lower_right);
insert_cw_after(base, lower_left_edge);
insert_ccw_after(reverse_base, lower_right_edge);
if (lower_left == upper_left) upper_left_edge = base;
if (lower_right == upper_right) upper_right_edge = reverse_base;
int left = lower_left;
int left_candidate = lower_left_edge;
int right = lower_right;
int right_candidate = lower_right_edge;
while (left != upper_left || right != upper_right) {
int next_left = edges[left_candidate].to;
int next_right = edges[right_candidate].to;
int next_left_candidate = edges[left_candidate].ccw;
int next_right_candidate = edges[right_candidate].cw;
if (left_candidate != upper_left_edge && next_left_candidate != base) {
int second_left = edges[next_left_candidate].to;
if (inside_circumcircle(left, right, next_left, second_left)) {
erase_edge(left_candidate);
left_candidate = next_left_candidate;
continue;
}
}
if (right_candidate != upper_right_edge && next_right_candidate != reverse_base) {
int second_right = edges[next_right_candidate].to;
if (inside_circumcircle(next_right, left, right, second_right)) {
erase_edge(right_candidate);
right_candidate = next_right_candidate;
continue;
}
}
bool choose_left = right_candidate == upper_right_edge;
if (left_candidate != upper_left_edge && right_candidate != upper_right_edge) {
if (orientation(left, right, next_right) < 0) {
choose_left = true;
} else if (orientation(next_left, left, right) < 0) {
choose_left = false;
} else {
choose_left = inside_circumcircle(left, right, next_right, next_left);
}
}
if (choose_left) {
next_left_candidate = edges[edges[left_candidate].reverse].ccw;
auto [new_base, new_reverse_base] = add_edge(next_left, right);
insert_cw_after(new_base, next_left_candidate);
insert_ccw_after(new_reverse_base, right_candidate);
left_candidate = next_left_candidate;
left = next_left;
} else {
next_right_candidate = edges[edges[right_candidate].reverse].cw;
auto [new_reverse_base, new_base] = add_edge(next_right, left);
insert_ccw_after(new_reverse_base, next_right_candidate);
insert_cw_after(new_base, left_candidate);
right_candidate = next_right_candidate;
right = next_right;
}
}
return {lower_left, base};
}
std::pair<int, int> solve_range(int left, int right) {
if (right - left == 2) {
auto [forward, backward] = add_edge(left, left + 1);
(void)backward;
return {left, forward};
}
if (right - left == 3) {
int middle = left + 1;
int last = left + 2;
auto [first_middle, middle_first] = add_edge(left, middle);
auto [middle_last, last_middle] = add_edge(middle, last);
int direction = orientation(left, middle, last);
if (direction == 0) {
insert_ccw_after(middle_first, middle_last);
return {left, first_middle};
}
auto [first_last, last_first] = add_edge(left, last);
if (direction > 0) {
insert_cw_after(first_middle, first_last);
insert_cw_after(middle_last, middle_first);
insert_cw_after(last_first, last_middle);
return {left, first_middle};
}
insert_ccw_after(first_middle, first_last);
insert_ccw_after(middle_last, middle_first);
insert_ccw_after(last_first, last_middle);
return {middle, middle_first};
}
int middle = (left + right) / 2;
auto [left_vertex, left_edge] = solve_range(left, middle);
auto [right_vertex, right_edge] = solve_range(middle, right);
return merge_triangulations(left_vertex, left_edge, right_vertex, right_edge);
}
void solve() {
int size = int(points.size());
if (size <= 1) return;
std::vector<int> order(size);
for (int i = 0; i < size; i++) order[i] = i;
std::stable_sort(order.begin(), order.end(), [&](int left, int right) {
if (points[left].x != points[right].x) {
return points[left].x < points[right].x;
}
return points[left].y < points[right].y;
});
std::vector<InternalPoint> original_points = points;
duplicate_representative.assign(size, 0);
int unique_size = 0;
for (int i = 0; i < size; i++) {
int vertex = order[i];
if (i == 0 || !(original_points[order[unique_size - 1]] == original_points[vertex])) {
order[unique_size] = vertex;
points[unique_size] = original_points[vertex];
unique_size++;
duplicate_representative[vertex] = vertex;
} else {
duplicate_representative[vertex] = order[unique_size - 1];
}
}
if (unique_size >= 2) solve_range(0, unique_size);
points.swap(original_points);
for (auto& edge : edges) edge.to = order[edge.to];
}
public:
explicit EuclideanDelaunay(const std::vector<Point<T>>& input_points) {
assert(input_points.size() <= std::size_t(std::numeric_limits<int>::max()));
points.reserve(input_points.size());
edges.reserve(std::size_t(6) * input_points.size());
for (const auto& point : input_points) {
points.push_back(InternalPoint{W(point.x), W(point.y)});
}
solve();
}
bool has_duplicates() const {
for (
int vertex = 0;
vertex < int(duplicate_representative.size());
++vertex
) {
if (duplicate_representative[vertex] != vertex) return true;
}
return false;
}
std::vector<std::pair<int, int>> get_edges() const {
std::vector<std::pair<int, int>> result;
result.reserve(edges.size() / 2 + duplicate_representative.size());
for (int edge = 0; edge < int(edges.size()); edge++) {
if (!edges[edge].enabled) continue;
int reverse = edges[edge].reverse;
if (edge < reverse) continue;
result.emplace_back(edges[edge].to, edges[reverse].to);
}
for (int vertex = 0; vertex < int(duplicate_representative.size()); vertex++) {
if (duplicate_representative[vertex] != vertex) {
result.emplace_back(vertex, duplicate_representative[vertex]);
}
}
return result;
}
};
} // namespace detail
// Returns O(n) Delaunay edges containing a Euclidean minimum spanning tree.
template <ExactCoordinate T>
std::vector<EuclideanMstEdge<wide_type<T>>> euclidean_mst_edges(
const std::vector<Point<T>>& points
) {
using W = wide_type<T>;
auto delaunay_edges = detail::EuclideanDelaunay<T>(points).get_edges();
std::vector<EuclideanMstEdge<W>> result;
result.reserve(delaunay_edges.size());
for (auto [from, to] : delaunay_edges) {
result.push_back(EuclideanMstEdge<W>{from, to, distance2(points[from], points[to])});
}
return result;
}
// Returns a Euclidean minimum spanning tree.
template <ExactCoordinate T>
EuclideanMst<wide_type<T>> euclidean_mst(const std::vector<Point<T>>& points) {
using W = wide_type<T>;
auto candidates = euclidean_mst_edges(points);
std::sort(candidates.begin(), candidates.end(), [](const auto& left, const auto& right) {
if (left.squared_distance != right.squared_distance) {
return left.squared_distance < right.squared_distance;
}
if (left.from != right.from) return left.from < right.from;
return left.to < right.to;
});
m1une::ds::Dsu dsu(int(points.size()));
EuclideanMst<W> result;
result.cost = 0;
result.edges.reserve(points.empty() ? 0 : points.size() - 1);
for (const auto& edge : candidates) {
if (dsu.same(edge.from, edge.to)) continue;
dsu.merge(edge.from, edge.to);
result.cost += std::sqrt(static_cast<long double>(edge.squared_distance));
result.edges.push_back(edge);
if (result.edges.size() + 1 == points.size()) break;
}
assert(points.empty() || result.edges.size() + 1 == points.size());
return result;
}
} // namespace geometry
} // namespace m1une