Sample code for closest pair of points

#include <climits>
#include <cmath>
#include <algorithm>
#include <chrono>
#include <iostream>
#include <numeric>
#include <random>
#include <utility>

template<typename CT>
struct Point {
    typedef CT CoordinateType;

    CoordinateType x;
    CoordinateType y;
};

template<typename CT>
std::ostream &operator<<(std::ostream &os, const Point<CT> &p)
{
    return os << "x:" << p.x << " y:" << p.y;
}

template<typename T1, typename T2>
std::ostream &operator<<(std::ostream &os, const std::pair<T1, T2> &pair)
{
    return os << " (" << pair.first << "," << pair.second << ") ";
}

template<typename CT>
CT squareRoot(CT n) {
    return static_cast<CT>(std::pow(n, 0.5));
}

template<typename CT>
CT distanceSquared(const Point<CT> &first, const Point<CT> &second) {
    const CT dx = (first.x - second.x), dy = (first.y - second.y);
    return (dx*dx) + (dy*dy);
}

template<typename CT>
CT distance(const Point<CT> &first, const Point<CT> &second) {
    return static_cast<CT>(std::pow(distanceSquared(first, second), 0.5));
}

template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
std::pair<PairType, CT> exhaustiveSearch(const PointType p[], int begin, int end)
{
    PairType pair;
    CT minSquared = std::numeric_limits<CT>::max();

    for (int i = begin; i < (end - 1); ++i) {
        const PointType &first(p[i]);

        for (int j = i + 1; j < end; ++j) {
            const PointType &second(p[j]);
            const CT distSquared = distanceSquared(second, first);
            if (distSquared < minSquared) {
                minSquared = distSquared;
                pair = std::make_pair(first, second);
            }
        }
    }

    return std::make_pair(pair, squareRoot(minSquared));
}

template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
std::pair<PairType, CT> recursiveSearch(const PointType p[], int begin, int end)
{
    const int span = end - begin;
    if (span < 64) {
        return exhaustiveSearch(p, begin, end);
    }

    const int mid = begin + (span + 1)/2;

    std::pair<PairType, CT> lresult = recursiveSearch(p, begin, mid);
    std::pair<PairType, CT> rresult = recursiveSearch(p, mid, end);

    PairType &lpair = lresult.first;
    CT &lmin = lresult.second;

    PairType &rpair = rresult.first;
    CT &rmin = rresult.second;

    {
        CT min = std::min(lmin, rmin);
        CT minSquared = min * min;
        CT upperBound = p[mid - 1].x + min;

        for (int i = mid; i < end; ++i) {
            const PointType &rp = p[i];
            if (rp.x >= upperBound) {
                break;
            } else {
                CT lowerBound = rp.x - min;
                for (int j = mid - 1; j >= begin; --j) {
                    const PointType &lp = p[j];
                    if (lp.x <= lowerBound) {
                        break;
                    } else {
                        const CT distSquared = distanceSquared(lp, rp);
                        if (distSquared < minSquared) {
                            minSquared = distSquared;
                            min = squareRoot(distSquared);
                            lmin = min;
                            lpair = std::make_pair(lp, rp);

                            upperBound = p[mid - 1].x + min;
                            lowerBound = rp.x - min;
                        }
                    }
                }
            }
        }
    }
    {
        CT min = std::min(lmin, rmin);
        CT minSquared = min * min;
        CT lowerBound = p[mid].x - min;

        for (int i = mid - 1; i >= begin; --i) {
            const PointType &lp = p[i];
            if (lp.x <= lowerBound) {
                break;
            } else {
                CT upperBound = lp.x + min;
                for (int j = mid; j < end; ++j) {
                    const PointType &rp = p[j];
                    if (rp.x >= upperBound) {
                        break;
                    } else {
                        const CT distSquared = distanceSquared(lp, rp);
                        if (distSquared < minSquared) {
                            minSquared = distSquared;
                            min = squareRoot(distSquared);
                            rmin = min;
                            rpair = std::make_pair(lp, rp);

                            lowerBound = p[mid].x - min;
                            upperBound = lp.x + min;
                        }
                    }
                }
            }
        }
    }

    return std::make_pair((rmin < lmin ? rpair : lpair), (rmin < lmin ? rmin : lmin));
}

template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
std::pair<PairType, CT> neighbourSearch(const PointType p[], int begin, int end)
{
    PairType pair;
    CT minSquared = distanceSquared(p[begin], p[begin + 1]);
    CT min = squareRoot(minSquared);

    for (int i = begin; i < (end - 1); ++i) {
        const PointType &first = p[i];
        CT upperBound = first.x + min;

        for (int j = i + 1; j < end; ++j) {
            const PointType &second = p[j];
            if (second.x >= upperBound) {
                break;
            } else {
                const CT distSquared = distanceSquared(first, second);
                if (distSquared < minSquared) {
                    minSquared = distSquared;
                    min = squareRoot(distSquared);
                    pair = std::make_pair(first, second);

                    upperBound = first.x + min;
                }
            }
        }
    }

    return std::make_pair(pair, min);
}

template<typename CT>
void timeOperations()
{
#ifdef USE_EXHAUSTIVE
    const int N = 40 * 1000;
#else
    const int N = 400 * 1000;
#endif

    std::vector<Point<CT> > points(N);
    std::vector<Point<CT> > sorted(N);

    std::default_random_engine rng(std::chrono::system_clock::to_time_t(std::chrono::system_clock::now()));
    std::uniform_real_distribution<CT> dst(0.0, 1.0);

    for (int i = 0; i < N; ++i) { points[i] = Point<CT>{ dst(rng), dst(rng) }; }

    std::pair<std::pair<Point<CT>, Point<CT> >, CT> min;

    std::chrono::system_clock::time_point tp = std::chrono::system_clock::now();
    std::chrono::system_clock::duration elapsed;

#ifdef USE_EXHAUSTIVE
    min = exhaustiveSearch(points.data(), 0, N);
    elapsed = std::chrono::system_clock::now() - tp;
    std::cout << "exhaustive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;

    tp = std::chrono::system_clock::now();
#endif

    for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
    std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
    min = recursiveSearch(sorted.data(), 0, N);
    elapsed = std::chrono::system_clock::now() - tp;
    std::cout << "recursive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;

    tp = std::chrono::system_clock::now();

    for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
    std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
    min = neighbourSearch(sorted.data(), 0, N);
    elapsed = std::chrono::system_clock::now() - tp;
    std::cout << "neighbour min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
}

int main(int,char**)
{
    timeOperations<float>();
    return 0;
}