C++

Closest Pair and Nearest-Neighbor Queries

Find the two closest points among a million in O(n log n) by divide and conquer or a sweep, or in expected O(n) with a randomized grid, and answer nearest-neighbor queries with a grid.

Two points too close together

A mesh generator must not place two vertices closer than some minimum spacing, or it produces slivers. A collision checker for a crowd of particles needs to know whether any two are touching. A data-cleaning step wants to know whether a point cloud contains near-duplicates. All three ask the same question: among n points, which two are closest to each other?

The obvious program compares every pair. That is n * (n - 1) / 2 distance computations: half a million for a thousand points, half a trillion for a million. The first number is nothing; the second is hours. This lesson builds three algorithms that answer the question for a million points in a few million steps, proves why the clever one is correct, and tests all of them against the obvious program on inputs full of ties.

The closest pair is a cousin of the nearest-neighbor query from the spatial indexes lesson. You could build a k-d tree and query it once per point, which costs O(n log n) on friendly data. The algorithms here are simpler, have guarantees that do not depend on the data being friendly, and are the standard examples of three design techniques: divide and conquer, plane sweep and randomized incremental construction.

As in the rest of the course, coordinates are integers with |x|, |y| <= 10^9 and we compare squared distances in long long, which are exact. Two different pairs can be equally close, so every test in this lesson compares the closest distance, never the pair: any of the tied pairs is a correct answer.

The baseline, by hand and in code

Here are eight points that will serve as the worked example for the whole lesson:

(1,9) (2,3) (5,6) (6,1) (7,7) (8,4) (11,8) (12,2)

The brute-force program looks at all 28 pairs and keeps the smallest squared distance. It is the reference that every faster algorithm will be checked against, so it should stay this simple.

main.cppc++23
#include <limits>
#include <print>
#include <vector>

struct Point { long long x, y; };

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

struct Pair { long long d2; Point a, b; };

// Every pair once: n * (n - 1) / 2 squared distances.
Pair bruteForce(const std::vector<Point>& p, long& count) {
    Pair best{std::numeric_limits<long long>::max(), {}, {}};
    for (size_t i = 0; i < p.size(); ++i)
        for (size_t j = i + 1; j < p.size(); ++j) {
            ++count;
            if (long long d = dist2(p[i], p[j]); d < best.d2) best = {d, p[i], p[j]};
        }
    return best;
}

int main() {
    std::vector<Point> pts{{1, 9}, {2, 3}, {5, 6}, {6, 1},
                           {7, 7}, {8, 4}, {11, 8}, {12, 2}};
    long count = 0;
    Pair best = bruteForce(pts, count);
    std::println("closest pair ({}, {}) and ({}, {}), squared distance {}",
                 best.a.x, best.a.y, best.b.x, best.b.y, best.d2);
    std::println("distances computed: {}", count);

    // Four pairs share the smallest distance: any of them is a correct answer.
    std::vector<Point> square{{0, 0}, {3, 0}, {3, 3}, {0, 3}};
    count = 0;
    best = bruteForce(square, count);
    std::println("square: ({}, {}) and ({}, {}), squared distance {}",
                 best.a.x, best.a.y, best.b.x, best.b.y, best.d2);
}
Output

Press Run to compile and run this example.

closest pair (5, 6) and (7, 7), squared distance 5
distances computed: 28
square: (0, 0) and (3, 0), squared distance 9

The closest pair is (5,6) and (7,7), at squared distance 2*2 + 1*1 = 5. The second call shows why a test must compare distances: in a square, four pairs share the smallest distance, and which one a program reports depends only on its loop order.

Divide and conquer: halves and a strip

Sort the points by x and cut them at the median into a left half and a right half. Solve each half recursively. The closest pair is now either inside the left half, inside the right half, or a cross pair with one point on each side. The first two come back from the recursion; call the smaller of their distances δ (delta). The whole difficulty is the cross pairs, because there are about n^2 / 4 of them.

In the example, the median is (7,7), so the dividing line is x = 7. The left half (1,9) (2,3) (5,6) (6,1) has its closest pair at squared distance 18, (2,3) and (5,6). The right half (7,7) (8,4) (11,8) (12,2) has squared distance 10, (7,7) and (8,4). So δ² = 10, δ ≈ 3.16.

Two observations cut the cross pairs down to a handful.

  1. A cross pair closer than δ has both points within δ of the dividing line, because the horizontal gap alone would otherwise be at least δ. So only the points in the vertical strip |x - 7| < δ matter: here (5,6) (6,1) (7,7) (8,4).
  2. Sort the strip by y. A point only needs to be compared with the points after it in that order whose y differs by less than δ, and there are never more than seven of them.

Scanning the strip in y order (6,1) (8,4) (5,6) (7,7) finds (5,6) and (7,7) at squared distance 5, the cross pair that beats both halves. The drawing shows the example: the dividing line, the strip of half-width δ shaded, the two halves' best pairs in blue, and the winning cross pair in red. As in every drawing in this course, the program flips the y-axis because SVG's y grows downward.

main.cppc++23
#include <algorithm>
#include <cmath>
#include <format>
#include <limits>
#include <print>
#include <string>
#include <vector>

struct Point { long long x, y; };

class Svg {
public:
    Svg(double width, double height) : w_(width), h_(height) {}
    void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
        add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
            a.x, fy(a.y), b.x, fy(b.y), color, width);
    }
    void dot(Point p, const char* color = "#1f2937", double r = 3) {
        add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
    }
    void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
                 const char* fill = "none") {
        std::string list;
        for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
        add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
    }
    void text(Point p, const std::string& s, double size = 8) {
        add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
    }
    void print() const {
        std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
        std::print("{}", body_);
        std::println("</svg>");
    }

private:
    double fy(double y) const { return h_ - y; }  // SVG's y axis points down
    template <class... Args>
    void add(std::format_string<Args...> fmt, Args&&... args) {
        body_ += std::format(fmt, std::forward<Args>(args)...);
        body_ += '\n';
    }
    double w_, h_;
    std::string body_;
};

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

struct Pair { long long d2; Point a, b; };

Pair bruteForce(const std::vector<Point>& p) {
    Pair best{std::numeric_limits<long long>::max(), {}, {}};
    for (size_t i = 0; i < p.size(); ++i)
        for (size_t j = i + 1; j < p.size(); ++j)
            if (long long d = dist2(p[i], p[j]); d < best.d2) best = {d, p[i], p[j]};
    return best;
}

Point screen(Point p) { return {p.x * 20 + 10, p.y * 20 + 10}; }  // 20 px per unit

int main() {
    std::vector<Point> pts{{1, 9}, {2, 3}, {5, 6}, {6, 1},
                           {7, 7}, {8, 4}, {11, 8}, {12, 2}};
    std::ranges::sort(pts, [](Point a, Point b) {
        return a.x != b.x ? a.x < b.x : a.y < b.y;
    });
    size_t mid = pts.size() / 2;
    long long midX = pts[mid].x;
    Pair left = bruteForce({pts.begin(), pts.begin() + mid});
    Pair right = bruteForce({pts.begin() + mid, pts.end()});
    long long delta2 = std::min(left.d2, right.d2);
    double delta = std::sqrt(double(delta2));  // only for drawing the strip

    Svg svg(270, 210);
    long long x0 = std::llround((midX - delta) * 20 + 10);
    long long x1 = std::llround((midX + delta) * 20 + 10);
    svg.polygon({{x0, 0}, {x1, 0}, {x1, 210}, {x0, 210}}, "none", "rgba(37,99,235,0.15)");
    svg.line(screen({midX, 0}), screen({midX, 10}), "#1f2937", 1);
    svg.line(screen(left.a), screen(left.b), "#2563eb", 1);
    svg.line(screen(right.a), screen(right.b), "#2563eb", 1);

    // Points in the strip, sorted by y, each compared with the next few.
    std::vector<Point> strip;
    for (Point p : pts)
        if ((p.x - midX) * (p.x - midX) < delta2) strip.push_back(p);
    Pair best{delta2, left.d2 < right.d2 ? left.a : right.a,
              left.d2 < right.d2 ? left.b : right.b};
    std::ranges::sort(strip, [](Point a, Point b) {
        return a.y != b.y ? a.y < b.y : a.x < b.x;
    });
    for (size_t i = 0; i < strip.size(); ++i)
        for (size_t j = i + 1; j < strip.size(); ++j) {
            long long dy = strip[j].y - strip[i].y;
            if (dy * dy >= best.d2) break;
            if (long long d = dist2(strip[i], strip[j]); d < best.d2)
                best = {d, strip[i], strip[j]};
        }
    svg.line(screen(best.a), screen(best.b), "#dc2626", 2.5);
    for (Point p : pts) svg.dot(screen(p), "#1f2937", 4);
    svg.text({screen({midX, 0}).x + 3, 196}, std::format("x = {}", midX), 10);
    svg.print();
}
Output

Press Run to compile and run this example.

Why seven is enough

The second observation is the heart of the algorithm, so here is the argument. Take a point p in the strip and the rectangle that starts at p's height, extends δ upward and δ to each side of the dividing line: 2δ wide, δ tall. Any point that p must still check lies in that rectangle.

Split the rectangle into eight squares of side δ/2, four on each side of the line. Two points inside one such square are at most its diagonal apart, δ/√2, which is less than δ. But two points from the same half are at least δ apart, because δ is the best distance within each half. So each square holds at most one point of its own half. (A point exactly on the dividing line belongs to one half, and the squares of that half's side include the line.) That leaves at most eight points in the rectangle, p included, so at most seven others. The loop that walks up the strip and stops at the first point with dy >= δ therefore does constant work per point.

Two details keep the whole algorithm at O(n log n).

  • The recursion must return its points sorted by y, merging the two halves on the way back up exactly as merge sort does. Sorting the strip from scratch at every level costs O(n log n) per level and O(n log^2 n) in total.
  • The strip loop must stop at the first point that is too high, not merely skip it.

The program counts the distance computations, the elements passed through merges (the merge sort's O(n log n) work), and the largest number of strip successors any point examined. Points are uniform in a square of side 10^8.

main.cppc++23
#include <algorithm>
#include <limits>
#include <print>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

constexpr long long kMaxCoord = 1'000'000'000;
static_assert(8 * kMaxCoord * kMaxCoord < std::numeric_limits<long long>::max(),
              "squared distances between points within kMaxCoord must fit");
constexpr long long kNone = std::numeric_limits<long long>::max();  // fewer than 2 points

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
bool byX(Point a, Point b) { return a.x != b.x ? a.x < b.x : a.y < b.y; }
bool byY(Point a, Point b) { return a.y != b.y ? a.y < b.y : a.x < b.x; }

struct Stats { long distances = 0, merged = 0; int maxAhead = 0; };

// p[lo, hi) arrives sorted by x and leaves sorted by y (merge sort on the way
// back up). Returns the smallest squared distance among those points.
long long closest(std::vector<Point>& p, size_t lo, size_t hi, std::vector<Point>& strip,
                  Stats& st) {
    if (hi - lo <= 3) {  // base case: brute force, then sort by y
        long long best = kNone;
        for (size_t i = lo; i < hi; ++i)
            for (size_t j = i + 1; j < hi; ++j, ++st.distances)
                best = std::min(best, dist2(p[i], p[j]));
        std::sort(p.begin() + lo, p.begin() + hi, byY);
        return best;
    }
    size_t mid = lo + (hi - lo) / 2;
    long long midX = p[mid].x;  // left half: x <= midX, right half: x >= midX
    long long best = std::min(closest(p, lo, mid, strip, st),
                              closest(p, mid, hi, strip, st));
    std::inplace_merge(p.begin() + lo, p.begin() + mid, p.begin() + hi, byY);
    st.merged += long(hi - lo);

    // Only points closer than sqrt(best) to the dividing line can do better.
    strip.clear();
    for (size_t i = lo; i < hi; ++i)
        if ((p[i].x - midX) * (p[i].x - midX) < best) strip.push_back(p[i]);
    for (size_t i = 0; i < strip.size(); ++i) {
        int ahead = 0;
        for (size_t j = i + 1; j < strip.size(); ++j) {
            long long dy = strip[j].y - strip[i].y;
            if (dy * dy >= best) break;  // sorted by y: everything further is too far
            ++ahead, ++st.distances;
            best = std::min(best, dist2(strip[i], strip[j]));
        }
        st.maxAhead = std::max(st.maxAhead, ahead);
    }
    return best;
}

long long closestPair(std::vector<Point> p, Stats& st) {
    std::ranges::sort(p, byX);
    std::vector<Point> strip;
    return closest(p, 0, p.size(), strip, st);
}

int main() {
    Lcg rng{1};
    std::println("{:>8} {:>10} {:>14} {:>10} {:>10} {:>10}", "n", "distances",
                 "brute force", "dist / n", "merged / n", "max ahead");
    for (int n : {1'000, 10'000, 100'000, 1'000'000}) {
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, 100'000'000), rng.range(0, 100'000'000)});
        Stats st;
        closestPair(pts, st);
        std::println("{:>8} {:>10} {:>14} {:>10.2f} {:>10.2f} {:>10}", n, st.distances,
                     1LL * n * (n - 1) / 2, double(st.distances) / n,
                     double(st.merged) / n, st.maxAhead);
    }
}
Output

Press Run to compile and run this example.

       n  distances    brute force   dist / n merged / n  max ahead
    1000       1101         499500       1.10       8.93          2
   10000      13036       49995000       1.30      12.00          2
  100000     142627     4999950000       1.43      15.07          2
 1000000    1216922   499999500000       1.22      18.85          2

Distance computations stay between 1.1 and 1.5 per point at every size; the brute-force column is the number of pairs the obvious program would compute. Merged elements per point grow by about 3 to 4 for every tenfold increase in n, close to log2(10) ≈ 3.3: the log n factor lives in the sorting, not in the geometry. On random data no point ever looks at more than two strip successors. Seven is a guarantee for the worst case, not a typical value.

Testing with ties on purpose

The program above never sees two points with the same x, the same y or the same distance. Real data does: grid-snapped CAD vertices, pixel coordinates, duplicated records. The property test below runs the divide-and-conquer algorithm on 6,000 inputs built to be awkward: full lattices where dozens of pairs tie for the minimum, coordinates between 0 and 2 where duplicates are certain, all points on one vertical line (so every split has the same x), and plain random sets for contrast.

main.cppc++23
#include <algorithm>
#include <limits>
#include <print>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

constexpr long long kMaxCoord = 1'000'000'000;
static_assert(8 * kMaxCoord * kMaxCoord < std::numeric_limits<long long>::max(),
              "squared distances between points within kMaxCoord must fit");
constexpr long long kNone = std::numeric_limits<long long>::max();  // fewer than 2 points

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
bool byX(Point a, Point b) { return a.x != b.x ? a.x < b.x : a.y < b.y; }
bool byY(Point a, Point b) { return a.y != b.y ? a.y < b.y : a.x < b.x; }

struct Stats { long distances = 0, merged = 0; int maxAhead = 0; };

// p[lo, hi) arrives sorted by x and leaves sorted by y (merge sort on the way
// back up). Returns the smallest squared distance among those points.
long long closest(std::vector<Point>& p, size_t lo, size_t hi, std::vector<Point>& strip,
                  Stats& st) {
    if (hi - lo <= 3) {  // base case: brute force, then sort by y
        long long best = kNone;
        for (size_t i = lo; i < hi; ++i)
            for (size_t j = i + 1; j < hi; ++j, ++st.distances)
                best = std::min(best, dist2(p[i], p[j]));
        std::sort(p.begin() + lo, p.begin() + hi, byY);
        return best;
    }
    size_t mid = lo + (hi - lo) / 2;
    long long midX = p[mid].x;  // left half: x <= midX, right half: x >= midX
    long long best = std::min(closest(p, lo, mid, strip, st),
                              closest(p, mid, hi, strip, st));
    std::inplace_merge(p.begin() + lo, p.begin() + mid, p.begin() + hi, byY);
    st.merged += long(hi - lo);

    // Only points closer than sqrt(best) to the dividing line can do better.
    strip.clear();
    for (size_t i = lo; i < hi; ++i)
        if ((p[i].x - midX) * (p[i].x - midX) < best) strip.push_back(p[i]);
    for (size_t i = 0; i < strip.size(); ++i) {
        int ahead = 0;
        for (size_t j = i + 1; j < strip.size(); ++j) {
            long long dy = strip[j].y - strip[i].y;
            if (dy * dy >= best) break;  // sorted by y: everything further is too far
            ++ahead, ++st.distances;
            best = std::min(best, dist2(strip[i], strip[j]));
        }
        st.maxAhead = std::max(st.maxAhead, ahead);
    }
    return best;
}

long long closestPair(std::vector<Point> p, Stats& st) {
    std::ranges::sort(p, byX);
    std::vector<Point> strip;
    return closest(p, 0, p.size(), strip, st);
}

long long bruteForce(const std::vector<Point>& p) {
    long long best = kNone;
    for (size_t i = 0; i < p.size(); ++i)
        for (size_t j = i + 1; j < p.size(); ++j)
            best = std::min(best, dist2(p[i], p[j]));
    return best;
}

int main() {
    Lcg rng{77};
    long cases = 0, failed = 0;
    int maxAhead = 0;
    for (int t = 0; t < 6'000; ++t) {
        std::vector<Point> pts;
        int kind = t % 6;
        if (kind == 0) {  // a full lattice: many pairs share the minimum distance
            int k = rng.range(2, 9);
            for (int a = 0; a < k; ++a)
                for (int b = 0; b < k; ++b) pts.push_back({3 * a, 3 * b});
        } else {
            int n = rng.range(2, 80);
            int side = kind == 1 ? 2 : kind == 2 ? 8 : kind == 3 ? 40 : 1'000'000;
            for (int i = 0; i < n; ++i)
                pts.push_back({rng.range(0, side), rng.range(0, side)});
            if (kind == 4) for (Point& p : pts) p.x = 5;  // all on one vertical line
        }
        Stats st;
        ++cases;
        failed += closestPair(pts, st) != bruteForce(pts);
        maxAhead = std::max(maxAhead, st.maxAhead);
    }
    std::println("{} cases, {} failed", cases, failed);
    std::println("most strip successors examined for one point: {}", maxAhead);
}
Output

Press Run to compile and run this example.

6000 cases, 0 failed
most strip successors examined for one point: 3

All 6,000 agree with brute force, and even on lattices and vertical lines the strip loop looked at no more than three successors.

A sweep line with an ordered set

The plane sweep of Hinrichs, Nievergelt and Schorn (1988) gets the same bound with one pass and a balanced tree. Sort the points by x and move a vertical line from left to right, stopping at each point. Keep an active set, ordered by y, of the points less than δ behind the line, where δ is the best distance found so far.

At each new point cur:

  1. Remove from the active set every point at least δ to the left of cur. They can never be part of a better pair again, because every later point is even further right.
  2. Look at the active points whose y is within δ of cur.y. By the same packing argument there are only a constant number of them. std::set::lower_bound finds the first one in O(log n).
  3. Update δ if one of them is closer, then insert cur.

Each point is inserted once and erased at most once, and each step does one lower_bound: 3n set operations at O(log n) each. To keep the window in integers, the program uses reach = ceil(sqrt(best)), which is at least δ, to find the start of the window, and then compares exact squared distances.

main.cppc++23
#include <algorithm>
#include <cmath>
#include <limits>
#include <print>
#include <set>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

constexpr long long kNone = std::numeric_limits<long long>::max();
constexpr long long kMin = std::numeric_limits<long long>::min();

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
long long ceilSqrt(long long d2) {
    long long s = std::llround(std::sqrt(double(d2)));
    while (s * s < d2) ++s;
    while (s > 0 && (s - 1) * (s - 1) >= d2) --s;
    return s;
}

bool byX(Point a, Point b) { return a.x != b.x ? a.x < b.x : a.y < b.y; }
struct ByY {
    bool operator()(Point a, Point b) const { return a.y != b.y ? a.y < b.y : a.x < b.x; }
};

struct Stats { long distances = 0, setOps = 0; };  // each set operation is O(log n)

// Sweep a vertical line from left to right. The active set holds the points
// less than sqrt(best) behind the line, ordered by y.
long long closestPair(std::vector<Point> p, Stats& st) {
    std::ranges::sort(p, byX);
    std::set<Point, ByY> active;
    long long best = kNone, reach = 0;  // reach = ceilSqrt(best), once best is finite
    size_t tail = 0;                    // p[tail] is the oldest point still active
    for (Point cur : p) {
        if (best != kNone) {
            // Points at least sqrt(best) behind the sweep line can never win again.
            for (; (cur.x - p[tail].x) * (cur.x - p[tail].x) >= best; ++tail) {
                active.erase(p[tail]);
                ++st.setOps;
            }
        }
        // Start at the lowest point that can be within sqrt(best) of cur in y.
        auto first = best == kNone ? active.begin()
                                   : active.lower_bound({kMin, cur.y - reach});
        ++st.setOps;
        for (auto it = first; it != active.end(); ++it) {
            long long dy = it->y - cur.y;
            if (best != kNone && dy > 0 && dy * dy >= best) break;  // above the window
            ++st.distances;
            long long d = dist2(*it, cur);
            if (d < best) {
                best = d;
                if (best == 0) return 0;  // duplicates: nothing can beat zero
                reach = ceilSqrt(best);
            }
        }
        active.insert(cur);
        ++st.setOps;
    }
    return best;
}

int main() {
    Lcg rng{1};
    std::println("{:>8} {:>10} {:>10} {:>12}", "n", "distances", "set ops",
                 "set ops / n");
    for (int n : {1'000, 10'000, 100'000, 1'000'000}) {
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, 100'000'000), rng.range(0, 100'000'000)});
        Stats st;
        closestPair(pts, st);
        std::println("{:>8} {:>10} {:>10} {:>12.2f}", n, st.distances, st.setOps,
                     double(st.setOps) / n);
    }
    long failed = 0;
    for (int t = 0; t < 3'000; ++t) {
        int n = rng.range(2, 40), side = t % 2 ? 4 : 1'000;
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, side), rng.range(0, side)});
        if (t % 5 == 0) for (Point& q : pts) q.x = 7;  // all on a vertical line
        long long brute = kNone;
        for (int i = 0; i < n; ++i)
            for (int j = i + 1; j < n; ++j)
                brute = std::min(brute, dist2(pts[i], pts[j]));
        Stats st;
        failed += closestPair(pts, st) != brute;
    }
    std::println("3000 random small cases against brute force: {} failed", failed);
}
Output

Press Run to compile and run this example.

       n  distances    set ops  set ops / n
    1000         10       2998         3.00
   10000         14      29999         3.00
  100000         25     299999         3.00
 1000000         33    2999997         3.00
3000 random small cases against brute force: 0 failed

The set operations are exactly 3n minus a few, each O(log n). The distance computations are tiny for uniform data, because once δ is small the window around each point is almost always empty. The whole algorithm is one sort, one std::set and one loop, with no recursion.

Expected linear time with a randomized grid

The n log n above comes from sorting. A different idea avoids sorting altogether. Suppose you knew δ. Cut the plane into square cells of side δ, rounded up to a whole number, and put each point in its cell with a hash table. A point's partner closer than δ must be in the same cell or one of the eight around it. And each cell holds only a few points: with side exactly δ, at most four, because each of its four quarters is too small to hold two points that are δ apart (rounding the side up to an integer can add a few more when δ is tiny, still a constant).

You do not know δ in advance, so learn it incrementally:

  1. Shuffle the points. Let δ be the distance between the first two and build the grid for them.
  2. Insert the points one at a time. For each, look in its 3 by 3 block of cells.
  3. If it is closer to some point than δ, the cells are now too big: set the new δ and rebuild the grid from scratch with all points inserted so far.

The rebuild costs O(i) when it happens at the i-th point. It happens only if the i-th point is part of the closest pair among the first i points, and in a random order that has probability at most 2/i. The expected rebuild cost at step i is therefore at most 2/i * i = 2, and the expected total is O(n), as long as each hash-table operation takes expected constant time. This is the randomized incremental algorithm analyzed by Golin, Raman, Schwarz and Smid (1995); Rabin's 1976 algorithm and the sieve of Khuller and Matias (1995) reach the same expected bound by other routes.

The program runs the algorithm several times on the same points, shuffled differently each run, and averages the counts. Cells use floorDiv so that negative coordinates go to the correct cell, and the cell side is the smallest integer s with s * s >= best.

main.cppc++23
#include <algorithm>
#include <cmath>
#include <limits>
#include <print>
#include <unordered_map>
#include <utility>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() { s = s * 6364136223846793005ULL + 1442695040888963407ULL; return unsigned(s >> 33); }
    int range(int lo, int hi) { return lo + int(next() % unsigned(hi - lo + 1)); }  // inclusive
};

constexpr long long kNone = std::numeric_limits<long long>::max();

long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

// Smallest integer s with s * s >= d2: a cell side that is at least sqrt(d2).
long long ceilSqrt(long long d2) {
    long long s = std::llround(std::sqrt(double(d2)));
    while (s * s < d2) ++s;
    while (s > 0 && (s - 1) * (s - 1) >= d2) --s;
    return s;
}

long long floorDiv(long long a, long long b) { return a / b - (a % b < 0); }

struct Stats { long distances = 0, rebuilds = 0, reinserted = 0; };

// Randomized incremental closest pair with a uniform grid (expected O(n)).
long long closestPair(std::vector<Point> p, Lcg& rng, Stats& st) {
    if (p.size() < 2) return kNone;
    for (size_t i = p.size() - 1; i > 0; --i)  // Fisher-Yates shuffle
        std::swap(p[i], p[rng.next() % (i + 1)]);
    long long best = dist2(p[0], p[1]), side = 0;
    ++st.distances;
    std::unordered_map<unsigned long long, std::vector<Point>> grid;
    auto key = [&](long long cx, long long cy) {
        return (static_cast<unsigned long long>(cx) << 32) ^
               static_cast<unsigned long long>(cy & 0xffffffff);
    };
    auto rebuild = [&](size_t count) {  // grid of the first count points
        ++st.rebuilds;
        st.reinserted += long(count);
        side = ceilSqrt(best);
        grid.clear();
        for (size_t i = 0; i < count; ++i)
            grid[key(floorDiv(p[i].x, side), floorDiv(p[i].y, side))].push_back(p[i]);
    };
    if (best == 0) return 0;
    rebuild(2);
    for (size_t i = 2; i < p.size(); ++i) {
        long long cx = floorDiv(p[i].x, side), cy = floorDiv(p[i].y, side);
        long long found = best;
        // A point closer than sqrt(best) <= side lies in one of the 9 cells around.
        for (long long dx = -1; dx <= 1; ++dx)
            for (long long dy = -1; dy <= 1; ++dy) {
                auto it = grid.find(key(cx + dx, cy + dy));
                if (it == grid.end()) continue;
                for (Point q : it->second) {
                    ++st.distances;
                    found = std::min(found, dist2(p[i], q));
                }
            }
        if (found < best) {
            best = found;
            if (best == 0) return 0;
            rebuild(i + 1);  // the cells are too big now: start over with i + 1 points
        } else {
            grid[key(cx, cy)].push_back(p[i]);
        }
    }
    return best;
}

int main() {
    Lcg rng{1};
    std::println("{:>7} {:>16} {:>16} {:>14}", "n", "rebuilds per run",
                 "reinserted / n", "distances / n");
    for (int n : {1'000, 10'000, 100'000}) {
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(0, 100'000'000), rng.range(0, 100'000'000)});
        // The same points, shuffled differently each run; average over the runs.
        const int runs = std::max(10, 200'000 / n);
        Stats st;
        for (int r = 0; r < runs; ++r) closestPair(pts, rng, st);
        std::println("{:>7} {:>16.1f} {:>16.2f} {:>14.3f}", n, double(st.rebuilds) / runs,
                     double(st.reinserted) / runs / n, double(st.distances) / runs / n);
    }
    // Check against brute force on small inputs where ties are common.
    long failed = 0;
    for (int t = 0; t < 3'000; ++t) {
        int n = rng.range(2, 40), sideLen = t % 2 ? 4 : 1'000;
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i)
            pts.push_back({rng.range(-sideLen, sideLen), rng.range(-sideLen, sideLen)});
        long long brute = kNone;
        for (int i = 0; i < n; ++i)
            for (int j = i + 1; j < n; ++j)
                brute = std::min(brute, dist2(pts[i], pts[j]));
        Stats st;
        failed += closestPair(pts, rng, st) != brute;
    }
    std::println("3000 random small cases against brute force: {} failed", failed);
}
Output

Press Run to compile and run this example.

      n rebuilds per run   reinserted / n  distances / n
   1000             13.0             2.01          0.034
  10000             18.1             2.19          0.005
 100000             21.8             1.68          0.001
3000 random small cases against brute force: 0 failed

The number of rebuilds grows like 2 ln n (13, 18 and 22 here), and the points reinserted by all rebuilds together stay around 2n (between 1.7n and 2.2n here), as the analysis predicts. Distance computations are a tiny fraction of n: with cells of side δ, a new point's 3 by 3 block is usually empty. The work per point is a constant number of hash lookups.

There is no contradiction with the lower bound. Finding the closest pair needs Ω(n log n) operations in the algebraic decision tree model, the model in which sorting needs n log n comparisons. The grid algorithm steps outside that model twice: it rounds coordinates down to cell indices and it uses randomization and hashing. Its O(n) is an expected bound; an unlucky shuffle, or a hash table that degrades, makes one run slower.

Previous
Spatial Indexes: k-d Trees and Quadtrees
Next
Polygon Clipping: Sutherland-Hodgman and Boolean Operations