library/geometry/closest-pair.hpp

• View this file on GitHub
• Last update: 2023-05-26 16:40:39+08:00
• Include: #include "library/geometry/closest-pair.hpp"

Depends on

• library/geometry/point.hpp

Code

#pragma once
#include <vector>
#include <algorithm>
#include <functional>
#include <limits>
#include "point.hpp"

namespace felix {

namespace geometry {

template<class T>

T closest_pair(std::vector<Point<T>> a) {
	std::sort(a.begin(), a.end(), [](Point<T> a, Point<T> b) {
		return a.x < b.x;
	});
	auto square = [&](T x) { return x * x; };
	std::function<T(int, int)> solve = [&](int l, int r) {
		if(l + 1 == r) {
			return std::numeric_limits<T>::max();
		}
		int mid = (l + r) / 2;
		T mid_x = a[mid].x;
		T ans = std::min(solve(l, mid), solve(mid, r));
		std::inplace_merge(a.begin() + l, a.begin() + mid, a.begin() + r, [](Point<T> a, Point<T> b) {
			return a.y < b.y;
		});
		std::vector<Point<T>> p;
		for(int i = l; i < r; i++) {
			if(square(a[i].x - mid_x) < ans) {
				p.push_back(a[i]);
			}
		}
		for(int i = 0; i < (int) p.size(); i++) {
			for(int j = i + 1; j < (int) p.size(); j++) {
				ans = std::min(ans, square(p[i].x - p[j].x) + square(p[i].y - p[j].y));
				if(square(p[i].y - p[j].y) > ans) {
					break;
				}
			}
		}
		return ans;
	};
	return solve(0, a.size());
}

} // namespace geometry

} // namespace felix

#line 2 "library/geometry/closest-pair.hpp"
#include <vector>
#include <algorithm>
#include <functional>
#include <limits>
#line 2 "library/geometry/point.hpp"
#include <iostream>
#include <cmath>

namespace felix {

namespace geometry {

template<class T>
struct Point {
	T x, y;

	Point(T a = 0, T b = 0) : x(a), y(b) {}
	Point(const std::pair<T, T>& p) : x(p.first), y(p.second) {}

	explicit constexpr operator std::pair<T, T>() const { return std::pair<T, T>(x, y); }

	constexpr Point& operator+=(const Point& rhs) & {
		x += rhs.x, y += rhs.y;
		return *this;
	}

	constexpr Point& operator-=(const Point& rhs) & {
		x -= rhs.x, y -= rhs.y;
		return *this;
	}

	constexpr Point& operator*=(const T& rhs) & {
		x *= rhs, y *= rhs;
		return *this;
	}

	constexpr Point& operator/=(const T& rhs) & {
		x /= rhs, y /= rhs;
		return *this;
	}

	constexpr Point operator+() const { return *this; }
	constexpr Point operator-() const { return Point(-x, -y); }
	friend constexpr Point operator+(Point lhs, Point rhs) { return lhs += rhs; }
	friend constexpr Point operator-(Point lhs, Point rhs) { return lhs -= rhs; }
	friend constexpr Point operator*(Point lhs, T rhs) { return lhs *= rhs; }
	friend constexpr Point operator/(Point lhs, T rhs) { return lhs /= rhs; }
	constexpr bool operator==(const Point& rhs) const { return x == rhs.x && y == rhs.y; }
	constexpr bool operator!=(const Point& rhs) const { return !(*this == rhs); }

	// rotate counter-clockwise
	constexpr Point rotate(T theta) const {
		T sin_t = std::sin(theta), cos_t = std::cos(theta);
		return Point(x * cos_t - y * sin_t, x * sin_t + y * cos_t);
	}

	friend constexpr T abs2(Point p) { return p.x * p.x + p.y * p.y; }
	friend constexpr long double abs(Point p) { return std::sqrt(abs2(p)); }
	friend constexpr long double angle(Point p) { return std::atan2(p.y, p.x); }
	friend constexpr T dot(Point lhs, Point rhs) { return lhs.x * rhs.x + lhs.y * rhs.y; }
	friend constexpr T cross(Point lhs, Point rhs) { return lhs.x * rhs.y - lhs.y * rhs.x; }

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

} // namespace geometry

} // namespace felix
#line 7 "library/geometry/closest-pair.hpp"

namespace felix {

namespace geometry {

template<class T>

T closest_pair(std::vector<Point<T>> a) {
	std::sort(a.begin(), a.end(), [](Point<T> a, Point<T> b) {
		return a.x < b.x;
	});
	auto square = [&](T x) { return x * x; };
	std::function<T(int, int)> solve = [&](int l, int r) {
		if(l + 1 == r) {
			return std::numeric_limits<T>::max();
		}
		int mid = (l + r) / 2;
		T mid_x = a[mid].x;
		T ans = std::min(solve(l, mid), solve(mid, r));
		std::inplace_merge(a.begin() + l, a.begin() + mid, a.begin() + r, [](Point<T> a, Point<T> b) {
			return a.y < b.y;
		});
		std::vector<Point<T>> p;
		for(int i = l; i < r; i++) {
			if(square(a[i].x - mid_x) < ans) {
				p.push_back(a[i]);
			}
		}
		for(int i = 0; i < (int) p.size(); i++) {
			for(int j = i + 1; j < (int) p.size(); j++) {
				ans = std::min(ans, square(p[i].x - p[j].x) + square(p[i].y - p[j].y));
				if(square(p[i].y - p[j].y) > ans) {
					break;
				}
			}
		}
		return ans;
	};
	return solve(0, a.size());
}

} // namespace geometry

} // namespace felix

Back to top page