Closest Pair of Points — Why Your Strip Takes O(n²)

The closest pair problem is a beautiful example of computational geometry enabling an otherwise-quadratic combine step to run in linear time. The naive O(n^2) solution checks every pair. The D&C solution achieves O(n log n) by proving a remarkable geometric fact: in the critical 'strip' around the dividing line, each point has at most 7 other points it needs to check.

This seven-point bound — derived from packing arguments in a 2δ×δ rectangle — is the kind of geometric reasoning that separates algorithm designers from algorithm users. It is also a template for many computational geometry optimisations: use the structure of space to bound the number of interactions.

Algorithm Overview

1. Sort points by x-coordinate: O(n log n)
2. Recursively find closest pair in left half (δ_L) and right half (δ_R)
3. δ = min(δ_L, δ_R)
4. Find all points within δ of the dividing line (the strip)
5. Check pairs in the strip sorted by y-coordinate
6. The key: each point in the strip has at most 7 other points to check

# Package: io.thecodeforge.closest_pair
from math import dist
from itertools import combinations

def closest_pair(points: list[tuple]) -> tuple[float, tuple, tuple]:
    points = sorted(points, key=lambda p: p[0])
    return _closest(points)

def _closest(pts):
    n = len(pts)
    if n <= 3:
        # Brute force for small input
        best = float('inf'), pts[0], pts[1]
        for a, b in combinations(pts, 2):
            d = dist(a, b)
            if d < best[0]: best = d, a, b
        return best
    mid = n // 2
    mid_x = pts[mid][0]
    d_left  = _closest(pts[:mid])
    d_right = _closest(pts[mid:])
    best = d_left if d_left[0] < d_right[0] else d_right
    delta = best[0]
    # Build strip: points within delta of dividing line
    strip = [p for p in pts if abs(p[0] - mid_x) < delta]
    strip.sort(key=lambda p: p[1])  # sort by y
    # Check strip pairs (at most 7 comparisons per point)
    for i in range(len(strip)):
        for j in range(i+1, min(i+8, len(strip))):
            if strip[j][1] - strip[i][1] >= delta:
                break
            d = dist(strip[i], strip[j])
            if d < delta:
                best = d, strip[i], strip[j]
                delta = d
    return best

points = [(2,3),(12,30),(40,50),(5,1),(12,10),(3,4)]
d, p1, p2 = closest_pair(points)
print(f'Closest: {p1} and {p2}, distance: {d:.4f}')

The 7-Neighbour Proof

Why does each point in the strip have at most 7 other points to check? Consider a δ×2δ rectangle straddling the midline. In each δ×δ half, the closest pair distance is ≥ δ (we already found the best in each half). So in each δ×δ quarter, we can pack at most 4 points at distance ≥ δ from each other (a 2×2 grid). Total: 8 points maximum in the 2δ×2δ strip rectangle including the reference point itself — so at most 7 others.

This bounds the inner loop to constant work per point → O(n) combine step.

Optimisation — Pre-Sort by Y

Re-sorting the strip each time adds a log factor. The optimisation: maintain two sorted arrays throughout — one sorted by x (for dividing) and one sorted by y (for strip checking). This gives T(n) = 2T(n/2) + O(n) → O(n log n) with smaller constants.

Implementation detail: during recursion, merge the y-sorted lists of the left and right halves (like merge sort) to produce the parent's y-sorted list. The strip can then be filtered in O(|strip|) from the y-sorted list without a full sort.

Complexity Analysis and Recurrence

The recurrence for the naive version (no pre-sort): T(n) = 2T(n/2) + O(n log n). Solving via Master Theorem: a=2, b=2, f(n)=O(n log n) = O(n^(log₂2) log n) = O(n log n). This falls under Case 2 of the Master Theorem (actually a special case). The solution is O(n log² n). Because the combine step sorts the strip, it's O(n log n), leading to T(n) = 2T(n/2) + O(n log n) = O(n log² n). Wait — careful: if we sort the strip each time, the combine is O(k log k) where k ≤ n. That sums to O(n log² n) in the worst case. Pre-sorting brings combine to O(n) and recurrence to O(n log n).

So pre-sorting doesn't just reduce constants — it changes the asymptotic complexity.

Implementation Details and Edge Cases

Base Case: n ≤ 3. Brute force is efficient and avoids recursion overhead. Use combinations for clarity, but a triple-nested loop is fine for n=3.

Median Selection: Sort by x once, then the median is simply p[mid]. No need for linear-time selection — the sort dominates anyway.

Duplicate Points: If two points have identical coordinates, the distance is 0. The algorithm should handle this correctly because δ becomes 0 and the strip width becomes 0 — only points exactly on the midline are considered, but duplicates on the same side will be caught in the recursive step. However, if duplicates exist across halves, they may be missed because the strip is empty (δ=0). Fix: before divide, deduplicate or handle zero-distance case explicitly.

Vertical and Horizontal Lines: Sorting by x works fine. If many points share the same x, the median may be ambiguous. Use a stable sort or choose the middle index after sorting.

Floating-Point Precision: Distance calculations can suffer from floating-point errors. Use math.dist which is stable, but for comparison, prefer squared distances to avoid sqrt overhead.

Integer Coordinates: If all coordinates are integers, use squared distances to avoid sqrt entirely, only taking sqrt at the end.

# Package: io.thecodeforge.closest_pair_optimized
from math import dist

def closest_pair_with_pre_sort(points: list[tuple]) -> tuple[float, tuple, tuple]:
    # Pre-sort by x for divide step
    px = sorted(points, key=lambda p: p[0])
    py = sorted(points, key=lambda p: p[1])
    return _closest(px, py)

def _closest(px: list, py: list) -> tuple[float, tuple, tuple]:
    n = len(px)
    if n <= 3:
        # Brute force
        best = float('inf'), px[0], px[1]
        for i in range(n):
            for j in range(i+1, n):
                d = dist(px[i], px[j])
                if d < best[0]:
                    best = d, px[i], px[j]
        return best
    mid = n // 2
    mid_x = px[mid][0]
    # Split y-sorted points into left and right based on x
    py_left = []
    py_right = []
    for p in py:
        if p[0] <= mid_x and len(py_left) < mid:
            py_left.append(p)
        else:
            py_right.append(p)
    # Recurse
    d_left = _closest(px[:mid], py_left)
    d_right = _closest(px[mid:], py_right)
    delta = min(d_left[0], d_right[0])
    # Build strip from py (already sorted by y)
    strip = [p for p in py if abs(p[0] - mid_x) < delta]
    # Check up to 7 neighbours
    best = d_left if d_left[0] < d_right[0] else d_right
    for i in range(len(strip)):
        for j in range(i+1, min(i+8, len(strip))):
            if strip[j][1] - strip[i][1] >= delta:
                break
            d = dist(strip[i], strip[j])
            if d < best[0]:
                best = d, strip[i], strip[j]
                delta = d
    return best

Extensions and Higher Dimensions

The D&C approach generalises to k dimensions. The strip becomes a slab of width 2δ, and the neighbour bound becomes O(3^k), which grows exponentially with dimension. For 3D, each point checks at most 31 neighbours (3^3 - 1 = 26? Actually the bound is more complex). But the combinatorial explosion makes D&C impractical for dimensions > 3. For high-dimensional data, use space-filling curves or KD-trees.

Related Problems: The packing argument also applies to finding the minimal distance between two convex polygons, or computing the diameter of a point set using rotating calipers.

Practical Note: In production, for 2D and 3D, the D&C algorithm is the gold standard. For higher dimensions, approximate methods (like LSH) are used.