This documentation is automatically generated by online-judge-tools/verification-helper
#include "geometry/closest_pair.hpp"
#pragma once
#include "point.hpp"
#include <tuple>
#include <limits>
#include <algorithm>
#include <numeric>
// 最近点対 (分割統治法)
// O(n log n)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_5_A
// return {index1, index2, distance}
template <class T> std::tuple<int, int, T> closest_pair(const std::vector<Point<T>>& p) {
const int n = (int)(p.size());
assert(n >= 2);
if (n == 2) {
return {0, 1, abs(p[0] - p[1])};
}
// may not be efficient due to indirect references ...
std::vector<int> ind(n);
std::iota(ind.begin(), ind.end(), 0);
std::sort(ind.begin(), ind.end(), [&](int i, int j) { return compare_x(p[i], p[j]); });
auto dfs = [&](auto f, int l, int r) -> std::tuple<int, int, T> {
if (r - l <= 1) return {-1, -1, std::numeric_limits<T>::max()};
const int md = (l + r) / 2;
T x = p[ind[md]].x;
// 分割統治法
auto [i1l, i2l, dl] = f(f, l, md);
auto [i1r, i2r, dr] = f(f, md, r);
int i1 = i1r, i2 = i2r;
T d = dr;
if (dl < dr) {
d = dl, i1 = i1l, i2 = i2l;
}
std::inplace_merge(ind.begin() + l, ind.begin() + md, ind.begin() + r, [&](int i, int j) { return compare_y(p[i], p[j]); });
// ind は y でソートされる
std::vector<int> near_x; // 直線 x からの距離が d 未満の頂点のインデックス
for (int i = l; i < r; i++) {
if (sign(std::abs(p[ind[i]].x - x) - d) >= 0) continue; // |p[ind[i]].x - x| >= d
const int sz = (int)(near_x.size());
// y座標との距離が d 以上になるまで繰り返す
for (int j = sz - 1; j >= 0; j--) {
Point cp = p[ind[i]] - p[near_x[j]];
if (sign(cp.y - d) >= 0) break; // cp.y >= d
T cd = abs(cp);
if (cd < d) {
d = cd, i1 = ind[i], i2 = near_x[j];
}
}
near_x.push_back(ind[i]);
}
return {i1, i2, d};
};
return dfs(dfs, 0, n);
}
#line 2 "geometry/closest_pair.hpp"
#line 2 "geometry/point.hpp"
#line 2 "geometry/geometry_template.hpp"
#include <type_traits>
// Constants (EPS, PI)
// EPS の変更は Constants<T>::set_eps(new_eps) で
template <class T> struct Constants {
static T EPS;
static void set_eps(const T e) { EPS = e; }
static constexpr T PI = 3.14159'26535'89793L;
};
template <> double Constants<double>::EPS = 1e-9;
template <> long double Constants<long double>::EPS = 1e-12;
template <> long long Constants<long long>::EPS = 0;
// 汎用関数
template <class T> inline int sign(const T x) { return x < -Constants<T>::EPS ? -1 : (x > Constants<T>::EPS ? 1 : 0); }
template <class T> inline bool equal(const T a, const T b) { return sign(a - b) == 0; }
template <class T> inline T radian_to_degree(const T r) { return r * 180.0 / Constants<T>::PI; }
template <class T> inline T degree_to_radian(const T d) { return d * Constants<T>::PI / 180.0; }
// type traits
template <class T> using is_geometry_floating_point = typename std::conditional<std::is_same<T, double>::value || std::is_same<T, long double>::value, std::true_type, std::false_type>::type;
template <class T> using is_geometry_integer = typename std::conditional<std::is_same<T, long long>::value, std::true_type, std::false_type>::type;
template <class T> using is_geometry = typename std::conditional<is_geometry_floating_point<T>::value || is_geometry_integer<T>::value, std::true_type, std::false_type>::type;
#line 4 "geometry/point.hpp"
#include <cmath>
#include <cassert>
// 点
template <class T> struct Point {
T x, y;
Point() = default;
Point(const T x, const T y) : x(x), y(y) {}
template <class U> Point(const Point<U> p) : x(p.x), y(p.y) {}
Point& operator+=(const Point& p) {
x += p.x, y += p.y;
return *this;
}
Point& operator-=(const Point& p) {
x -= p.x, y -= p.y;
return *this;
}
Point& operator*=(const Point& p) {
static_assert(is_geometry_floating_point<T>::value == true);
return *this = Point(x * p.x - y * p.y, x * p.y + y * p.x);
}
Point& operator/=(const Point& p) {
static_assert(is_geometry_floating_point<T>::value == true);
return *this = Point(x * p.x + y * p.y, -x * p.y + y * p.x) / (p.x * p.x + p.y * p.y);
}
Point& operator*=(const T k) {
x *= k, y *= k;
return *this;
}
Point& operator/=(const T k) {
static_assert(is_geometry_floating_point<T>::value == true);
x /= k, y /= k;
return *this;
}
Point operator+() const { return *this; }
Point operator-() const { return Point(-x, -y); }
friend Point operator+(const Point& a, const Point& b) { return Point(a) += b; }
friend Point operator-(const Point& a, const Point& b) { return Point(a) -= b; }
friend Point operator*(const Point& a, const Point& b) { return Point(a) *= b; }
friend Point operator/(const Point& a, const Point& b) { return Point(a) /= b; }
friend Point operator*(const Point& p, const T k) { return Point(p) *= k; }
friend Point operator/(const Point& p, const T k) { return Point(p) /= k; }
// 辞書式順序
friend bool operator<(const Point& a, const Point& b) { return a.x == b.x ? a.y < b.y : a.x < b.x; }
friend bool operator>(const Point& a, const Point& b) { return a.x == b.x ? a.y > b.y : a.x > b.x; }
friend bool operator==(const Point& a, const Point& b) { return a.x == b.x and a.y == b.y; }
// I/O
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 << ')'; }
};
// 汎用関数
// 点の一致判定
template <class T> inline bool equal(const Point<T>& a, const Point<T>& b) { return equal(a.x, b.x) and equal(a.y, b.y); }
// 内積
template <class T> inline T dot(const Point<T>& a, const Point<T>& b) { return a.x * b.x + a.y * b.y; }
// 外積
template <class T> inline T cross(const Point<T>& a, const Point<T>& b) { return a.x * b.y - a.y * b.x; }
// rad ラジアンだけ反時計回りに回転
template <class T> inline Point<T> rotate(const Point<T>& p, const T theta) {
static_assert(is_geometry_floating_point<T>::value == true);
return p * Point<T>(std::cos(theta), std::sin(theta));
}
// (x, y) の辞書式順序 (誤差許容)
template <class T> inline bool compare_x(const Point<T>& a, const Point<T>& b) { return equal(a.x, b.x) ? sign(a.y - b.y) < 0 : sign(a.x - b.x) < 0; }
// (y, x) の辞書式順序 (誤差許容)
template <class T> inline bool compare_y(const Point<T>& a, const Point<T>& b) { return equal(a.y, b.y) ? sign(a.x - b.x) < 0 : sign(a.y - b.y) < 0; }
// 整数のまま行う偏角ソート
// 無限の精度をもつ arg(p) = atan2(y, x) で比較し, 同じ場合は norm(p) で比較 (atan2(0, 0) = 0 とする)
// 基本的に (-PI, PI] でソートされ, 点 (0, 0) は (-PI, 0) と [0, PI] の間に入る
// https://ngtkana.hatenablog.com/entry/2021/11/13/202103
// https://judge.yosupo.jp/problem/sort_points_by_argument
template <class T> inline bool compare_atan2(const Point<T>& a, const Point<T>& b) {
static_assert(is_geometry_integer<T>::value == true);
if ((Point<T>(a.y, -a.x) > Point<T>(0, 0)) == (Point<T>(b.y, -b.x) > Point<T>(0, 0))) { // a, b in (-PI, 0] or a, b in (0, PI]
if (a.x * b.y != a.y * b.x) return a.x * b.y > a.y * b.x; // cross(a, b) != 0
return a == Point<T>(0, 0) ? b.x > 0 and b.y == 0 : (b == Point<T>(0, 0) ? a.y < 0 : norm(a) < norm(b));
}
return Point<T>(a.y, -a.x) < Point<T>(b.y, -b.x);
}
// 絶対値の 2 乗
template <class T> inline T norm(const Point<T>& p) { return p.x * p.x + p.y * p.y; }
// 絶対値
template <class T> inline T abs(const Point<T>& p) {
static_assert(is_geometry_floating_point<T>::value == true);
return std::sqrt(norm(p));
}
// 偏角
template <class T> inline T arg(const Point<T>& p) {
static_assert(is_geometry_floating_point<T>::value == true);
return std::atan2(p.y, p.x); // (-PI, PI]
}
// 共役複素数 (x 軸について対象な点)
template <class T> inline Point<T> conj(const Point<T>& p) { return Point(p.x, -p.y); }
// 極座標系 -> 直交座標系 (絶対値が rho で偏角が theta ラジアン)
template <class T> inline Point<T> polar(const T rho, const T theta = T(0)) {
static_assert(is_geometry_floating_point<T>::value == true);
assert(sign(rho) >= 0);
return Point<T>(std::cos(theta), std::sin(theta)) * rho;
}
// ccw の戻り値
enum class Ccw {
COUNTER_CLOCKWISE = 1, // a->b->c 反時計回り
CLOCKWISE = -1, // a->b->c 時計回り
ONLINE_BACK = 2, // c->a->b 直線
ONLINE_FRONT = -2, // a->b->c 直線
ON_SEGMENT = 0 // a->c->b 直線
};
// 点 a, b, c の位置関係
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C
template <class T> Ccw ccw(const Point<T>& a, const Point<T>& b, const Point<T>& c) {
Point<T> ab = b - a;
Point<T> ac = c - a;
if (sign(cross(ab, ac)) == 1) return Ccw::COUNTER_CLOCKWISE;
if (sign(cross(ab, ac)) == -1) return Ccw::CLOCKWISE;
if (sign(dot(ab, ac)) == -1) return Ccw::ONLINE_BACK;
if (sign(norm(ab) - norm(ac)) == -1) return Ccw::ONLINE_FRONT;
return Ccw::ON_SEGMENT;
}
// 線分 a -> b から 線分 a -> c までの角度 (ラジアンで -PI より大きく PI 以下)
template <class T> T get_angle(const Point<T>& a, const Point<T>& b, const Point<T>& c) {
Point<T> ab = b - a;
Point<T> ac = c - a;
// a-bが x 軸になるように回転
ac *= conj(ab) / norm(ab);
return arg(ac); // (-PI, PI]
}
#line 4 "geometry/closest_pair.hpp"
#include <tuple>
#include <limits>
#include <algorithm>
#include <numeric>
// 最近点対 (分割統治法)
// O(n log n)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_5_A
// return {index1, index2, distance}
template <class T> std::tuple<int, int, T> closest_pair(const std::vector<Point<T>>& p) {
const int n = (int)(p.size());
assert(n >= 2);
if (n == 2) {
return {0, 1, abs(p[0] - p[1])};
}
// may not be efficient due to indirect references ...
std::vector<int> ind(n);
std::iota(ind.begin(), ind.end(), 0);
std::sort(ind.begin(), ind.end(), [&](int i, int j) { return compare_x(p[i], p[j]); });
auto dfs = [&](auto f, int l, int r) -> std::tuple<int, int, T> {
if (r - l <= 1) return {-1, -1, std::numeric_limits<T>::max()};
const int md = (l + r) / 2;
T x = p[ind[md]].x;
// 分割統治法
auto [i1l, i2l, dl] = f(f, l, md);
auto [i1r, i2r, dr] = f(f, md, r);
int i1 = i1r, i2 = i2r;
T d = dr;
if (dl < dr) {
d = dl, i1 = i1l, i2 = i2l;
}
std::inplace_merge(ind.begin() + l, ind.begin() + md, ind.begin() + r, [&](int i, int j) { return compare_y(p[i], p[j]); });
// ind は y でソートされる
std::vector<int> near_x; // 直線 x からの距離が d 未満の頂点のインデックス
for (int i = l; i < r; i++) {
if (sign(std::abs(p[ind[i]].x - x) - d) >= 0) continue; // |p[ind[i]].x - x| >= d
const int sz = (int)(near_x.size());
// y座標との距離が d 以上になるまで繰り返す
for (int j = sz - 1; j >= 0; j--) {
Point cp = p[ind[i]] - p[near_x[j]];
if (sign(cp.y - d) >= 0) break; // cp.y >= d
T cd = abs(cp);
if (cd < d) {
d = cd, i1 = ind[i], i2 = near_x[j];
}
}
near_x.push_back(ind[i]);
}
return {i1, i2, d};
};
return dfs(dfs, 0, n);
}