Geometric Divide-and-Conquer: Planar Closest Pair of Points in O(N log N)

The Naive Bottleneck vs. Shamos-Hoey (1975) Divide-and-Conquer

Given $N$ points $P = \{p_0, p_1, \dots, p_{N-1}\}$ in the Euclidean plane $\mathbb{R}^2$, the Closest Pair of Points problem asks for the minimum Euclidean distance $\min_{i \neq j} \|p_i - p_j\|_2$. Exhaustive pairwise inspection requires $\binom{N}{2} = \Theta(N^2)$ distance calculations. Can we apply spatial Divide-and-Conquer like Merge Sort to reach $O(N \log N)$?

1. Divide: Presort all points by $X$-coordinate ($P_X$) and by $Y$-coordinate ($P_Y$) once in $O(N \log N)$ time. At each recursive step, split $P_X$ at the median index $m = \lfloor N / 2 \rfloor$ with vertical dividing line $x = x_m$ into $P_L$ and $P_R$. Partition $P_Y$ in $O(N)$ time into $Y$-sorted subsets $P_{YL}$ and $P_{YR}$.
2. Conquer: Recursively find the minimum distance $\delta_L$ in $P_L$ and $\delta_R$ in $P_R$. Let $\delta = \min(\delta_L, \delta_R)$.
3. Combine (The $2\delta$ Strip): Any pair $(p_l, p_r)$ with $p_l \in P_L, p_r \in P_R$ that beats $\delta$ must have both points within horizontal distance $\delta$ of the dividing line $x = x_m$. Filter $P_Y$ in $O(N)$ time (preserving $Y$-sorted order) into a strip array $S = \{p \in P_Y : |p.x - x_m| < \delta\}$.

Production C++20 Implementation: $O(N \log N)$ Closest Pair of Points

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <vector>

struct Point2D {
    double x;
    double y;
    std::uint32_t id;
};

struct ClosestPairResult {
    double min_dist;
    std::uint32_t id1;
    std::uint32_t id2;
};

inline double dist_sq(const Point2D& a, const Point2D& b) noexcept {
    const double dx = a.x - b.x;
    const double dy = a.y - b.y;
    return dx * dx + dy * dy;
}

ClosestPairResult solve_closest_pair(std::vector<Point2D>& pts) {
    const std::size_t n = pts.size();
    if (n < 2) return {std::numeric_limits<double>::infinity(), 0, 0};

    // 1. Sort initially by X-coordinate (tie-break by Y)
    std::sort(pts.begin(), pts.end(), [](const Point2D& a, const Point2D& b) {
        return (a.x != b.x) ? (a.x < b.x) : (a.y < b.y);
    });

    std::vector<Point2D> scratch(n);
    std::vector<Point2D> strip(n);
    double best_d2 = std::numeric_limits<double>::infinity();
    std::uint32_t best_u = 0, best_v = 0;

    auto update_best = [&](const Point2D& a, const Point2D& b) noexcept {
        const double d2 = dist_sq(a, b);
        if (d2 < best_d2) {
            best_d2 = d2;
            best_u = a.id;
            best_v = b.id;
        }
    };

    auto rec = [&](auto& self, std::size_t lo, std::size_t hi) -> void {
        if (hi - lo <= 8) {
            // Small base case: O(1) brute force + insertion sort by Y
            for (std::size_t i = lo; i < hi; ++i) {
                for (std::size_t j = i + 1; j < hi; ++j) {
                    update_best(pts[i], pts[j]);
                }
            }
            std::sort(pts.begin() + lo, pts.begin() + hi,
                      [](const Point2D& a, const Point2D& b) { return a.y < b.y; });
            return;
        }

        const std::size_t mid = lo + (hi - lo) / 2;
        const double mid_x = pts[mid].x;

        self(self, lo, mid);
        self(self, mid, hi);

        // Merge [lo, mid) and [mid, hi) by Y-coordinate in O(N) time
        std::merge(pts.begin() + lo, pts.begin() + mid,
                   pts.begin() + mid, pts.begin() + hi,
                   scratch.begin() + lo,
                   [](const Point2D& a, const Point2D& b) { return a.y < b.y; });
        std::copy(scratch.begin() + lo, scratch.begin() + hi, pts.begin() + lo);

        // Collect points within horizontal distance sqrt(best_d2) of mid_x
        std::size_t s_len = 0;
        for (std::size_t i = lo; i < hi; ++i) {
            const double dx = pts[i].x - mid_x;
            if (dx * dx < best_d2) {
                strip[s_len++] = pts[i];
            }
        }

        // Scan Y-sorted strip: inner loop executes at most 7 times per point!
        for (std::size_t i = 0; i < s_len; ++i) {
            for (std::size_t j = i + 1; j < s_len; ++j) {
                const double dy = strip[j].y - strip[i].y;
                if (dy * dy >= best_d2) break; // Early exit on vertical distance
                update_best(strip[i], strip[j]);
            }
        }
    };

    rec(rec, 0, n);
    return {std::sqrt(best_d2), best_u, best_v};
}