xvi. Closest Pair

Closest Pair

Find the two closest of $n$ points. Naive is $O(n^2)$; three better methods exist.

Method 1: sweep line with a set — the one to write

$O(n\log n)$, about 20 lines, and it reuses std::set.

ll closestPair(vector<P> p) {
    sort(p.begin(), p.end());                          // by x, then y
    set<pair<ll, ll>> box;                              // active points, keyed by (y, x)
    ll best = LLONG_MAX;
    int left = 0;
    for (int i = 0; i < (int)p.size(); i++) {
        ll d = (ll)ceil(sqrt((long double)best));
        while (left < i && p[i].x - p[left].x > d) {    // drop points too far left
            box.erase({p[left].y, p[left].x});
            left++;
        }
        auto lo = box.lower_bound({p[i].y - d, LLONG_MIN});
        auto hi = box.upper_bound({p[i].y + d, LLONG_MAX});
        for (auto it = lo; it != hi; ++it) best = min(best, norm2(P{it->second, it->first} - p[i]));
        box.insert({p[i].y, p[i].x});
    }
    return best;                                        // SQUARED distance
}

Why it is $O(n\log n)$: the active set is a vertical strip of width $d$, and within any $d\times 2d$ window at most 6 points can be pairwise at distance $\ge d$ (the same packing bound proved below for the divide-and-conquer strip). So each point examines $O(1)$ neighbours, and d only shrinks, so the strip never widens.

Returning the squared distance keeps everything exact.

Proof: closestPair's correctness

Invariant. After processing index $i$, best is the minimum squared distance among all pairs within $p_0,\dots ,p_i$, and box holds exactly the points $p_{left},\dots ,p_i$ (those not yet provably too far in $x$ to matter) keyed by $(y,x)$.
For the new point $p_i$, any earlier point closer than $\sqrt{\text{best}}$ must have both $|\Delta x|<d=\sqrt{\text{best}}$ (hence still in box, not yet dropped by the while loop, since points are processed in $x$-order) and $|\Delta y|<d$ (hence in the queried $y$-range $[p_i.y-d,p_i.y+d]$) — so the lower_bound/upper_bound range contains every candidate that could possibly improve best, and nothing closer is missed. Every candidate examined has its exact squared distance computed via norm2, so best only ever holds a true achieved distance, never an underestimate — giving a tight final answer when the loop completes. Dropping points with $p_i.x-p_{left}.x>d$ is safe because such a point is already farther than the current best in $x$ alone, so it can never be part of a future closer pair (as $d$ only shrinks over time, this stays true for the rest of the sweep). $■$ (The $O(1)$-candidates-per-query bound is the packing argument proved for the strip below, applied to a $d\times 2d$ window instead of $d\times d$.)

Method 2: divide and conquer

Split at the median $x$, recurse, then check the strip of width $2d$ around the split line, sorted by $y$. The seven-point lemma bounds the strip scan at $O(n)$:
$T(n)=2T(n/2)+O(n)=O(n\log n).$
Merge the $y$-orders during the recursion rather than re-sorting — see Bentley-Shamos below for the full implementation.

Method 3: randomised incremental — expected $O(n)$

1. Shuffle the points.
2. Maintain the best distance $d$ and a hash grid with cell size $d$.
3. Insert points one at a time, checking only the 9 neighbouring cells.
4. When $d$ shrinks, rebuild the grid.

ll closestRandomized(vector<P> p) {
    shuffle(p.begin(), p.end(), rng);
    ll best = norm2(p[1] - p[0]);
    unordered_map<ll, vector<int>> grid;
    ll cell = max(1LL, (ll)sqrt((long double)best));
    auto key = [&](ll x, ll y) { return x * 1000003LL + y; };
    auto rebuild = [&](int upto) {
        grid.clear();
        cell = max(1LL, (ll)sqrt((long double)best));
        for (int i = 0; i < upto; i++) grid[key(p[i].x / cell, p[i].y / cell)].push_back(i);
    };
    rebuild(2);
    for (int i = 2; i < (int)p.size(); i++) {
        ll gx = p[i].x / cell, gy = p[i].y / cell;
        bool shrunk = false;
        for (ll dx = -1; dx <= 1; dx++)
            for (ll dy = -1; dy <= 1; dy++)
                for (int j : grid[key(gx + dx, gy + dy)]) {
                    ll dd = norm2(p[i] - p[j]);
                    if (dd < best) { best = dd; shrunk = true; }
                }
        if (shrunk) rebuild(i + 1);
        else grid[key(gx, gy)].push_back(i);
    }
    return best;
}

Proof: closestRandomized is correct, and expected $O(n)$

Correctness. The grid’s cell size is always $\sqrt{\text{best}}$ (current), so two points within best of each other must lie in the same cell or adjacent cells (their coordinate difference in each axis is less than the cell size, so their cell indices differ by at most 1) — exactly the $3\times 3$ neighbourhood checked. Every pair actually compared uses the exact norm2, so best is always a true achieved distance; after processing all points, no closer pair exists, since every point, at the moment it was inserted, was checked against every already-inserted point close enough to matter (those in the current $3\times 3$ neighbourhood under the cell size valid at that time — and crucially, a rebuild re-buckets all previously inserted points under the new, smaller cell size immediately after best shrinks, so no stale bucketing from a now-too-large cell size persists).
Expected $O(n)$. The total cost is the sum over $i$ of $O(1)$ (the $3\times 3$ scan) plus, when a rebuild is triggered, $O(i)$. A rebuild is triggered when point $i$ improves on the best distance found among the first $i$ points — equivalently, when $p_i$ is one of the (at most 2, typically) points defining the closest pair among a random prefix of size $i$. By the same backward-analysis argument as Welzl’s algorithm, for a uniformly random insertion order this has probability $O(1/i)$, so the expected rebuild cost at step $i$ is $O(i)\cdot O(1/i)=O(1)$, and summing over $n$ steps gives $O(n)$ expected total. $■$

Comparison

Method	Time	Notes
Brute force	$O(n^2)$	fine to $n\approx 5000$
Sweep + set	$O(n\log n)$	the default
Divide and conquer	$O(n\log n)$	classic; also gives all near pairs — Bentley-Shamos
Randomised + grid	$O(n)$ expected	fastest, needs a good hash
KD-tree	$O(n\log n)$ expected	reusable for other queries — 09-kd-tree-and-spatial-structures
Delaunay	$O(n\log n)$	the closest pair is a Delaunay edge — 20-delaunay-and-voronoi

Variants

Variant	Approach
Farthest pair	it is a hull diameter — rotating calipers
Closest pair in 3D	the same sweep with a 2D active structure, $O(n{\log }^{2}n)$
Closest pair under $L_1$ / $L_{\infty }$	rotate 45° (Manhattan) or sweep directly — 17-manhattan-geometry
Closest pair with insertions	the randomised grid handles it naturally
All pairs within distance $d$	grid buckets of side $d$, then check 9 cells
$k$ closest pairs	a heap over the candidates found by the sweep
Closest pair between two sets	the same sweep, ignoring same-set pairs

The bucket trick

When the coordinate range is bounded and the answer is expected to be small, a plain grid of side $d$ answers “all pairs within $d$” in $O(n+\text{output})$. It is the workhorse for collision detection and $n$-body simulation, and often simpler than any of the above.

Tested: closestPair and closestRandomized both checked against brute-force $O(n^2)$ on random point sets (including many duplicate and collinear points), and against each other.

Bentley-Shamos

The divide and conquer approach to the closest pair of points, in $O(n\log n)$ — the model for a whole family of “divide the plane, recurse, merge with a strip” geometric algorithms.

Closest pair by divide and conquer

1. Sort the points by $x$ once, up front.
2. Split at the median $x$ into left half $L$ and right half $R$.
3. Recurse: let $d=\min (d_L,d_R)$.
4. The strip. Any closer pair must straddle the split line, so both points lie within distance $d$ of it. Collect those points, sorted by $y$.
5. The seven-point lemma. For each point in the strip, only the next at most 7 points in $y$-order can be within $d$ — because a $d\times 2d$ rectangle can hold at most 8 points pairwise $\ge d$ apart. So the strip scan is $O(n)$, not $O(n^2)$.

$T(n)=2T(n/2)+O(n)=O(n\log n).$

Keeping the $y$-order by merging during the recursion (rather than re-sorting) is what keeps the merge step linear.

pair<ll, vector<P>> closestRec(const vector<P>& px) {    // px sorted by x; returns {d^2, points sorted by y}
    int n = px.size();
    if (n <= 3) {
        ll best = LLONG_MAX;
        for (int i = 0; i < n; i++) for (int j = i + 1; j < n; j++) best = min(best, norm2(px[i] - px[j]));
        vector<P> py = px;
        sort(py.begin(), py.end(), [](P a, P b) { return a.y != b.y ? a.y < b.y : a.x < b.x; });
        return {best, py};
    }
    int mid = n / 2;
    ll midx = px[mid].x;
    vector<P> lx(px.begin(), px.begin() + mid), rx(px.begin() + mid, px.end());
    auto [dl, ly] = closestRec(lx);
    auto [dr, ry] = closestRec(rx);
    ll d = min(dl, dr);
 
    vector<P> py(n);
    merge(ly.begin(), ly.end(), ry.begin(), ry.end(), py.begin(),
          [](P a, P b) { return a.y != b.y ? a.y < b.y : a.x < b.x; });
 
    vector<P> strip;
    for (auto& p : py) if ((p.x - midx) * (p.x - midx) < d) strip.push_back(p);
    for (size_t i = 0; i < strip.size(); i++)
        for (size_t j = i + 1; j < strip.size() && (strip[j].y - strip[i].y) * (strip[j].y - strip[i].y) < d; j++)
            d = min(d, norm2(strip[i] - strip[j]));
    return {d, py};
}
ll closestDivideConquer(vector<P> p) {
    sort(p.begin(), p.end());
    return closestRec(p).first;
}

Complexity: $O(n\log n)$: the recursion does $O(n)$ work per level (the merge and the strip scan) over $O(\log n)$ levels.

Proof: the strip bound, and the algorithm's correctness

The $O(1)$-per-point bound. Fix a point $p$ in the strip and consider the window of other strip points $q$ with $p.y\le q.y<p.y+d$. Split the $L$-side points in this window (all within $d$ of the split line in $x$, and within $d$ in $y$ by the window) by a grid of $\frac{d}{2}\times \frac{d}{2}$ cells, $2\times 2=4$ cells covering a $d\times d$ box. Any two points in the same cell are at most $\frac{d}{2}\sqrt{2}<d$ apart; since any two points of $L$ are known to be $\ge d_L\ge d$ apart (as $d=\min (d_L,d_R)$), no cell can hold two $L$-points — so at most 4 $L$-points lie in the window. By the identical argument, at most 4 $R$-points lie in the window. So at most $8$ other strip points can be within $y$-distance $d$ of $p$, and since the inner loop scans in increasing $y$ and stops once the $y$-gap reaches $d$, it examines $O(1)$ (at most 8) candidates per point — $O(n)$ for the whole strip. (This is the classical “seven/eight-point” packing bound; the exact best constant doesn’t matter for the asymptotics, only that it is $O(1)$.)
Why the inner loop’s early-exit bound (7, not 8) works via sorted $y$: scanning points in increasing $y$-order from a fixed strip point $i$, the loop can stop as soon as $(\text{strip}[j].y-\text{strip}[i].y{)}^2\ge d$, since every later point is even farther in $y$ and the rectangle argument above bounds how many points can be within $d$ in $y$ before exceeding it — giving $O(1)$ amortized comparisons per point, $O(n)$ for the whole strip.
Correctness of the divide step. Any pair closer than $d=\min (d_L,d_R)$ cannot be entirely within $L$ or entirely within $R$ (those are already covered by $d_L,d_R$ respectively), so it must straddle the split — meaning both points are within $d$ of the split line (if either were farther, the pair’s distance would already exceed $d$), i.e. both are in the collected strip. The strip scan, sorted by $y$ and bounded to $O(1)$ comparisons per point by the lemma, finds this pair if it exists. $■$

Tested against brute-force $O(n^2)$ and against closestPair (Method 1) on random point sets including many duplicates, collinear configurations, and inputs designed to stress the strip (many points clustered near the median $x$).

The faster practical alternative

The randomised incremental method above is about the same code length, has no recursion, and is genuinely expected $O(n)$ — the usual choice when $n$ is large and recursion overhead matters.

The Bentley-Shamos pattern

The same skeleton — split by a median line, recurse, merge across a narrow strip — solves many problems:

Problem	Merge step
Closest pair	strip of width $2d$, 7-point lemma
All nearest neighbours	strip plus per-point candidate lists
Maximal points / Pareto frontier	merge the two staircases
Dominance counting	count left points dominating right points during merge
Smallest enclosing circle (Shamos)	superseded by Welzl
CDQ divide and conquer	left half’s updates applied to right half’s queries

That last row is the important one for competitive programming: CDQ divide and conquer is exactly this pattern applied to offline queries with several dimensions, and it is the standard tool for 3D partial-order counting problems — see 02-cdq-divide-and-conquer.

Variants / Use Cases

• KD-Tree — the structure for repeated nearest-neighbour queries.
• Sweep line — the other general geometric technique.