cp-library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub rniya/cp-library

:heavy_check_mark: Closest Pair
(src/geometry/closest_pair.hpp)

Depends on

Verified with

Code

#pragma once
#include <algorithm>
#include <limits>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>
#include "Point.hpp"

namespace geometry {

template <typename T> std::tuple<int, int, T> closest_pair(std::vector<Point<T>> P) {
    int n = P.size();
    if (n == 1) return {0, 0, 0};
    std::vector<int> ord(n);
    std::iota(begin(ord), end(ord), 0);
    std::sort(begin(ord), end(ord), [&](int p, int q) { return P[p] < P[q]; });
    T inf = std::numeric_limits<T>::max(), min_dist = inf;
    std::pair<int, int> ans{-1, -1};

    auto dfs = [&](auto self, int l, int r) -> T {
        if (r - l <= 1) return inf;
        int mid = (l + r) >> 1;
        auto x_mid = P[ord[mid]].x;
        auto res = std::min(self(self, l, mid), self(self, mid, r));
        std::inplace_merge(begin(ord) + l, begin(ord) + mid, begin(ord) + r,
                           [&](int p, int q) { return P[p].y < P[q].y; });
        std::vector<int> tmp;
        for (int i = l; i < r; i++) {
            if (sgn(std::abs(P[ord[i]].x - x_mid) - res) >= 0) continue;
            for (int j = int(tmp.size()) - 1; j >= 0; j--) {
                auto d = P[ord[i]] - P[tmp[j]];
                if (sgn(d.y - res) >= 0) break;
                if (d.norm() < res) {
                    res = d.norm();
                    if (d.norm() < min_dist) {
                        min_dist = d.norm();
                        ans = {ord[i], tmp[j]};
                    }
                }
            }
            tmp.emplace_back(ord[i]);
        }
        return res;
    };
    dfs(dfs, 0, n);
    return {ans.first, ans.second, min_dist};
}

}  // namespace geometry
#line 2 "src/geometry/closest_pair.hpp"
#include <algorithm>
#include <limits>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>
#line 2 "src/geometry/Point.hpp"
#include <cassert>
#include <cmath>
#include <iostream>
#include <type_traits>

namespace geometry {

template <typename T> struct Point {
    static T EPS;

    static void set_eps(T eps) { EPS = eps; }

    T x, y;

    Point() {}

    Point(T x, T y) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }

    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }

    Point operator*(T t) const { return Point(x * t, y * t); }

    Point operator/(T t) const { return Point(x / t, y / t); }

    bool operator==(const Point& p) const { return x == p.x and y == p.y; }

    bool operator!=(const Point& p) const { return not((*this) == p); }

    bool operator<(const Point& p) const { return x != p.x ? x < p.x : y < p.y; }

    friend std::istream& operator>>(std::istream& is, Point& p) { return is >> p.x >> p.y; }

    friend std::ostream& operator<<(std::ostream& os, const Point& p) { return os << p.x << ' ' << p.y; }

    T norm() { return std::sqrt(x * x + y * y); }

    T norm2() { return x * x + y * y; }

    T arg() { return std::atan2(y, x); }

    T dot(const Point& p) { return x * p.x + y * p.y; }

    T det(const Point& p) { return x * p.y - y * p.x; }

    Point perp() { return Point(-y, x); }

    Point unit() { return *this / norm(); }

    Point normal() { return perp().unit(); }

    Point rotate(T rad) { return Point(std::cos(rad) * x - std::sin(rad) * y, std::sin(rad) * x + std::cos(rad) * y); }
};

template <> double Point<double>::EPS = 1e-9;
template <> long double Point<long double>::EPS = 1e-12;
template <> int Point<int>::EPS = 0;
template <> long long Point<long long>::EPS = 0;

template <typename T> int sgn(T x) { return x < -Point<T>::EPS ? -1 : x > Point<T>::EPS ? 1 : 0; }

}  // namespace geometry
#line 9 "src/geometry/closest_pair.hpp"

namespace geometry {

template <typename T> std::tuple<int, int, T> closest_pair(std::vector<Point<T>> P) {
    int n = P.size();
    if (n == 1) return {0, 0, 0};
    std::vector<int> ord(n);
    std::iota(begin(ord), end(ord), 0);
    std::sort(begin(ord), end(ord), [&](int p, int q) { return P[p] < P[q]; });
    T inf = std::numeric_limits<T>::max(), min_dist = inf;
    std::pair<int, int> ans{-1, -1};

    auto dfs = [&](auto self, int l, int r) -> T {
        if (r - l <= 1) return inf;
        int mid = (l + r) >> 1;
        auto x_mid = P[ord[mid]].x;
        auto res = std::min(self(self, l, mid), self(self, mid, r));
        std::inplace_merge(begin(ord) + l, begin(ord) + mid, begin(ord) + r,
                           [&](int p, int q) { return P[p].y < P[q].y; });
        std::vector<int> tmp;
        for (int i = l; i < r; i++) {
            if (sgn(std::abs(P[ord[i]].x - x_mid) - res) >= 0) continue;
            for (int j = int(tmp.size()) - 1; j >= 0; j--) {
                auto d = P[ord[i]] - P[tmp[j]];
                if (sgn(d.y - res) >= 0) break;
                if (d.norm() < res) {
                    res = d.norm();
                    if (d.norm() < min_dist) {
                        min_dist = d.norm();
                        ans = {ord[i], tmp[j]};
                    }
                }
            }
            tmp.emplace_back(ord[i]);
        }
        return res;
    };
    dfs(dfs, 0, n);
    return {ans.first, ans.second, min_dist};
}

}  // namespace geometry
Back to top page