Closest Pair of Points
(geometry/closest_pair.hpp)
- View this file on GitHub
- Last update: 2026-10-05 22:23:07+09:00
- Include:
#include "geometry/closest_pair.hpp"
Overview
closest_pair(points) returns two distinct input points whose Euclidean
distance is minimum. It uses divide and conquer after sorting the points by
their x-coordinate.
Duplicate points are supported and produce squared distance zero.
Interface
template <Coordinate T>
struct ClosestPair {
int first;
int second;
wide_type<T> distance_squared;
};
template <Coordinate T>
std::optional<ClosestPair<T>> closest_pair(
const std::vector<Point<T>>& points
);
Result
| Member | Description |
|---|---|
first |
Smaller original index of the selected pair. |
second |
Larger original index of the selected pair. |
distance_squared |
Squared Euclidean distance in wide_type<T>. |
The function returns std::nullopt when fewer than two points are supplied.
When several pairs have the same minimum distance, the lexicographically
smallest index pair is returned.
For integral T, squared distances are calculated with signed 128-bit
arithmetic. They must fit in wide_type<T>. Floating-point coordinates are
compared without an epsilon.
Operations
| Function | Description | Complexity |
|---|---|---|
closest_pair(const std::vector<Point<T>>& points) |
Finds the closest pair without modifying points. |
$O(N\log N)$ time and $O(N)$ memory |
Example
#include "geometry/closest_pair.hpp"
#include <cassert>
#include <vector>
int main() {
using Point = m1une::geometry::Point<long long>;
std::vector<Point> points;
points.emplace_back(0, 0);
points.emplace_back(5, 0);
points.emplace_back(2, 1);
auto answer = m1une::geometry::closest_pair(points);
assert(answer.has_value());
assert(answer->first == 0);
assert(answer->second == 2);
assert(answer->distance_squared == 5);
}
Depends on
Required by
Verified with
- verify/geometry/centroid.test.cpp
- verify/geometry/closest_pair.test.cpp
- verify/geometry/geometry_algorithms.test.cpp
- verify/geometry/rational.test.cpp
Code
#ifndef M1UNE_GEOMETRY_CLOSEST_PAIR_HPP
#define M1UNE_GEOMETRY_CLOSEST_PAIR_HPP 1
#include <algorithm>
#include <optional>
#include <utility>
#include <vector>
#include "point.hpp"
namespace m1une {
namespace geometry {
template <Coordinate T>
struct ClosestPair {
int first;
int second;
wide_type<T> distance_squared;
};
// Returns two distinct original indices with minimum Euclidean distance.
template <Coordinate T>
std::optional<ClosestPair<T>> closest_pair(
const std::vector<Point<T>>& points
) {
const int n = int(points.size());
if (n < 2) return std::nullopt;
struct IndexedPoint {
Point<T> point;
int index;
};
std::vector<IndexedPoint> ordered;
ordered.reserve(n);
for (int index = 0; index < n; index++) {
ordered.push_back(IndexedPoint{points[index], index});
}
std::sort(
ordered.begin(),
ordered.end(),
[](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point < second.point) return true;
if (second.point < first.point) return false;
return first.index < second.index;
}
);
std::optional<ClosestPair<T>> duplicate_result;
for (int first = 0; first < n;) {
int last = first + 1;
while (last < n && ordered[last].point == ordered[first].point) last++;
if (last - first >= 2) {
int first_index = ordered[first].index;
int second_index = ordered[first + 1].index;
std::pair<int, int> candidate(first_index, second_index);
if (
!duplicate_result ||
candidate < std::pair(
duplicate_result->first,
duplicate_result->second
)
) {
duplicate_result = ClosestPair<T>{
first_index,
second_index,
wide_type<T>(0)
};
}
}
first = last;
}
if (duplicate_result) return duplicate_result;
std::optional<ClosestPair<T>> result;
auto consider = [&result, &points](int first, int second) {
if (second < first) std::swap(first, second);
wide_type<T> squared = distance2(points[first], points[second]);
std::pair<int, int> candidate(first, second);
if (
!result ||
squared < result->distance_squared ||
(
squared == result->distance_squared &&
candidate < std::pair(result->first, result->second)
)
) {
result = ClosestPair<T>{first, second, squared};
}
};
consider(ordered[0].index, ordered[1].index);
auto by_y = [](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point.y != second.point.y) {
return first.point.y < second.point.y;
}
if (first.point.x != second.point.x) {
return first.point.x < second.point.x;
}
return first.index < second.index;
};
std::vector<IndexedPoint> buffer(n);
auto solve = [&](auto&& self, int left, int right) -> void {
if (right - left <= 3) {
for (int first = left; first < right; first++) {
for (int second = first + 1; second < right; second++) {
consider(ordered[first].index, ordered[second].index);
}
}
std::sort(ordered.begin() + left, ordered.begin() + right, by_y);
return;
}
int middle = (left + right) / 2;
T middle_x = ordered[middle].point.x;
self(self, left, middle);
self(self, middle, right);
std::merge(
ordered.begin() + left,
ordered.begin() + middle,
ordered.begin() + middle,
ordered.begin() + right,
buffer.begin() + left,
by_y
);
std::copy(
buffer.begin() + left,
buffer.begin() + right,
ordered.begin() + left
);
int strip_size = 0;
for (int index = left; index < right; index++) {
using W = wide_type<T>;
W dx = W(ordered[index].point.x) - W(middle_x);
if (result->distance_squared < dx * dx) continue;
for (int previous = strip_size - 1; previous >= 0; previous--) {
W dy = W(ordered[index].point.y) - W(buffer[previous].point.y);
if (result->distance_squared < dy * dy) break;
consider(ordered[index].index, buffer[previous].index);
}
buffer[strip_size++] = ordered[index];
}
};
solve(solve, 0, n);
return result;
}
} // namespace geometry
} // namespace m1une
#endif // M1UNE_GEOMETRY_CLOSEST_PAIR_HPP#line 1 "geometry/closest_pair.hpp"
#include <algorithm>
#include <optional>
#include <utility>
#include <vector>
#line 1 "geometry/point.hpp"
#include <cmath>
#include <concepts>
#include <cassert>
#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 10 "geometry/closest_pair.hpp"
namespace m1une {
namespace geometry {
template <Coordinate T>
struct ClosestPair {
int first;
int second;
wide_type<T> distance_squared;
};
// Returns two distinct original indices with minimum Euclidean distance.
template <Coordinate T>
std::optional<ClosestPair<T>> closest_pair(
const std::vector<Point<T>>& points
) {
const int n = int(points.size());
if (n < 2) return std::nullopt;
struct IndexedPoint {
Point<T> point;
int index;
};
std::vector<IndexedPoint> ordered;
ordered.reserve(n);
for (int index = 0; index < n; index++) {
ordered.push_back(IndexedPoint{points[index], index});
}
std::sort(
ordered.begin(),
ordered.end(),
[](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point < second.point) return true;
if (second.point < first.point) return false;
return first.index < second.index;
}
);
std::optional<ClosestPair<T>> duplicate_result;
for (int first = 0; first < n;) {
int last = first + 1;
while (last < n && ordered[last].point == ordered[first].point) last++;
if (last - first >= 2) {
int first_index = ordered[first].index;
int second_index = ordered[first + 1].index;
std::pair<int, int> candidate(first_index, second_index);
if (
!duplicate_result ||
candidate < std::pair(
duplicate_result->first,
duplicate_result->second
)
) {
duplicate_result = ClosestPair<T>{
first_index,
second_index,
wide_type<T>(0)
};
}
}
first = last;
}
if (duplicate_result) return duplicate_result;
std::optional<ClosestPair<T>> result;
auto consider = [&result, &points](int first, int second) {
if (second < first) std::swap(first, second);
wide_type<T> squared = distance2(points[first], points[second]);
std::pair<int, int> candidate(first, second);
if (
!result ||
squared < result->distance_squared ||
(
squared == result->distance_squared &&
candidate < std::pair(result->first, result->second)
)
) {
result = ClosestPair<T>{first, second, squared};
}
};
consider(ordered[0].index, ordered[1].index);
auto by_y = [](const IndexedPoint& first, const IndexedPoint& second) {
if (first.point.y != second.point.y) {
return first.point.y < second.point.y;
}
if (first.point.x != second.point.x) {
return first.point.x < second.point.x;
}
return first.index < second.index;
};
std::vector<IndexedPoint> buffer(n);
auto solve = [&](auto&& self, int left, int right) -> void {
if (right - left <= 3) {
for (int first = left; first < right; first++) {
for (int second = first + 1; second < right; second++) {
consider(ordered[first].index, ordered[second].index);
}
}
std::sort(ordered.begin() + left, ordered.begin() + right, by_y);
return;
}
int middle = (left + right) / 2;
T middle_x = ordered[middle].point.x;
self(self, left, middle);
self(self, middle, right);
std::merge(
ordered.begin() + left,
ordered.begin() + middle,
ordered.begin() + middle,
ordered.begin() + right,
buffer.begin() + left,
by_y
);
std::copy(
buffer.begin() + left,
buffer.begin() + right,
ordered.begin() + left
);
int strip_size = 0;
for (int index = left; index < right; index++) {
using W = wide_type<T>;
W dx = W(ordered[index].point.x) - W(middle_x);
if (result->distance_squared < dx * dx) continue;
for (int previous = strip_size - 1; previous >= 0; previous--) {
W dy = W(ordered[index].point.y) - W(buffer[previous].point.y);
if (result->distance_squared < dy * dy) break;
consider(ordered[index].index, buffer[previous].index);
}
buffer[strip_size++] = ordered[index];
}
};
solve(solve, 0, n);
return result;
}
} // namespace geometry
} // namespace m1une