幾何の点・ベクトル・偏角順

geometry/point.hpp · GitHub で見る

共通ページ

#include "lib/geometry/point.hpp"
using namespace geometry;

IPoint a{0, 0}, b{3, 0}, c{1, 2};
int turn = orientation(a, b, c); // 1
auto squared = distance_squared(b, c); // __int128_t の 8
FPoint p = rotate(b, pi / 2);

geometry 名前空間を使います。Point<T> は座標 x,y を持ち、T は符号付き組み込み整数型または浮動小数点型です。IPoint は Point<std::int64_t>、FPoint は Point<long double> です。Real は long double、pi は円周率、epsilon は 1e-12L です。

整数の内積・外積・距離二乗は __int128_t、実数の場合は Real で計算します。この返り値型を wide_t<T> として使えます。integer<T> は組み込み整数型かを表します。 cross(a,b,c) と distance_squared(a,b) は、座標を拡張してから減算します。Point<T> 自体の四則演算は T のままで、整数の除算は0方向に丸めます。

操作	内容・返り値	計算量
Point<T>() / Point<T>(x,y) / Point<T>(p)	点の構築。既定値は原点。異なる座標型からの変換は明示的	\(O(1)\)
p.x, p.y	座標	\(O(1)\)
p == q, p < q	完全一致、x・y の辞書順。EPS は使用しない	\(O(1)\)
p + q, p - q, -p	ベクトルの加減算	\(O(1)\)
p * c, c * p, p / c	座標型 T のスカラーとの演算。除数は0以外	\(O(1)\)
p += q, p -= q, p *= c, p /= c	対応する代入演算	\(O(1)\)
in >> p, out << p	x y 形式の入出力。ストリームが座標型を扱える場合	\(O(1)\)
sign(x,eps = epsilon)	負・0・正を -1,0,1 で返す。整数では EPS を無視	\(O(1)\)
dot(a,b)	内積	\(O(1)\)
cross(a,b)	外積	\(O(1)\)
cross(a,b,c)	(b-a) と (c-a) の外積	\(O(1)\)
norm2(p)	原点からの距離二乗	\(O(1)\)
distance_squared(a,b)	2点間の距離二乗	\(O(1)\)
orientation(a,b,c,eps = epsilon)	a→b から見て c が左・同一直線・右なら 1,0,-1	\(O(1)\)
ccw(a,b,c,eps = epsilon)	下記の5種類の位置関係	\(O(1)\)
to_real(p)	FPoint に変換	\(O(1)\)
length(p), distance(a,b)	ベクトルの長さ、2点間の距離。Real を返す	\(O(1)\)
almost_equal(a,b,eps = epsilon)	実数は各座標差が EPS 以下、整数は完全一致	\(O(1)\)
rotate90(p)	原点を中心に反時計回りに90度回転。型を保つ	\(O(1)\)
argument(p)	偏角。ラジアン、範囲は -π から π	\(O(1)\)
angle(a,b)	a から b への符号付き角度。両方0以外とする	\(O(1)\)
rotate(p,theta)	原点を中心に反時計回りに theta ラジアン回転。FPoint	\(O(1)\)
unit(p)	長さ1の optional<FPoint>。零ベクトルなら nullopt	\(O(1)\)
PolarLess<T>{}(a,b)	0以上2π未満の偏角順。原点が最初、同方向は距離順	\(O(1)\)
polar_sort(points)	vector<Point<T>>& を PolarLess<T> の順に並べ替える。返り値は void。空列・重複を扱える	\(O(n\log n)\)

ccw(a,b,c) は、左が 1、右が -1、a より後方が 2、b より先が -2、線分上が 0 です。a と b が同じ場合、その点自身は 0、それ以外は -2 です。

偏角ソートは polar_sort(points); とします。比較に EPS を入れると整合的なソート順にならないため、辞書順・偏角順とも EPS は使いません。

整数は全ての中間値が符号付き128 bitに収まる範囲で使います。例えば各座標の絶対値が \(10^{18}\) 以下なら、ここでの外積・距離二乗は収まります。点同士の加減算などは、別途 T に収まる必要があります。実数は有限値で、中間値も有限となる範囲とします。 実数の sign・orientation・ccw の EPS は、比較する値に対する絶対許容誤差です。orientation では外積の単位になります。必要なら引数で調整し、完全な比較には eps = 0 を指定します。浮動小数点の丸め誤差そのものがなくなるわけではありません。

polar_sort は \(O(n\log n)\) 時間・\(O(n)\) 追加空間です。それ以外の各操作の時間・追加空間は \(O(1)\) です。

依存ライブラリ

• 64 bit キーの安定基数ソート

検証コード

• test/dynamic_convex_hull.test.cpp
• test/geometry_argument_sort.test.cpp
• test/geometry_circle_line.test.cpp
• test/geometry_circle_polygon.test.cpp
• test/geometry_closest_pair.test.cpp
• test/geometry_convex_cut.test.cpp
• test/geometry_convex_hull.test.cpp
• test/geometry_furthest_pair.test.cpp
• test/geometry_polygon_contains.test.cpp
• test/geometry_projection.test.cpp
• test/geometry_segment_distance.test.cpp

コード

#ifndef GEOMETRY_POINT_HPP
#define GEOMETRY_POINT_HPP

#include "../components/radix_sort.hpp"
#include <cmath>
#include <cstdint>
#include <istream>
#include <numbers>
#include <optional>
#include <ostream>
#include <type_traits>

namespace geometry {
using Real = long double;
inline constexpr Real epsilon = 1e-12L;
inline constexpr Real pi = std::numbers::pi_v<Real>;
enum class Location { outside = -1, boundary = 0, inside = 1 };

template<class T>
inline constexpr bool integer = std::is_integral_v<T> || std::is_same_v<T, __int128_t> || std::is_same_v<T, __uint128_t>;
template<class T>
using wide_t = std::conditional_t<integer<T>, __int128_t, Real>;

template<class T>
struct Point {
	static_assert((std::is_arithmetic_v<T> && std::is_signed_v<T>) || std::is_same_v<T, __int128_t>);
	T x{}, y{};
	constexpr Point() = default;
	constexpr Point(T x, T y) : x(x), y(y) {}
	template<class U> explicit constexpr Point(const Point<U>& p) : x(static_cast<T>(p.x)), y(static_cast<T>(p.y)) {}
	constexpr bool operator==(const Point&) const = default;
	constexpr bool operator<(const Point& p) const { return x < p.x || (x == p.x && y < p.y); }
	constexpr Point operator-() const { return {-x, -y}; }
	constexpr Point& operator+=(const Point& p) { x += p.x; y += p.y; return *this; }
	constexpr Point& operator-=(const Point& p) { x -= p.x; y -= p.y; return *this; }
	constexpr Point& operator*=(T c) { x *= c; y *= c; return *this; }
	constexpr Point& operator/=(T c) { x /= c; y /= c; return *this; }
	friend constexpr Point operator+(Point p, const Point& q) { p += q; return p; }
	friend constexpr Point operator-(Point p, const Point& q) { p -= q; return p; }
	friend constexpr Point operator*(Point p, T c) { p *= c; return p; }
	friend constexpr Point operator*(T c, Point p) { p *= c; return p; }
	friend constexpr Point operator/(Point p, T c) { p /= c; return p; }
	friend std::istream& operator>>(std::istream& in, Point& p) { return in >> p.x >> p.y; }
	friend std::ostream& operator<<(std::ostream& out, const Point& p) { return out << p.x << ' ' << p.y; }
};

using IPoint = Point<std::int64_t>;
using FPoint = Point<Real>;

template<class T>
constexpr int sign(T value, Real eps = epsilon) {
	if constexpr (integer<T>) return (value > 0) - (value < 0);
	else return (value > eps) - (value < -eps);
}
template<class T>
constexpr wide_t<T> dot(Point<T> a, Point<T> b) {
	return wide_t<T>(a.x) * b.x + wide_t<T>(a.y) * b.y;
}
template<class T>
constexpr wide_t<T> cross(Point<T> a, Point<T> b) {
	return wide_t<T>(a.x) * b.y - wide_t<T>(a.y) * b.x;
}
template<class T>
constexpr wide_t<T> cross(Point<T> a, Point<T> b, Point<T> c) {
	using W = wide_t<T>;
	return (W(b.x) - a.x) * (W(c.y) - a.y) - (W(b.y) - a.y) * (W(c.x) - a.x);
}
template<class T>
constexpr wide_t<T> norm2(Point<T> p) { return dot(p, p); }
template<class T>
constexpr wide_t<T> distance_squared(Point<T> a, Point<T> b) {
	using W = wide_t<T>;
	W x = W(a.x) - b.x, y = W(a.y) - b.y;
	return x * x + y * y;
}
template<class T>
constexpr int orientation(Point<T> a, Point<T> b, Point<T> c, Real eps = epsilon) { return sign(cross(a, b, c), eps); }
template<class T>
constexpr int ccw(Point<T> a, Point<T> b, Point<T> c, Real eps = epsilon) {
	if (int turn = orientation(a, b, c, eps)) return turn;
	using W = wide_t<T>;
	W x = W(b.x) - a.x, y = W(b.y) - a.y, u = W(c.x) - a.x, v = W(c.y) - a.y;
	if (sign(x * u + y * v, eps) < 0) return 2;
	if (sign(u * u + v * v - x * x - y * y, eps) > 0) return -2;
	return 0;
}
template<class T>
constexpr FPoint to_real(Point<T> p) { return FPoint(p); }
template<class T>
Real length(Point<T> p) { return std::hypot(Real(p.x), Real(p.y)); }
template<class T>
Real distance(Point<T> a, Point<T> b) {
	return std::hypot(Real(wide_t<T>(a.x) - b.x), Real(wide_t<T>(a.y) - b.y));
}
template<class T>
constexpr bool almost_equal(Point<T> a, Point<T> b, Real eps = epsilon) {
	if constexpr (integer<T>) return a == b;
	else return std::abs(Real(a.x) - b.x) <= eps && std::abs(Real(a.y) - b.y) <= eps;
}
template<class T>
constexpr Point<T> rotate90(Point<T> p) { return {-p.y, p.x}; }
template<class T>
Real argument(Point<T> p) { return std::atan2(Real(p.y), Real(p.x)); }
template<class T>
Real angle(Point<T> a, Point<T> b) { return std::atan2(Real(cross(a, b)), Real(dot(a, b))); }
template<class T>
FPoint rotate(Point<T> p, Real theta) {
	Real c = std::cos(theta), s = std::sin(theta);
	return {Real(p.x) * c - Real(p.y) * s, Real(p.x) * s + Real(p.y) * c};
}
template<class T>
std::optional<FPoint> unit(Point<T> p) {
	Real n = length(p);
	if (!n) return std::nullopt;
	return to_real(p) / n;
}
template<class T>
struct PolarLess {
	bool operator()(Point<T> a, Point<T> b) const {
		if (a == Point<T>{} || b == Point<T>{}) return a == Point<T>{} && b != Point<T>{};
		auto half = [](Point<T> p) { return p.y < 0 || (p.y == 0 && p.x < 0); };
		if (half(a) != half(b)) return half(a) < half(b);
		if constexpr (integer<T> && sizeof(T) <= 4) {
			auto c = std::int64_t(a.x) * b.y - std::int64_t(a.y) * b.x;
			if (c) return c > 0;
			auto norm = [](Point<T> p) { return std::uint64_t(std::int64_t(p.x) * p.x) + std::uint64_t(std::int64_t(p.y) * p.y); };
			return norm(a) < norm(b);
		} else {
			if (auto c = cross(a, b); c != 0) return c > 0;
			return norm2(a) < norm2(b);
		}
	}
};

template<class T>
void polar_sort(std::vector<Point<T>> &points) {
	if constexpr (integer<T> && sizeof(T) <= 4) {
		if (points.size() >= 256) {
			struct Entry { std::uint64_t key; Point<T> point; };
			std::vector<Entry> order;
			order.reserve(points.size());
			constexpr std::uint64_t scale = (std::uint64_t(1) << 61) - 1;
			for (auto p : points) {
				std::int64_t x = p.x, y = p.y;
				std::uint64_t a = x < 0 ? -x : x, b = y < 0 ? -y : y;
				unsigned quadrant = y < 0 ? (x < 0 ? 2 : 3) : (x < 0 ? 1 : 0);
				unsigned octant = quadrant * 2 + ((a < b) ^ (quadrant & 1));
				std::uint64_t ratio = a || b ? std::min(scale, std::uint64_t(double(std::min(a, b)) / double(std::max(a, b)) * double(scale))) : 0;
				if (octant & 1) ratio = scale - ratio;
				order.push_back({(std::uint64_t(octant) << 61) | ratio, p});
			}
			stable_radix_sort<11>(order, [](const Entry &entry) { return entry.key; });
			for (std::size_t first = 0; first < order.size();) {
				std::size_t last = first + 1;
				while (last < order.size() && order[last].key == order[first].key) ++last;
				auto less = [](const Entry &a, const Entry &b) { return PolarLess<T>{}(a.point, b.point); };
				if (last - first <= 8 || !std::is_sorted(order.begin() + first, order.begin() + last, less))
					std::sort(order.begin() + first, order.begin() + last, less);
				first = last;
			}
			for (std::size_t i = 0; i < points.size(); ++i) points[i] = order[i].point;
			return;
		}
	}
	std::sort(points.begin(), points.end(), PolarLess<T>{});
}

}

#endif

Back to top page