# Closest points using Rabin randomizing approach

I was told to use the following Rabin algorithm to find the shortest distance between 2 points in 2D:

• Randomly choose sqrt(n) and brute force to find the closest distance in this set. Let this distance be delta.

• Make a delta grid on top of all n points, which partitions all n points into different bins/buckets.

• For each bin k, find the shortest distance between 2 points in this bin. Also, find the shortest distance between points in bin k and points in its 8 neighbor bins (4 should be enough to avoid double counting).

• The result is the smallest distance we compute.

The expected running time for this algorithm is said to be $O(n)$. However, when I tried to implement it and compared with the divide-and-conquer (DC) algorithm (supposed to be O(n log n)), although the results were the same, it was not as fast as expected.

In particular, testing on my i7-4700mq with MS Visual Studio Community 2015, when n <= 5 000 000, it was only about 0.1 second faster than the DC algorithm. Then when n was larger, the DC algorithm began to overtake it, and when n = 50 000 000, it took the DC algorithm about 70 seconds and this randomizing algorithm about 100 seconds.

Is there anything wrong in my code? I was told to use this algorithm with the "being told" expected running time. I found very little information about this algorithm online.

#include <algorithm>
#include <vector>
#include <climits>
#include <unordered_map>
#include <time.h>
#include <random>

#define PAIR std::pair<double, double> // x y coordinates
#define EPS 1e-8
#define MARGIN 1e-3

// Eucledian distance
double distance(PAIR p, PAIR q) {
return sqrt((p.first - q.first) * (p.first - q.first) + (p.second - q.second) * (p.second - q.second));
}

// brute force for small data set
double closest_2d_points_slow(std::vector<PAIR> &a, int L, int R, std::vector<PAIR> &answer) {
if (R - L < 1) return DBL_MAX;

for (int i = L; i < R; ++i)
for (int j = i + 1; j <= R; ++j) {
double dd = distance(a[i], a[j]);
if (dd < d - EPS) {
d = dd;
}
}
return d;
}

// estimate 'delta' by brute force on sqrt(n) points.
double estimate_delta(std::vector<PAIR> &a) {
int n = (int)a.size(), k = (int)sqrt(n), k1 = k;
std::vector<PAIR> sample(k);

// choose k = sqrt(n) points randomly.
std::mt19937 rng((unsigned int)time(0));
for (int i = 0; i < n; ++i) {
std::uniform_int_distribution<int> gen(0, n - i - 1);
if (gen(rng) < k1) {
sample[--k1] = a[i];
}
}
return closest_2d_points_slow(sample, 0, k - 1, std::vector<PAIR>());
}

// find x_min, x_max, y_min, y_max
void minXY_maxXY(std::vector<PAIR> &a, double * answer) {
int left = 0;
if (a.size() % 2 > 0) { // odd number of points
++left;
}
for (int i = left; i < (int)a.size(); i += 2) {
// x-coordinates
if (a[i].first < a[i + 1].first - EPS) {
if (answer[0] > a[i].first + EPS)
if (answer[1] < a[i + 1].first - EPS)
}
else {
if (answer[0] > a[i + 1].first + EPS)
if (answer[1] < a[i].first - EPS)
}
// y-coordinates
if (a[i].second < a[i + 1].second - EPS) {
if (answer[2] > a[i].second + EPS)
if (answer[3] < a[i + 1].second - EPS)
}
else {
if (answer[2] > a[i + 1].second + EPS)
if (answer[3] < a[i].second - EPS)
}
}
}

// convert xy-coordinate to a number representing a cell in a rows x cols grid.
int sub2ind(PAIR p, int rows, int cols, double x_min, double y_min, double delta) {
return int(floor((p.second - y_min) / delta)) * cols + int(floor((p.first - x_min) / delta));
}

// convert a cell number to (row, col)
void ind2sub(int bin, int rows, int cols, int &row, int &col) {
row = bin / cols;
col = bin % cols;
}

// return the neighbors of a cell
void bin_neighbors(int bin, int rows, int cols, std::vector<int> &neighbors) {
int row, col;
ind2sub(bin, rows, cols, row, col);
neighbors.clear();
if (row == 0) {
neighbors.push_back(bin + cols);
if (col > 0)
neighbors.push_back(bin + cols - 1);
if (col < cols - 1) {
neighbors.push_back(bin + 1);
neighbors.push_back(bin + cols + 1);
}
return;
}
if (row == rows - 1) {
if (col < cols - 1)
neighbors.push_back(bin + 1);
return;
}
neighbors.push_back(bin + cols);
if (col > 0)
neighbors.push_back(bin + cols - 1);
if (col < cols - 1) {
neighbors.push_back(bin + 1);
neighbors.push_back(bin + cols + 1);
}
}

// Main function to compute the cloest 2d distance between 2 points using Rabin randomized algorithm.
double closest_2d_points_rabin_randomized(std::vector<PAIR> &a, std::vector<PAIR> &answer) {
if (a.size() < 25)
return closest_2d_points_slow(a, 0, (int)a.size() - 1, answer);

// compute estimated delta and use it as a grid.
double delta = estimate_delta(a);
double mm[4];
minXY_maxXY(a, mm);

// compute the number of rows and cols in a grid.
int rows = (int)ceil((mm[3] - mm[2]) / delta);
int cols = (int)ceil((mm[1] - mm[0]) / delta);

// project each point into cells according to its coordinate and delta.
std::unordered_map<int, std::vector<PAIR>> point2bin;
point2bin.reserve(a.size());

std::vector<int> bins;
for (PAIR p : a) {
int bin = sub2ind(p, rows, cols, mm[0], mm[2], delta);
if (point2bin.find(bin) == point2bin.end()) {
bins.push_back(bin);
std::vector<PAIR> list(1, p);
point2bin.insert(std::pair<int, std::vector<PAIR>>(bin, list));
}
else
point2bin[bin].push_back(p);
}

// look at each cell and its 4 adjacent cells.
double d = DBL_MAX, dd = DBL_MAX;
for (int bin : bins) {
std::vector<PAIR> pts = point2bin[bin];
int bin_size = (int)pts.size();
std::vector<int> neighbors;
bin_neighbors(bin, rows, cols, neighbors);
for (int i = 0; i < bin_size; ++i) {
// points in one cell
for (int j = i + 1; j < bin_size; ++j) {
dd = distance(pts[i], pts[j]);
if (dd < d - EPS) {
d = dd;
}
}
// points in its 4 neighbors
for (int another_bin : neighbors)
if (point2bin.find(another_bin) != point2bin.end()) {
for (PAIR q : point2bin[another_bin]) {
dd = distance(pts[i], q);
if (dd < d - EPS) {
d = dd;