Closest pair of points

Why this matters

Many programs need the nearest neighbour of a point. A map app looks for the nearest shop. A game checks whether two objects touch. A data tool groups points that sit close together. The simplest way to find the closest pair is to measure the distance between every two points and keep the smallest. That works, but it gets slow fast. With 1,000 points there are 499,500 pairs to measure. With 100,000 points there are about 5 billion.

This lesson shows a divide and conquer method that finds the same answer in O(n log n) steps. O(n log n) (big O notation) means the work grows about as fast as n times log n as n grows. The idea is the same as merge sort: cut the problem in half, solve each half, and combine the answers. The new part is the combining step, which has to check pairs that cross the cut. The reader should already know recursion and the base case from the recursion module. The earlier lessons in this module, including the binary search ones, use the same halving pattern.

The problem in plain words

A point is a pair of numbers, written (x, y). Think of it as a spot on a sheet of graph paper. The x number tells how far right the spot is from the left edge. The y number tells how far up it is from the bottom edge.

The distance between two points is the length of the straight line between them. For points (x1, y1) and (x2, y2), the distance is the square root of the sum of two squared differences: the difference in x squared, plus the difference in y squared. "Squared" means multiplied by itself. "Square root" is the number that, multiplied by itself, gives the value inside.

Take the points (0, 0) and (3, 4). The x difference is 3, and 3 squared is 9. The y difference is 4, and 4 squared is 16. The sum is 25, and the square root of 25 is 5. So the distance is 5. You can check this with the Pythagorean rule for right triangles.

The task is this. You get n points. You must return the smallest distance between any two of them. You do not need to return which pair it is, although the code below could be changed to do so.

The simple method: check every pair

Before the faster method, look at the simple one, because the fast method uses it on small inputs.

Take the points one at a time. For each point, measure the distance to every later point in the list. Keep the smallest distance seen. Every pair is checked exactly once, so the answer is correct.

The cost is the number of pairs. With n points there are n times (n minus 1) divided by 2 pairs. For n = 4 that is 6 pairs. For n = 1,000 it is 499,500. Doubling n makes the work about four times as large. Work that grows with the square of the input size is called quadratic time, and it is the reason the simple method does not scale.

The idea: split, solve, then check the seam

The divide and conquer method follows four steps.

Sort the points by x, so that they run from left to right.
Cut the list in the middle. The left half holds the points with smaller x. The right half holds the points with larger x. Call the x value of the first point in the right half the cut line.
Solve each half by the same method, recursively. You get the smallest distance inside the left half, and the smallest distance inside the right half. Call the smaller of these two numbers d. So far, d is the best distance found.
Check for pairs that cross the cut. A pair with one point on each side might be closer than d. Only points that lie within d of the cut line can take part in such a pair, because a point farther away than d cannot be within d of anything on the other side. Those points form a strip, a narrow band around the cut line. Check the points in the strip against each other.

Step 4 is where the real work is. A strip can hold many points, so checking every pair in it could be slow. The method avoids that with one more ordering trick, shown in the next section.

Why only a few points in the strip need checking

Sort the points in the strip by y, so that they run from bottom to top. Now take one point P in the strip. The only points that can be within d of P are the points that are above P by less than d. Stop as soon as the y difference reaches d. Any point further up is too far from P, and so are all the points after it. Since the strip is sorted by y, you can stop early.

How many points can sit inside the box that is checked for each P? The box is d wide on each side of the cut, so 2d wide in total, and d tall. Split that box into 8 small squares, each d/2 wide and d/2 tall. Each small square lies entirely on one side of the cut. Two points on the same side are never closer than d, because d is no larger than the best distance inside either half. But two points inside one small square are closer than d, because the diagonal of the square is shorter than d. So each small square holds at most one point. Eight squares means at most 8 points. Those 8 include P, so each point needs to be compared with at most 7 points after it. That is a fixed number, no matter how many points there are.

This is why the strip step takes time proportional to the number of points in the strip, not to the square of that number.

Why the time is O(n log n)

Each call splits its list in half. Splitting a list of n items in half again and again takes about log2(n) levels before every piece holds one item. The log here means base 2 logarithm: the number of times you can halve n before reaching 1. For n = 1,024 the answer is 10, because 2 to the power 10 is 1,024.

At each level, the work adds up to about n steps. The cut merges the two sorted halves, which takes n steps. Building the strip takes n steps. Checking the strip takes n steps, since each point is compared with at most 7 others. So each level costs about n steps, and there are log n levels. The total is n times log n, written O(n log n). The first sort by x adds another n log n steps, which does not change the total.

Worked example 1: six points by hand

Here are six points. The letters are names for them.

A is (2, 3)
B is (3, 4)
C is (5, 1)
D is (12, 10)
E is (12, 30)
F is (40, 50)

Sorted by x, the order is A, B, C, D, E, F. D and E share x = 12. Either can come first after sorting, and the answer does not change. This trace puts D first.

The list has six points, so the cut is in the middle. The left half is A, B, C. The right half is D, E, F. The cut line is at x = 12, the x value of the first point of the right half.

Solve the left half. There are three points, so measure all three pairs.

A to B: the x difference is 1 and the y difference is 1. The distance is the square root of 1 plus 1, which is about 1.414.
A to C: the differences are 3 and 2. The distance is the square root of 9 plus 4, which is about 3.606.
B to C: the differences are 2 and 3. The distance is again the square root of 13, about 3.606.

The smallest is 1.414, for the pair A and B. So the left answer is d-left = 1.414.

Solve the right half the same way.

D to E: the y difference is 20 and the x difference is 0. The distance is 20.
D to F: the differences are 28 and 40. The distance is about 48.8.
E to F: the differences are 28 and 20. The distance is about 34.4.

The smallest is 20, so the right answer is d-right = 20.

The best so far is d = 1.414.

Now check the strip. Keep the points whose x is within 1.414 of the cut line at x = 12. That range runs from about 10.6 to 13.4. Only D and E are in it. Sort them by y: D has y = 10 and E has y = 30. The y gap is 20, which is larger than 1.414. The inner loop stops at once, so no pairs are measured.

The final answer is 1.414, from the pair A and B. The brute force method gives the same answer. The program below prints both. Its output is shown after the code.

The complete program

The program below puts the method in code. It keeps a record of how many distance checks it makes, so you can compare with the simple method later. It checks the fast method against brute force on 500 random sets of points, and it counts the checks on 100,000 points.

C++ · editable

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

struct Point {
    double x;
    double y;
};

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

bool by_y(const Point& a, const Point& b) {
    return a.y < b.y;
}

std::size_t checks = 0;

// Points p[lo..hi) are sorted by x when this starts.
// When it returns, p[lo..hi) is sorted by y.
double solve(std::vector<Point>& p, std::size_t lo, std::size_t hi) {
    std::size_t n = hi - lo;

if (n <= 3) {
        double best = std::numeric_limits<double>::infinity();
        for (std::size_t i = lo; i < hi; ++i) {
            for (std::size_t j = i + 1; j < hi; ++j) {
                ++checks;
                best = std::min(best, dist(p[i], p[j]));
            }
        }
        std::sort(p.begin() + lo, p.begin() + hi, by_y);
        return best;
    }

std::size_t mid = lo + n / 2;
    double mid_x = p[mid].x;

double d_left = solve(p, lo, mid);
    double d_right = solve(p, mid, hi);
    double d = std::min(d_left, d_right);

std::inplace_merge(p.begin() + lo, p.begin() + mid, p.begin() + hi, by_y);

std::vector<Point> strip;
    for (std::size_t i = lo; i < hi; ++i) {
        if (std::abs(p[i].x - mid_x) < d) {
            strip.push_back(p[i]);
        }
    }

for (std::size_t i = 0; i < strip.size(); ++i) {
        for (std::size_t j = i + 1; j < strip.size() && strip[j].y - strip[i].y < d; ++j) {
            ++checks;
            d = std::min(d, dist(strip[i], strip[j]));
        }
    }
    return d;
}

double closest_pair(std::vector<Point> pts) {
    if (pts.size() < 2) {
        return std::numeric_limits<double>::infinity();
    }
    std::sort(pts.begin(), pts.end(), [](const Point& a, const Point& b) { return a.x < b.x; });
    return solve(pts, 0, pts.size());
}

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

int main() {
    std::vector<Point> example = {
        {2, 3}, {12, 30}, {40, 50}, {5, 1}, {12, 10}, {3, 4}
    };
    std::cout << "example 1 fast:  " << closest_pair(example) << "\n";
    std::cout << "example 1 brute: " << brute_force(example) << "\n";

std::vector<Point> example2 = {
        {0, 0}, {1, 5}, {2, 1}, {3, 0}, {4, 6}, {5, 9}
    };
    std::cout << "example 2 fast:  " << closest_pair(example2) << "\n";
    std::cout << "example 2 brute: " << brute_force(example2) << "\n";

// Random check against brute force.
    unsigned seed = 12345;
    auto next = [&seed]() {
        seed = seed * 1103515245u + 12345u;
        return static_cast<double>((seed >> 8) % 100000) / 100.0;
    };
    int mismatches = 0;
    for (int trial = 0; trial < 500; ++trial) {
        int n = 2 + trial % 60;
        std::vector<Point> pts;
        for (int i = 0; i < n; ++i) {
            pts.push_back({next(), next()});
        }
        double a = closest_pair(pts);
        double b = brute_force(pts);
        if (std::abs(a - b) > 1e-9) {
            ++mismatches;
        }
    }
    std::cout << "random trials mismatches: " << mismatches << "\n";

// Operation count on a larger input.
    std::vector<Point> big;
    for (int i = 0; i < 100000; ++i) {
        big.push_back({next(), next()});
    }
    checks = 0;
    double best = closest_pair(big);
    std::cout << "n = 100000 closest distance: " << best << "\n";
    std::cout << "n = 100000 distance checks: " << checks << "\n";
    std::cout << "n * n / 2 for brute force: " << (100000.0 * 99999.0 / 2.0) << "\n";
}

Save the program as closest.cpp, then compile and run it with:

g++ -std=c++20 -Wall -o closest closest.cpp
./closest

The output on g++ 15.2.0 is:

example 1 fast:  1.41421
example 1 brute: 1.41421
example 2 fast:  1.41421
example 2 brute: 1.41421
random trials mismatches: 0
n = 100000 closest distance: 0.0141421
n = 100000 distance checks: 138765
n * n / 2 for brute force: 4.99995e+09

The distance check count is the number of times the program measures a distance inside the base case or the strip. For 100,000 points the brute force method needs 4,999,950,000 checks. The fast method needs 138,765. That gap grows with n, and it is the reason the method exists.

Reading the code

The Point struct holds two numbers. A struct is a named group of values. dist uses std::hypot, a standard function that returns the square root of the sum of two squares. It also avoids overflow for very large numbers, which a plain multiply could cause.

The function solve works on a range of the vector, from index lo up to but not including hi. This avoids copying the vector for each call. On entry the range is sorted by x. On exit it is sorted by y. That second fact is what lets the parent call merge the two halves cheaply.

Small ranges with three or fewer points use the brute force loop. Every pair is measured, and then the range is sorted by y, so the parent can use it.

For larger ranges, mid is the index of the first point in the right half. mid_x is the x value of that point, which is the cut line. Each half is solved by a recursive call. The two answers are combined in d.

std::inplace_merge joins the two sorted halves into one list sorted by y. It works in place, meaning it reorders the same vector. The library may use a temporary buffer internally. It needs each half to be sorted by y already, and the recursive calls guarantee that.

The strip is built by keeping points whose x is closer than d to the cut line. Because the range is now sorted by y, the strip is also sorted by y. The nested loop compares each strip point with the points after it. The loop stops when the y gap reaches d, as explained earlier.

The closest_pair function is the entry point. It sorts the points by x, then calls solve on the whole list. It returns infinity when there are fewer than two points, because a single point has no pair. The std::numeric_limits value is a standard way to write infinity for a double.

A second example: the answer crosses the cut

In the first example the answer lay entirely inside the left half. This example shows the strip doing real work, because the closest pair has one point on each side of the cut.

Take six points, already sorted by x.

P is (0, 0)
Q is (1, 5)
R is (2, 1)
S is (3, 0)
T is (4, 6)
U is (5, 9)

The left half is P, Q, R. The right half is S, T, U. The cut line is at x = 3, the x of S.

The left half has three points, so measure all pairs.

P to Q: the distance is the square root of 1 plus 25, about 5.10.
P to R: the distance is the square root of 4 plus 1, about 2.236.
Q to R: the distance is the square root of 1 plus 16, about 4.12.

The left answer is d-left = 2.236.

The right half is treated the same way.

S to T: the square root of 1 plus 36, about 6.08.
S to U: the square root of 4 plus 81, about 9.22.
T to U: the square root of 1 plus 9, about 3.16.

The right answer is d-right = 3.16. The best so far is d = 2.236.

Now the strip. Keep points whose x is within 2.236 of x = 3. That range runs from about 0.76 to 5.24. Point P at x = 0 is out. Q at x = 1, R at x = 2, S at x = 3, T at x = 4, and U at x = 5 are in.

Sort the strip by y. The y values are S at 0, R at 1, Q at 5, T at 6, and U at 9. So the order is S, R, Q, T, U.

Now run the inner loop.

Start with S. Compare with R. The y gap is 1, which is less than 2.236, so measure. The distance from S (3, 0) to R (2, 1) is the square root of 1 plus 1, about 1.414. That is smaller than d, so d becomes 1.414.
Compare S with Q next. The y gap is 5, which is not less than d (now 1.414). The loop stops for S.
Start with R. Compare with Q. The y gap is 4, which is more than 1.414. Stop.
Start with Q. Compare with T. The y gap is 1, which is less than 1.414, so measure. The distance is the square root of 9 plus 1, about 3.16. That is not smaller, so d stays 1.414.
Compare Q with U. The y gap is 4. Stop.
Start with T. Compare with U. The y gap is 3. Stop.

The answer is 1.414, from the pair R and S. R is on the left of the cut and S is on the right. Neither half alone contains this pair, so the strip check is the only way the method finds it. The program prints 1.41421 for this set.

A third example: why the strip must be sorted

The third example shows what goes wrong when the strip is not sorted by y. Take three strip points: A at (3, 0), B at (3, 10), and C at (2, 1). Suppose d is 2.0, and the strip is in the order A, B, C.

Start with A. The loop checks the next point, B. The y gap is 10 minus 0, which is 10. That is not less than d, so the loop stops at once. C is never compared with A.

The pair A and C is at distance about 1.414, which is less than d. So the correct answer is 1.414, but the loop returns 2.0. The program below runs this loop on the three points and prints the comparisons it makes.

C++ · editable

Its output is:

compare 1 with 2
result 2.00000

The loop compared B with C, not A with C. The stop test assumes that the points after the current one have larger y values, so the first point that is far enough up ends the loop. In a sorted strip that is always true. In this unsorted strip it is not, and the pair A and C is skipped. The sort by y in the base case and the merge step exist to make that assumption true.

Common mistakes

Forgetting the strip check

If the combining step returns only the smaller of the two half answers, the program misses pairs that cross the cut. The second example shows this. Remove the strip loop from the program above, and run it again. The output for example 2 becomes:

example 2 fast:  2.23607
example 2 brute: 1.41421

The fast answer is wrong because the pair R and S was never checked. The result of the wrong version is a number that looks valid, so the mistake is easy to miss. The random test against brute force catches it, which is why the program includes that test.

Leaving out the include for std::sort

The program uses std::sort and std::inplace_merge. These come from the <algorithm> header. If that include line is removed, g++ stops with an error. The first lines of the g++ 15.2 output are (line numbers depend on your file):

no_algorithm.cpp: In function ‘double solve(std::vector<Point>&, std::size_t, std::size_t)’:
no_algorithm.cpp:36:14: error: ‘sort’ is not a member of ‘std’; did you mean ‘sqrt’?
   36 |         std::sort(p.begin() + lo, p.begin() + hi, by_y);
      |              ^~~~
      |              sqrt

The message names the missing piece, sort, and suggests sqrt, which is a close name in <cmath>. The fix is to add #include <algorithm> at the top.

Sorting only by x inside the recursion

A common error is to sort by x at every level, or to forget to sort by y at all. The merge step needs both halves sorted by y. The strip step needs the strip sorted by y. If you sort by x at each level, the merge does not produce a y-sorted list, and the third example shows what can go wrong.

Measuring the strip without the stop test

Some people write the strip loop as two nested loops that check every pair. The result is still correct, but the time is no longer O(n log n) in the worst case, because the strip can hold many points. The stop test on the y gap is what limits each point to at most 7 comparisons.

Using the wrong value for the cut

The cut line must be the x value of a point, or any value between the two halves. Using the x of the last point of the left half also works, as long as the strip test uses the same value. Mixing two different values in the strip test is a common source of skipped pairs.

Practice

Exercise 1: trace four points

Take the points (0, 0), (5, 5), (2, 1), and (4, 0). Sort them by x. Split them into two halves of two points each. Find d-left and d-right, then build the strip around the cut and find the answer. Give the exact distance and the pair.

Show the answer

Sorted by x the points are (0, 0), (2, 1), (4, 0), (5, 5). The cut is at the x of the third point, which is 4. The left half is (0, 0) and (2, 1). Its distance is the square root of 4 plus 1, about 2.236. The right half is (4, 0) and (5, 5). Its distance is the square root of 1 plus 25, about 5.099. So d is 2.236.

The strip keeps points whose x is within 2.236 of 4. That range runs from about 1.76 to 6.24. Points (2, 1), (4, 0), and (5, 5) are in it. Sorted by y, the strip is (4, 0), (2, 1), (5, 5).

From (4, 0), the next point (2, 1) has a y gap of 1, so measure it. The distance is the square root of 4 plus 1, which is about 2.236. That equals d, so d stays the same. The next point (5, 5) has a y gap of 5, so stop.

From (2, 1), the next point (5, 5) has a y gap of 4, which is more than d, so stop.

The answer is about 2.236, which is the square root of 5. Two pairs reach this distance: (0, 0) with (2, 1), and (4, 0) with (2, 1). The program's answer for this input is 2.23607, and brute force gives the same.

Exercise 2: why 7 points

The text says each point in the strip needs at most 7 comparisons. Explain in your own words why 8 points is the most that can fit in one box.

Show the answer

The box is d wide on each side of the cut and d tall, so it is 2d wide and d tall. Cut it into 8 squares, each d/2 by d/2. The diagonal of a square with side d/2 is d divided by the square root of 2, which is about 0.71 d. Each square lies on one side of the cut. Two points on the same side are never closer than d, because d is no larger than the best distance inside either half. But any two points in one square are closer than that diagonal, so closer than d. So each square holds at most one point. Eight squares means at most 8 points. The point being checked is one of them, so it is compared with at most 7 others.

Exercise 3: count the checks

Use the program from this lesson as a starting point. Write a second program that makes random points with integer x and y values from 0 to 99,999, and counts the distance checks for 1,000 points and for 10,000 points. Compare each count with n times log2(n), where n is the number of points.

Show the answer

In one run with that kind of input, the counts were 1,167 for 1,000 points and 12,983 for 10,000 points. For n = 1,000, log2(1,000) is about 9.97, so n times log2(n) is about 9,966. The count is about 0.12 of that. For n = 10,000, log2(10,000) is about 13.29, so n times log2(n) is about 132,877. The count is about 0.10 of that. The ratio stays roughly the same as n grows, which matches the O(n log n) bound. The brute force count for the same sizes is 499,500 and 49,995,000, so the gap grows as n grows.

Exercise 4: why the base case sorts

In the base case, the code sorts a list of three or fewer points by y, after measuring the pairs. Explain why this sort is needed, even though the base case already has its answer.

Show the answer

The base case returns only a distance, but its parent needs the list in y order. The parent calls std::inplace_merge on the two halves, and that call assumes each half is sorted by y. If the base case did not sort, the merged list would not be sorted by y, the strip would not be sorted by y, and the early stop in the strip loop could skip a pair. So the sort in the base case is part of the work the parent depends on.

Exercise 5: a change to the program

Change closest_pair so that it also returns the two points that are closest. Describe the change you would make to solve, and say what it returns.

Show the answer

Let solve return a small struct that holds the distance and the two points. Each place that sets best or d also stores the two points that produced that distance. The base case stores the pair it measured. The strip loop stores the pair it measured. When the function combines the two halves, it returns the struct with the smaller distance. The closest_pair function returns that struct. The main cost of the change is that every update to d also copies two points, which is small compared with the distance measurement itself.

Summary

The simple method measures every pair, which takes time that grows with the square of n.
The divide and conquer method sorts by x, solves each half, and checks a strip around the cut line.
Points in the strip are sorted by y, and each point needs at most 7 comparisons with later points.
The total work is O(n log n), because each of the log n levels does about n steps.
Sorting by y in the base case and in the merge step is what makes the early stop in the strip loop safe.

Closest pair of points

Why this matters

The problem in plain words

The simple method: check every pair

The idea: split, solve, then check the seam

Why only a few points in the strip need checking

Why the time is O(n log n)

Worked example 1: six points by hand

The complete program

Reading the code

A second example: the answer crosses the cut

A third example: why the strip must be sorted

Common mistakes

Forgetting the strip check

Leaving out the include for std::sort

Sorting only by x inside the recursion

Measuring the strip without the stop test

Using the wrong value for the cut

Practice

Exercise 1: trace four points

Exercise 2: why 7 points

Exercise 3: count the checks

Exercise 4: why the base case sorts

Exercise 5: a change to the program

Summary

Stuck? Ask the tutor