Farthest Pair of Points
(geometry/farthest_pair.hpp)

• View this file on GitHub
• Last update: 2026-10-05 22:23:07+09:00
• Include: #include "geometry/farthest_pair.hpp"

Overview

farthest_pair(points) returns two distinct input points whose Euclidean distance is maximum. It first constructs the convex hull, then examines antipodal hull vertices with rotating calipers.

furthest_pair(points) is an alias with identical behavior.

Result

template <Coordinate T>
struct FarthestPair {
    int first;
    int second;
    wide_type<T> distance_squared;
};

first and second are indices in the original input, ordered so that first < second. The distance is squared, avoiding a square root and retaining exact arithmetic for integral coordinates.

Both functions return std::nullopt when fewer than two points are supplied. Duplicate points are supported. If all points coincide, two distinct input indices are returned with squared distance zero.

When several pairs have the same maximum distance, the lexicographically smallest index pair is returned.

Complexity

The time complexity is $O(N\log N)$ and the memory usage is $O(N)$.

Example

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

    auto result = m1une::geometry::farthest_pair(points);
    if (result) {
        std::cout << result->first << ' ' << result->second << "\n";
        std::cout << static_cast<long long>(result->distance_squared) << "\n";
    }
}

Depends on

• Convex Hull (geometry/convex_hull.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/farthest_pair.test.cpp
• verify/geometry/geometry_algorithms.test.cpp
• verify/geometry/rational.test.cpp

Code

#ifndef M1UNE_GEOMETRY_FARTHEST_PAIR_HPP
#define M1UNE_GEOMETRY_FARTHEST_PAIR_HPP 1

#include <algorithm>
#include <cstddef>
#include <map>
#include <optional>
#include <utility>
#include <vector>

#include "convex_hull.hpp"

namespace m1une {
namespace geometry {

template <Coordinate T>
struct FarthestPair {
    int first;
    int second;
    wide_type<T> distance_squared;
};

// Returns two distinct original indices with maximum Euclidean distance.
template <Coordinate T>
std::optional<FarthestPair<T>> farthest_pair(
    const std::vector<Point<T>>& points
) {
    if (points.size() < 2) return std::nullopt;

    std::vector<Point<T>> hull = convex_hull(points);
    if (hull.size() == 1) {
        FarthestPair<T> result;
        result.first = 0;
        result.second = 1;
        result.distance_squared = 0;
        return result;
    }

    std::map<Point<T>, int> original_index;
    for (int index = 0; index < int(points.size()); ++index) {
        original_index.emplace(points[index], index);
    }

    std::vector<int> hull_index;
    hull_index.reserve(hull.size());
    for (const Point<T>& point : hull) {
        hull_index.push_back(original_index.find(point)->second);
    }

    std::optional<FarthestPair<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]);
        if (
            !result.has_value() ||
            result->distance_squared < squared ||
            (
                result->distance_squared == squared &&
                std::pair(first, second) <
                    std::pair(result->first, result->second)
            )
        ) {
            result = FarthestPair<T>{first, second, squared};
        }
    };

    if (hull.size() == 2) {
        consider(hull_index[0], hull_index[1]);
        return result;
    }

    std::size_t opposite = 1;
    for (std::size_t index = 0; index < hull.size(); ++index) {
        std::size_t next = (index + 1) % hull.size();
        while (true) {
            std::size_t candidate = (opposite + 1) % hull.size();
            auto current_area = cross(
                hull[index],
                hull[next],
                hull[opposite]
            );
            auto candidate_area = cross(
                hull[index],
                hull[next],
                hull[candidate]
            );
            if (candidate_area <= current_area) break;
            opposite = candidate;
        }
        consider(hull_index[index], hull_index[opposite]);
        consider(hull_index[next], hull_index[opposite]);

        std::size_t candidate = (opposite + 1) % hull.size();
        auto current_area = cross(hull[index], hull[next], hull[opposite]);
        auto candidate_area = cross(hull[index], hull[next], hull[candidate]);
        if (candidate_area == current_area) {
            consider(hull_index[index], hull_index[candidate]);
            consider(hull_index[next], hull_index[candidate]);
        }
    }
    return result;
}

template <Coordinate T>
std::optional<FarthestPair<T>> furthest_pair(
    const std::vector<Point<T>>& points
) {
    return farthest_pair(points);
}

}  // namespace geometry
}  // namespace m1une

#endif  // M1UNE_GEOMETRY_FARTHEST_PAIR_HPP

#line 1 "geometry/farthest_pair.hpp"



#include <algorithm>
#include <cstddef>
#include <map>
#include <optional>
#include <utility>
#include <vector>

#line 1 "geometry/convex_hull.hpp"



#line 8 "geometry/convex_hull.hpp"

#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/convex_hull.hpp"

namespace m1une {
namespace geometry {

// Returns the convex hull counterclockwise from its lexicographically smallest
// point. The first point is not repeated at the end.
template <Coordinate T>
std::vector<Point<T>> convex_hull(
    std::vector<Point<T>> points,
    bool include_collinear = false
) {
    std::sort(points.begin(), points.end());
    points.erase(std::unique(points.begin(), points.end()), points.end());
    std::size_t size = points.size();
    if (size <= 1) return points;

    std::vector<Point<T>> hull;
    hull.reserve(2 * size);
    auto should_pop = [include_collinear](
        const Point<T>& first,
        const Point<T>& second,
        const Point<T>& third
    ) {
        int turn = orientation(first, second, third);
        return include_collinear ? turn < 0 : turn <= 0;
    };

    for (const Point<T>& point : points) {
        while (
            hull.size() >= 2 &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }

    std::size_t lower_size = hull.size();
    for (std::size_t index = size - 1; index-- > 0;) {
        const Point<T>& point = points[index];
        while (
            hull.size() > lower_size &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }
    hull.pop_back();

    if (include_collinear && hull.size() == 2 * points.size() - 2) {
        hull = std::move(points);
    }
    return hull;
}

}  // namespace geometry
}  // namespace m1une


#line 12 "geometry/farthest_pair.hpp"

namespace m1une {
namespace geometry {

template <Coordinate T>
struct FarthestPair {
    int first;
    int second;
    wide_type<T> distance_squared;
};

// Returns two distinct original indices with maximum Euclidean distance.
template <Coordinate T>
std::optional<FarthestPair<T>> farthest_pair(
    const std::vector<Point<T>>& points
) {
    if (points.size() < 2) return std::nullopt;

    std::vector<Point<T>> hull = convex_hull(points);
    if (hull.size() == 1) {
        FarthestPair<T> result;
        result.first = 0;
        result.second = 1;
        result.distance_squared = 0;
        return result;
    }

    std::map<Point<T>, int> original_index;
    for (int index = 0; index < int(points.size()); ++index) {
        original_index.emplace(points[index], index);
    }

    std::vector<int> hull_index;
    hull_index.reserve(hull.size());
    for (const Point<T>& point : hull) {
        hull_index.push_back(original_index.find(point)->second);
    }

    std::optional<FarthestPair<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]);
        if (
            !result.has_value() ||
            result->distance_squared < squared ||
            (
                result->distance_squared == squared &&
                std::pair(first, second) <
                    std::pair(result->first, result->second)
            )
        ) {
            result = FarthestPair<T>{first, second, squared};
        }
    };

    if (hull.size() == 2) {
        consider(hull_index[0], hull_index[1]);
        return result;
    }

    std::size_t opposite = 1;
    for (std::size_t index = 0; index < hull.size(); ++index) {
        std::size_t next = (index + 1) % hull.size();
        while (true) {
            std::size_t candidate = (opposite + 1) % hull.size();
            auto current_area = cross(
                hull[index],
                hull[next],
                hull[opposite]
            );
            auto candidate_area = cross(
                hull[index],
                hull[next],
                hull[candidate]
            );
            if (candidate_area <= current_area) break;
            opposite = candidate;
        }
        consider(hull_index[index], hull_index[opposite]);
        consider(hull_index[next], hull_index[opposite]);

        std::size_t candidate = (opposite + 1) % hull.size();
        auto current_area = cross(hull[index], hull[next], hull[opposite]);
        auto candidate_area = cross(hull[index], hull[next], hull[candidate]);
        if (candidate_area == current_area) {
            consider(hull_index[index], hull_index[candidate]);
            consider(hull_index[next], hull_index[candidate]);
        }
    }
    return result;
}

template <Coordinate T>
std::optional<FarthestPair<T>> furthest_pair(
    const std::vector<Point<T>>& points
) {
    return farthest_pair(points);
}

}  // namespace geometry
}  // namespace m1une

Back to top page