Manhattan Minimum Spanning Tree
(geometry/manhattan_mst.hpp)

• View this file on GitHub
• Last update: 2026-10-05 22:23:07+09:00
• Include: #include "geometry/manhattan_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.

This header constructs a minimum spanning tree of two-dimensional exact-coordinate points under Manhattan distance:

\[d(i,j)=|x_i-x_j|+|y_i-y_j|.\]

It uses four coordinate sweeps to generate only $O(N)$ candidate edges, then applies Kruskal’s algorithm.

Types

template <class T>
struct ManhattanMstEdge {
    int from;
    int to;
    T cost;
};

template <class T>
struct ManhattanMst {
    T cost;
    std::vector<ManhattanMstEdge<T>> edges;
};

For integral input coordinates, the result uses wide_type<T>, which is __int128_t. This prevents overflow during coordinate transformations and distance calculations when the input type is a standard integer type.

Functions

Function	Description	Complexity
manhattan_mst_edges(points)	Returns $O(N)$ candidate edges containing at least one Manhattan MST.	$O(N\log N)$
manhattan_mst(points)	Returns the MST cost and its selected edges.	$O(N\log N)$

Vertices are identified by their indices in points. Duplicate points are supported and may be joined by zero-cost edges. For zero or one point, the result has cost zero and no edges.

Example

#include "geometry/manhattan_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(2, 1);
    points.emplace_back(-1, 3);

    auto mst = m1une::geometry::manhattan_mst(points);
    long long cost = static_cast<long long>(mst.cost);
    std::cout << cost << "\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)

Verified with

• verify/geometry/centroid.test.cpp
• verify/geometry/geometry_algorithms.test.cpp
• verify/geometry/manhattan_mst.test.cpp
• verify/geometry/rational.test.cpp

Code

#ifndef M1UNE_GEOMETRY_MANHATTAN_MST_HPP
#define M1UNE_GEOMETRY_MANHATTAN_MST_HPP 1

#include <algorithm>
#include <cassert>
#include <concepts>
#include <limits>
#include <map>
#include <numeric>
#include <utility>
#include <vector>

#include "../ds/dsu/dsu.hpp"
#include "point.hpp"

namespace m1une {
namespace geometry {

template <class T>
struct ManhattanMstEdge {
    int from;
    int to;
    T cost;
};

template <class T>
struct ManhattanMst {
    T cost;
    std::vector<ManhattanMstEdge<T>> edges;
};

// Returns O(n) edges containing a Manhattan minimum spanning tree.
template <ExactCoordinate T>
std::vector<ManhattanMstEdge<wide_type<T>>> manhattan_mst_edges(const std::vector<Point<T>>& points) {
    using W = wide_type<T>;
    assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
    int n = int(points.size());
    std::vector<Point<W>> transformed;
    transformed.reserve(points.size());
    for (const auto& point : points) {
        transformed.emplace_back(W(point.x), W(point.y));
    }

    std::vector<int> order(n);
    std::iota(order.begin(), order.end(), 0);
    std::vector<ManhattanMstEdge<W>> edges;
    edges.reserve(std::size_t(4) * points.size());

    for (int direction = 0; direction < 4; direction++) {
        std::sort(order.begin(), order.end(), [&transformed](int i, int j) {
            W first = transformed[i].x + transformed[i].y;
            W second = transformed[j].x + transformed[j].y;
            if (first != second) return first < second;
            if (transformed[i].x != transformed[j].x) {
                return transformed[i].x < transformed[j].x;
            }
            return i < j;
        });

        std::map<W, int> sweep;
        for (int i : order) {
            auto it = sweep.lower_bound(-transformed[i].y);
            while (it != sweep.end()) {
                int j = it->second;
                if (transformed[i].x - transformed[j].x < transformed[i].y - transformed[j].y) {
                    break;
                }

                W dx = W(points[i].x) - W(points[j].x);
                W dy = W(points[i].y) - W(points[j].y);
                if (dx < 0) dx = -dx;
                if (dy < 0) dy = -dy;
                edges.push_back(ManhattanMstEdge<W>{i, j, dx + dy});
                it = sweep.erase(it);
            }
            sweep[-transformed[i].y] = i;
        }

        for (auto& point : transformed) {
            if (direction & 1) {
                point.x = -point.x;
            } else {
                std::swap(point.x, point.y);
            }
        }
    }
    return edges;
}

// Returns a Manhattan minimum spanning tree.
template <ExactCoordinate T>
ManhattanMst<wide_type<T>> manhattan_mst(const std::vector<Point<T>>& points) {
    using W = wide_type<T>;
    auto candidates = manhattan_mst_edges(points);
    std::sort(candidates.begin(), candidates.end(), [](const auto& a, const auto& b) { return a.cost < b.cost; });

    m1une::ds::Dsu dsu(int(points.size()));
    ManhattanMst<W> result;
    result.cost = W(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 += edge.cost;
        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_MANHATTAN_MST_HPP

#line 1 "geometry/manhattan_mst.hpp"



#include <algorithm>
#include <cassert>
#include <concepts>
#include <limits>
#include <map>
#include <numeric>
#include <utility>
#include <vector>

#line 1 "ds/dsu/dsu.hpp"



#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"



#include <cmath>
#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 15 "geometry/manhattan_mst.hpp"

namespace m1une {
namespace geometry {

template <class T>
struct ManhattanMstEdge {
    int from;
    int to;
    T cost;
};

template <class T>
struct ManhattanMst {
    T cost;
    std::vector<ManhattanMstEdge<T>> edges;
};

// Returns O(n) edges containing a Manhattan minimum spanning tree.
template <ExactCoordinate T>
std::vector<ManhattanMstEdge<wide_type<T>>> manhattan_mst_edges(const std::vector<Point<T>>& points) {
    using W = wide_type<T>;
    assert(points.size() <= std::size_t(std::numeric_limits<int>::max()));
    int n = int(points.size());
    std::vector<Point<W>> transformed;
    transformed.reserve(points.size());
    for (const auto& point : points) {
        transformed.emplace_back(W(point.x), W(point.y));
    }

    std::vector<int> order(n);
    std::iota(order.begin(), order.end(), 0);
    std::vector<ManhattanMstEdge<W>> edges;
    edges.reserve(std::size_t(4) * points.size());

    for (int direction = 0; direction < 4; direction++) {
        std::sort(order.begin(), order.end(), [&transformed](int i, int j) {
            W first = transformed[i].x + transformed[i].y;
            W second = transformed[j].x + transformed[j].y;
            if (first != second) return first < second;
            if (transformed[i].x != transformed[j].x) {
                return transformed[i].x < transformed[j].x;
            }
            return i < j;
        });

        std::map<W, int> sweep;
        for (int i : order) {
            auto it = sweep.lower_bound(-transformed[i].y);
            while (it != sweep.end()) {
                int j = it->second;
                if (transformed[i].x - transformed[j].x < transformed[i].y - transformed[j].y) {
                    break;
                }

                W dx = W(points[i].x) - W(points[j].x);
                W dy = W(points[i].y) - W(points[j].y);
                if (dx < 0) dx = -dx;
                if (dy < 0) dy = -dy;
                edges.push_back(ManhattanMstEdge<W>{i, j, dx + dy});
                it = sweep.erase(it);
            }
            sweep[-transformed[i].y] = i;
        }

        for (auto& point : transformed) {
            if (direction & 1) {
                point.x = -point.x;
            } else {
                std::swap(point.x, point.y);
            }
        }
    }
    return edges;
}

// Returns a Manhattan minimum spanning tree.
template <ExactCoordinate T>
ManhattanMst<wide_type<T>> manhattan_mst(const std::vector<Point<T>>& points) {
    using W = wide_type<T>;
    auto candidates = manhattan_mst_edges(points);
    std::sort(candidates.begin(), candidates.end(), [](const auto& a, const auto& b) { return a.cost < b.cost; });

    m1une::ds::Dsu dsu(int(points.size()));
    ManhattanMst<W> result;
    result.cost = W(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 += edge.cost;
        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

Back to top page