cp-algorithms · izanbf1803 · Jun 24, 2025 · Jun 24, 2025 · Jun 26, 2025 · Jun 26, 2025
diff --git a/src/geometry/nearest_points.md b/src/geometry/nearest_points.md
@@ -162,6 +162,166 @@ mindist = 1E20;
 rec(0, n);
 ```
 
+## Linear time randomized algorithms
+
+### A randomized algorithm with linear expected time
+
+An alternative method arises from a very simple idea to heuristically improve the runtime: We can divide the plane into a grid of $d \times d$ squares, then it is only required to test distances between same-block or adjacent-block points (unless all squares are disconnected from each other, but we will avoid this by design), since any other pair has a larger distance than the two points in the same square.
+
+<div style="text-align: center;">
+    <img src="nearest_points_blocks_example.png" alt="Example of the squares strategy" width="350px">
+</div>
+
+
+We will consider only the squares containing at least one point. Denote by $n_1, n_2, \dots, n_k$ the number of points in each of the $k$ remaining squares. Assuming at least two points are in the same or in adjacent squares,  the time complexity is $\Theta(\sum_{i=1}^{k} n_i^2)$.
+
+??? info "Proof"
+	For the $i$-th square containing $n_i$ points, the number of pairs inside is $\Theta(n_i^2)$. If the $i$-th square is adjacent to the $j$-th square, then we also perform $n_i n_j \le \max(n_i, n_j)^2 \le n_i^2 + n_j^2$ distance comparisons. Notice that each cube has at most $8$ adjacent cubes, so we can bound the sum of all comparisons by $\Theta(\sum_{i=1}^{k} n_i^2)$. $\quad \blacksquare$
+
+Now we need to decide on how to set $d$ so that it minimizes $\Theta(\sum_{i=1}^{k} n_i^2)$.
+
+####  Choosing d
+
+We need $d$ to be an approximation of the minimum distance $d$, and the trick is to just sample $n$ distances randomly and choose $d$ to be the smallest of these distances. We now prove that the expected running time is linear.
+
+??? info "Proof"
+	Imagine the disposition of points in squares with a particular choice of $d$, say $x$. Consider $d$ a random variable, resulting from our sampling of distances. Let's define $C(x) := \sum_{i=1}^{k(x)} n_i(x)^2$ as the cost estimation for a particular disposition when we choose $d=x$. Now, let's define $\lambda(x)$ such that $C(x) = \lambda(x) \, n$. What is the probability that such choice $x$ survives the sampling of $n$ independent distances? If a single pair among the sampled ones has distance smaller than $x$, this arrangement will be replaced by the smaller $d$. Inside a square, at least a quarter of the pairs would raise a smaller distance (imagine four subsquares in every square, and use the pigeonhole principle), so we have $\sum_{i=1}^{k} \frac{1}{4} {n_i \choose 2}$ pairs which yield a smaller final $d$. This is, approximately, $\frac{1}{8} \sum_{i=1}^{k} n_i^2 = \frac{1}{8} \lambda(x) n$. On the other hand, there are about $\frac{1}{2} n^2$ pairs that can be sampled. We have that the probability of sampling a pair with distance smaller than $x$ is at least (approximately) 
+
+	$$\frac{\lambda(x) n / 8}{n^2 / 2} = \frac{\lambda(x)/4}{n}$$
+
+	so the probability of at least one such pair being chosen during the $n$ rounds (and therefore finding a smaller $d$) is 
+
+	$$1 - \left(1 - \frac{\lambda(x)/4}{n}\right)^n \ge 1 - e^{-\lambda(x)/4}$$
+
+	(we have used that $(1 + x)^n \le e^{xn}$ for any real number $x$, check [Bernoulli inequalities](https://en.wikipedia.org/wiki/Bernoulli%27s_inequality#Related_inequalities)). <br> Notice this goes to $1$ exponentially as $\lambda(x)$ increases. This hints that $\lambda$ will be small usually.
+
+
+	We have shown that $\Pr(d \le x) \ge 1 - e^{-\lambda(x)/4}$, or equivalently, $\Pr(d \ge x) \le e^{-\lambda(x)/4}$. We need to know $\Pr(\lambda(d) \ge \text{something})$ to be able to estimate its expected value. We notice that $\lambda(d) \ge \lambda(x) \iff d \ge x$. This is because making the squares smaller only reduces the number of points in each square (splits the points into other squares), and this keeps reducing the sum of squares. Therefore,
+
+	$$\Pr(\lambda(d) \ge \lambda(x)) = \Pr(d \ge x) \le e^{-\lambda(x)/4} \implies \Pr(\lambda(d) \ge t) \le e^{-t/4} \implies \mathbb{E}[\lambda(d)] \le \int_{0}^{+\infty} e^{-t/4} \, \mathrm{d}t = 4$$
+
+	(we have used that $E[X] = \int_0^{+\infty} \Pr(X \ge x) \, \mathrm{d}x$, check [Stackexchange proof](https://math.stackexchange.com/a/1690829)).
+
+	Finally, $\mathbb{E}[C(d)] = \mathbb{E}[\lambda(d) \, n] \le 4n$, and the expected running time is $O(n)$, with a reasonable constant factor. $\quad \blacksquare$
+
+#### Implementation of the algorithm
+
+The advantage of this algorithm is that it is straightforward to implement, but still has good performance in practise. We first sample $n$ distances and set $d$ as the minimum of the distances. Then we insert points into the "blocks" by using a hash table from 2D coordinates to a vector of points. Finally, just compute distances between same-block pairs and adjacent-block pairs. Hash table operations have $O(1)$ expected time cost, and therefore our algorithm retains the $O(n)$ expected time cost with an increased constant.
+
+```{.cpp file=nearest_pair_randomized}
+using ll = long long;
+using ld = long double;
+
+struct RealPoint {
+    ld x, y;
+    RealPoint() {}
+    RealPoint(ld x_, ld y_) : x(x_), y(y_) {}
+};
+using pt = RealPoint;
+
+struct CustomHash {
+    size_t operator()(const pair<ll,ll>& p) const {
+        static const uint64_t C = chrono::steady_clock::now().time_since_epoch().count();
+        return C ^ ((p.first << 32) ^ p.second);
+    }
+};
+
+ld dist(const pt& a, const pt& b) {
+    ld dx = a.x - b.x;
+    ld dy = a.y - b.y;
+    return sqrt(dx*dx + dy*dy);
+}
+
+pair<pt,pt> closest_pair_of_points_rand_reals(vector<pt> P) {
+    const ld eps = 1e-9;
+
+    int n = int(P.size());
+    assert(n >= 2);
+    unordered_map<pair<ll,ll>,vector<pt>,CustomHash> grid;
+    grid.reserve(n);
+
+    mt19937 prng(chrono::system_clock::now().time_since_epoch().count());
+    uniform_int_distribution<int> uniform(0, n-1);
+
+    ld d = dist(P[0], P[1]);
+    pair<pt,pt> closest = {P[0], P[1]};
+
+    auto consider_pair = [&](const pt& a, const pt& b) -> void {
+        ld ab = dist(a, b);
+        if (ab + eps < d) {
+            d = ab;
+            closest = {a, b};
+        }
+    };
+
+    for (int i = 0; i < n; ++i) {
+        int j = uniform(prng);
+        int k = uniform(prng);
+        while (j == k)
+            k = uniform(prng);
+        consider_pair(P[j], P[k]);
+    }
+
+    for (const pt& p : P)
+        grid[{ll(p.x/d), ll(p.y/d)}].push_back(p);
+
+    for (const auto& it : grid) { // same block
+        int k = int(it.second.size());
+        for (int i = 0; i < k; ++i) {
+            for (int j = i+1; j < k; ++j)
+                consider_pair(it.second[i], it.second[j]);
+        }
+    }
+
+    for (const auto& it : grid) { // adjacent blocks
+        auto coord = it.first;
+        for (int dx = 0; dx <= 1; ++dx) {
+            for (int dy = -1; dy <= 1; ++dy) {
+                if (dx == 0 and dy == 0) continue;
+                pair<ll,ll> neighbour = {
+                    coord.first  + dx, 
+                    coord.second + dy
+                };
+                for (const pt& p : it.second) {
+                    if (not grid.count(neighbour)) continue;
+                    for (const pt& q : grid.at(neighbour))
+                        consider_pair(p, q);
+                }
+            }
+        }
+    }
+
+    return closest;
+}
+```
+
+
+### An alternative randomized linear expected time algorithm
+
+Now we introduce a different randomized algorithm which is less practical but very easy to show that it runs in expected linear time.
+
+- Permute the $n$ points randomly
+- Take $\delta := \operatorname{dist}(p_1, p_2)$
+- Partition the plane in squares of side $\delta/2$
+- For $i = 1,2,\dots,n$:
+	- Take the square corresponding to $p_i$
+	- Interate over the $25$ squares within two steps to our square in the grid of squares partitioning the plane
+	- If some $p_j$ in those squares has $\operatorname{dist}(p_j, p_i) < \delta$, then
+		- Recompute the partition and squares with $\delta := \operatorname{dist}(p_j, p_i)$
+		- Store points $p_1, \dots, p_i$ in the corresponding squares
+	- else, store $p_i$ in the corresponding square
+- output $\delta$
+
+While this algorithm may look slow, because of recomputing everything multiple times, we can show the total expected cost is linear. 
+
+??? info "Proof"
+	Let $X_i$ the random variable that is $1$ when point $p_i$ causes a change of $\delta$ and a recomputation of the data structures, and $0$ if not. It is easy to show that the cost is $O(n + \sum_{i=1}^{n} i X_i)$, since on the $i$-th step we are considering only the first $i$ points. However, turns out that $\Pr(X_i = 1) \le \frac{2}{i}$. This is because on the $i$-th step, $\delta$ is the distance of the closest pair in $\{p_1,\dots,p_i\}$, and $\Pr(X_i = 1)$ is the probability of $p_i$ belonging to the closest pair, which only happens in $2(i-1)$ pairs out of the $i(i-1)$ possible pairs (assuming all distances are different), so the probability is at most $\frac{2(i-1)}{i(i-1)} = \frac{2}{i}$, since we previously shuffled the points uniformly.
+
+	We can therefore see that the expected cost is 
+
+	$$O\left(n + \sum_{i=1}^{n} i \Pr(X_i = 1)\right) \le O\left(n + \sum_{i=1}^{n} i \frac{2}{i}\right) = O(3n) = O(n) \quad \quad \blacksquare$$ 
+
+
 ## Generalization: finding a triangle with minimal perimeter
 
 The algorithm described above is interestingly generalized to this problem: among a given set of points, choose three different points so that the sum of pairwise distances between them is the smallest.

diff --git a/src/geometry/nearest_points_blocks_example.png b/src/geometry/nearest_points_blocks_example.png