#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;

    answer.resize(2);
    answer[0] = a[L];
    answer[1] = a[L + 1];
    double d = distance(answer[0], answer[1]);
    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;
                answer[0] = a[i];
                answer[1] = a[j];
            }
        }
    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
// answer[0] = x_min, answer[1] = x_max
// answer[2] = y_min, answer[3] = y_max
void minXY_maxXY(std::vector<PAIR> &a, double * answer) {
    answer[0] = answer[2] = DBL_MAX;
    answer[1] = answer[3] = DBL_MIN;
    int left = 0;
    if (a.size() % 2 > 0) { // odd number of points
        answer[0] = answer[1] = a[0].first;
        answer[2] = answer[3] = a[0].second;
        ++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)
                answer[0] = a[i].first;
            if (answer[1] < a[i + 1].first - EPS)
                answer[1] = a[i + 1].first;
        }
        else {
            if (answer[0] > a[i + 1].first + EPS)
                answer[0] = a[i + 1].first;
            if (answer[1] < a[i].first - EPS)
                answer[1] = a[i].first;
        }
        // y-coordinates
        if (a[i].second < a[i + 1].second - EPS) {
            if (answer[2] > a[i].second + EPS)
                answer[2] = a[i].second;
            if (answer[3] < a[i + 1].second - EPS)
                answer[3] = a[i + 1].second;
        }
        else {
            if (answer[2] > a[i + 1].second + EPS)
                answer[2] = a[i + 1].second;
            if (answer[3] < a[i].second - EPS)
                answer[3] = a[i].second;
        }
    }
    answer[0] -= MARGIN;
    answer[1] += MARGIN;
    answer[2] -= MARGIN;
    answer[3] += MARGIN;
}

// 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.
    answer.resize(2);
    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;
    answer.resize(2);
    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;
                    answer[0] = pts[i];
                    answer[1] = pts[j];
                }
            }
            // 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;
                            answer[0] = pts[i];
                            answer[1] = q;
                        }
                    }
                }
        }
    }
    return d;
}