Closest Pair of Points

The Implementation

This follows the derivation in the parent part directly: presort by $x$, recurse on halves, and combine with a strip check that only ever compares each point against a small constant number of neighbors in $y$-order.

import math


def dist(a, b):
    return math.hypot(a[0] - b[0], a[1] - b[1])


def brute_force(points):
    best = float("inf")
    best_pair = None
    for i in range(len(points)):
        for j in range(i + 1, len(points)):
            d = dist(points[i], points[j])
            if d < best:
                best, best_pair = d, (points[i], points[j])
    return best, best_pair


def closest_pair(points_sorted_by_x):
    n = len(points_sorted_by_x)
    if n <= 3:
        return brute_force(points_sorted_by_x)

    mid = n // 2
    mid_x = points_sorted_by_x[mid][0]
    left, right = points_sorted_by_x[:mid], points_sorted_by_x[mid:]

    d_left, pair_left = closest_pair(left)
    d_right, pair_right = closest_pair(right)
    delta, best_pair = (d_left, pair_left) if d_left < d_right else (d_right, pair_right)

    # Strip: points within delta of the dividing line, sorted by y.
    strip = sorted((p for p in points_sorted_by_x if abs(p[0] - mid_x) < delta), key=lambda p: p[1])

    for i in range(len(strip)):
        j = i + 1
        # Packing argument: at most a small constant number of candidates
        # can ever be closer than delta -- 7 is a safe, common slack bound.
        while j < len(strip) and j < i + 8:
            d = dist(strip[i], strip[j])
            if d < delta:
                delta, best_pair = d, (strip[i], strip[j])
            j += 1

    return delta, best_pair


pts = [(2, 3), (12, 30), (40, 50), (5, 1), (12, 10), (3, 4)]
print(closest_pair(sorted(pts)))

#include <algorithm>
#include <cmath>
#include <vector>

struct Point { double x, y; };

double dist(const Point& a, const Point& b) {
    return std::hypot(a.x - b.x, a.y - b.y);
}

double bruteForce(const std::vector<Point>& pts) {
    double best = 1e18;
    for (size_t i = 0; i < pts.size(); ++i)
        for (size_t j = i + 1; j < pts.size(); ++j)
            best = std::min(best, dist(pts[i], pts[j]));
    return best;
}

double closestPair(std::vector<Point> pts) {  // pts must be pre-sorted by x
    if (pts.size() <= 3) return bruteForce(pts);

    size_t mid = pts.size() / 2;
    double midX = pts[mid].x;
    std::vector<Point> left(pts.begin(), pts.begin() + mid);
    std::vector<Point> right(pts.begin() + mid, pts.end());

    double delta = std::min(closestPair(left), closestPair(right));

    std::vector<Point> strip;
    for (const auto& p : pts)
        if (std::abs(p.x - midX) < delta) strip.push_back(p);
    std::sort(strip.begin(), strip.end(),
              [](const Point& a, const Point& b) { return a.y < b.y; });

    for (size_t i = 0; i < strip.size(); ++i)
        for (size_t j = i + 1; j < strip.size() && j < i + 8; ++j)
            delta = std::min(delta, dist(strip[i], strip[j]));

    return delta;
}

Common mistake: re-sorting the strip (or the whole point set) by $y$ inside every recursive call instead of maintaining a presorted-by-$y$ list across the recursion. It still produces a correct answer, but costs an extra $O(\log n)$ factor — exactly the gap between $O(n \log n)$ and $O(n \log^2 n)$ described in the parent part.

Worked Example

Trace It Yourself

Run the algorithm on [(2,3), (12,30), (40,50), (5,1), (12,10), (3,4)]. The brute-force base case fires almost immediately here since $n$ is small — trace what happens with a larger synthetic input (e.g. 20 random points) to see the strip check actually engage, and count how many comparisons the strip step performs versus $\binom{n}{2}$ for brute force.

Cookie Consent

Closest Pair of Points

The Implementation

Trace It Yourself

Related Reading

Closest Pair and Geometric Divide and Conquer