NicheLibrary
This documentation is automatically generated by NotLeonian/competitive-verifier (forked from competitive-verifier/competitive-verifier)
verify/standalone-line-convex-polygon-intersection.test.cpp
- View this file on GitHub
- Last update: 2026-07-11 16:54:42+09:00
Depends on
- 凸包と直線の交点（ $O(\log N)$ ） (geometry/line-convex-polygon-intersection.hpp)
- 128 bit 整数型 (internal/int128.hpp)
Code
// competitive-verifier: STANDALONE
#include <cassert>
#include <cmath>
#include <complex>
#include <vector>
#include "../geometry/line-convex-polygon-intersection.hpp"
int main() {
using Point = std::complex<long long>;
using NicheLibrary::Int128;
{
std::vector<Point> hull = {Point(2, 3)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 3), Point(1, 3));
assert(result.size() == 1);
assert(result[0].x_numerator == Int128(2));
assert(result[0].y_numerator == Int128(3));
assert(result[0].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(2, 3)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 4), Point(1, 4));
assert(result.empty());
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 2), Point(1, 2));
assert(result.size() == 1);
assert(result[0].x_numerator == Int128(2));
assert(result[0].y_numerator == Int128(2));
assert(result[0].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 0), Point(1, 1));
assert(result.size() == 2);
assert(result[0].denominator == Int128(1));
assert(result[1].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 5), Point(1, 5));
assert(result.empty());
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 0), Point(0, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 2), Point(1, 2));
assert(result.size() == 2);
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 0), Point(0, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 1), Point(3, 2));
assert(result.size() == 2);
auto is_point = [&](int index, Int128 x_numerator, Int128 y_numerator,
Int128 denominator) {
return result[index].x_numerator == x_numerator &&
result[index].y_numerator == y_numerator &&
result[index].denominator == denominator;
};
const bool expected_order =
is_point(0, Int128(0), Int128(1), Int128(1)) &&
is_point(1, Int128(9), Int128(7), Int128(4));
const bool reverse_order =
is_point(1, Int128(0), Int128(1), Int128(1)) &&
is_point(0, Int128(9), Int128(7), Int128(4));
assert(expected_order || reverse_order);
const int rational_index = result[0].denominator == Int128(4) ? 0 : 1;
const auto real_point =
result[rational_index]
.template to_point<std::complex<long double>>();
assert(std::abs(real_point.real() - 2.25L) < 1e-15L);
assert(std::abs(real_point.imag() - 1.75L) < 1e-15L);
}
return 0;
}
#line 1 "verify/standalone-line-convex-polygon-intersection.test.cpp"
// competitive-verifier: STANDALONE
#include <cassert>
#include <cmath>
#include <complex>
#include <vector>
#line 1 "geometry/line-convex-polygon-intersection.hpp"
// 凸包と直線の共通部分を求める。
// 凸包が空、1 点、線分、面積正の狭義凸多角形のどれであっても扱える。
// 共通部分が空集合ならば 0 個、1 点ならば 1 個、線分ならばその端点 2 個を返す。
// 前提: line_a != line_b。
// 前提: hull のサイズが 2 ならば 2 点は相異なる。
// 前提: hull のサイズが 3 以上ならば hull は反時計回りであり、
// 連続する 3 頂点が一直線に並んでおらず、面積が正の狭義凸多角形である。
// 座標型が整数ならば整数演算で処理し、交点を有理表現で返す。
// 標準の 64 bit 以下の整数座標では、既定で 128 bit 整数型を内部計算に用いる。
// 計算途中に必要な値が内部計算の型に収まることを仮定する。
// 計算量 O(log N)。
#line 19 "geometry/line-convex-polygon-intersection.hpp"
#include <limits>
#include <type_traits>
#include <utility>
#line 23 "geometry/line-convex-polygon-intersection.hpp"
#line 1 "internal/int128.hpp"
// 128 bit 符号なし整数型 UInt128 と 128 bit 符号付き整数型 Int128 を提供する。
// GCC 拡張には依存しない。
#include <bit>
#line 9 "internal/int128.hpp"
#include <cstdint>
#include <istream>
#line 12 "internal/int128.hpp"
#include <ostream>
#include <string>
#line 15 "internal/int128.hpp"
namespace NicheLibrary {
class UInt128 {
public:
constexpr UInt128() = default;
template <
class T,
std::enable_if_t<std::is_integral_v<T> && std::is_signed_v<T>, int> = 0>
constexpr UInt128(T value)
: high_(value < 0 ? ~std::uint64_t{} : std::uint64_t{}),
low_(static_cast<std::uint64_t>(value)) {}
template <class T,
std::enable_if_t<std::is_integral_v<T> && !std::is_signed_v<T>,
int> = 0>
constexpr UInt128(T value)
: high_(0), low_(static_cast<std::uint64_t>(value)) {}
static constexpr UInt128 from_words(std::uint64_t high, std::uint64_t low) {
UInt128 value;
value.high_ = high;
value.low_ = low;
return value;
}
constexpr std::uint64_t high() const { return high_; }
constexpr std::uint64_t low() const { return low_; }
template <class T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
explicit constexpr operator T() const {
return static_cast<T>(low_);
}
explicit constexpr operator bool() const { return high_ != 0 || low_ != 0; }
explicit constexpr operator long double() const {
if (high_ == 0) {
return static_cast<long double>(low_);
}
return static_cast<long double>(high_) * 0x1p64L +
static_cast<long double>(low_);
}
friend constexpr bool operator==(UInt128 lhs, UInt128 rhs) {
return lhs.high_ == rhs.high_ && lhs.low_ == rhs.low_;
}
friend constexpr bool operator!=(UInt128 lhs, UInt128 rhs) {
return !(lhs == rhs);
}
friend constexpr bool operator<(UInt128 lhs, UInt128 rhs) {
if (lhs.high_ != rhs.high_) {
return lhs.high_ < rhs.high_;
}
return lhs.low_ < rhs.low_;
}
friend constexpr bool operator>(UInt128 lhs, UInt128 rhs) {
return rhs < lhs;
}
friend constexpr bool operator<=(UInt128 lhs, UInt128 rhs) {
return !(rhs < lhs);
}
friend constexpr bool operator>=(UInt128 lhs, UInt128 rhs) {
return !(lhs < rhs);
}
constexpr UInt128 operator+() const { return *this; }
constexpr UInt128 operator-() const { return UInt128{} - *this; }
constexpr UInt128 &operator+=(UInt128 rhs) {
const std::uint64_t old_low = low_;
low_ += rhs.low_;
high_ += rhs.high_ + (low_ < old_low ? 1 : 0);
return *this;
}
constexpr UInt128 &operator-=(UInt128 rhs) {
const std::uint64_t old_low = low_;
low_ -= rhs.low_;
high_ -= rhs.high_ + (old_low < rhs.low_ ? 1 : 0);
return *this;
}
constexpr UInt128 &operator*=(UInt128 rhs) {
*this = *this * rhs;
return *this;
}
constexpr UInt128 &operator/=(UInt128 rhs) {
*this = *this / rhs;
return *this;
}
constexpr UInt128 &operator%=(UInt128 rhs) {
*this = *this % rhs;
return *this;
}
friend constexpr UInt128 operator+(UInt128 lhs, UInt128 rhs) {
lhs += rhs;
return lhs;
}
friend constexpr UInt128 operator-(UInt128 lhs, UInt128 rhs) {
lhs -= rhs;
return lhs;
}
friend constexpr UInt128 operator*(UInt128 lhs, UInt128 rhs) {
const UInt128 p00 = multiply_u64(lhs.low_, rhs.low_);
return from_words(
p00.high_ + lhs.low_ * rhs.high_ + lhs.high_ * rhs.low_, p00.low_);
}
static constexpr void div_mod(UInt128 lhs, UInt128 rhs, UInt128 "ient,
UInt128 &remainder) {
assert(rhs != UInt128{});
div_mod_unchecked(lhs, rhs, quotient, remainder);
}
friend constexpr UInt128 operator/(UInt128 lhs, UInt128 rhs) {
UInt128 quotient;
UInt128 remainder;
div_mod(lhs, rhs, quotient, remainder);
return quotient;
}
friend constexpr UInt128 operator%(UInt128 lhs, UInt128 rhs) {
UInt128 quotient;
UInt128 remainder;
div_mod(lhs, rhs, quotient, remainder);
return remainder;
}
private:
friend class Int128;
static constexpr std::uint64_t word_base() {
return std::uint64_t{1} << 32;
}
static constexpr std::uint64_t word_mask() { return word_base() - 1; }
static constexpr void to_words(UInt128 value, std::uint32_t words[4]) {
words[0] = static_cast<std::uint32_t>(value.low_ & word_mask());
words[1] = static_cast<std::uint32_t>(value.low_ >> 32);
words[2] = static_cast<std::uint32_t>(value.high_ & word_mask());
words[3] = static_cast<std::uint32_t>(value.high_ >> 32);
}
static constexpr UInt128 from_32bit_words(const std::uint32_t words[4]) {
return UInt128::from_words(
(static_cast<std::uint64_t>(words[3]) << 32) | words[2],
(static_cast<std::uint64_t>(words[1]) << 32) | words[0]);
}
static constexpr int word_length(const std::uint32_t words[4]) {
int length = 4;
while (length > 0 && words[length - 1] == 0) {
--length;
}
return length;
}
// divisor ≠ 0 が呼び出し元で保証されていることを仮定する。
static constexpr UInt128 div_mod_uint32(UInt128 value,
std::uint32_t divisor,
std::uint32_t &remainder) {
const std::uint32_t words[4] = {
static_cast<std::uint32_t>(value.low_ & word_mask()),
static_cast<std::uint32_t>(value.low_ >> 32),
static_cast<std::uint32_t>(value.high_ & word_mask()),
static_cast<std::uint32_t>(value.high_ >> 32),
};
std::uint32_t quotient_words[4] = {};
std::uint64_t rem = 0;
for (int i = 3; i >= 0; --i) {
const std::uint64_t current = (rem << 32) | words[i];
quotient_words[i] = static_cast<std::uint32_t>(current / divisor);
rem = current % divisor;
}
remainder = static_cast<std::uint32_t>(rem);
return from_32bit_words(quotient_words);
}
static constexpr void shift_left_words(const std::uint32_t source[4],
int source_length, int shift,
std::uint32_t result[5]) {
std::uint64_t carry = 0;
for (int i = 0; i < source_length; ++i) {
const std::uint64_t current =
(static_cast<std::uint64_t>(source[i]) << shift) | carry;
result[i] = static_cast<std::uint32_t>(current & word_mask());
carry = current >> 32;
}
result[source_length] = static_cast<std::uint32_t>(carry);
}
static constexpr void shift_right_words(const std::uint32_t source[5],
int source_length, int shift,
std::uint32_t result[4]) {
if (shift == 0) {
for (int i = 0; i < source_length; ++i) {
result[i] = source[i];
}
return;
}
for (int i = 0; i < source_length; ++i) {
const std::uint32_t next = i + 1 < 5 ? source[i + 1] : 0;
result[i] = static_cast<std::uint32_t>(
(source[i] >> shift) |
(static_cast<std::uint64_t>(next) << (32 - shift)));
}
}
static constexpr bool subtract_product(std::uint32_t words[5], int offset,
const std::uint32_t divisor[5],
int divisor_length,
std::uint32_t multiplier) {
std::uint64_t borrow = 0;
std::uint64_t carry = 0;
for (int i = 0; i < divisor_length; ++i) {
const std::uint64_t product =
static_cast<std::uint64_t>(multiplier) * divisor[i] + carry;
carry = product >> 32;
const std::uint64_t subtrahend = (product & word_mask()) + borrow;
if (static_cast<std::uint64_t>(words[offset + i]) >= subtrahend) {
words[offset + i] = static_cast<std::uint32_t>(
static_cast<std::uint64_t>(words[offset + i]) - subtrahend);
borrow = 0;
} else {
words[offset + i] = static_cast<std::uint32_t>(
word_base() + words[offset + i] - subtrahend);
borrow = 1;
}
}
const std::uint64_t subtrahend = carry + borrow;
if (static_cast<std::uint64_t>(words[offset + divisor_length]) >=
subtrahend) {
words[offset + divisor_length] = static_cast<std::uint32_t>(
static_cast<std::uint64_t>(words[offset + divisor_length]) -
subtrahend);
return false;
}
words[offset + divisor_length] = static_cast<std::uint32_t>(
word_base() + words[offset + divisor_length] - subtrahend);
return true;
}
static constexpr void add_back(std::uint32_t words[5], int offset,
const std::uint32_t divisor[5],
int divisor_length) {
std::uint64_t carry = 0;
for (int i = 0; i < divisor_length; ++i) {
const std::uint64_t sum =
static_cast<std::uint64_t>(words[offset + i]) + divisor[i] +
carry;
words[offset + i] = static_cast<std::uint32_t>(sum & word_mask());
carry = sum >> 32;
}
words[offset + divisor_length] = static_cast<std::uint32_t>(
static_cast<std::uint64_t>(words[offset + divisor_length]) + carry);
}
static constexpr void div_mod_knuth(const std::uint32_t lhs_words[4],
int lhs_length,
const std::uint32_t rhs_words[4],
int rhs_length,
std::uint32_t quotient_words[4],
std::uint32_t remainder_words[4]) {
const int shift = std::countl_zero(rhs_words[rhs_length - 1]);
std::uint32_t normalized_lhs[5] = {};
std::uint32_t normalized_rhs[5] = {};
shift_left_words(lhs_words, lhs_length, shift, normalized_lhs);
shift_left_words(rhs_words, rhs_length, shift, normalized_rhs);
const int quotient_length = lhs_length - rhs_length + 1;
for (int j = quotient_length - 1; j >= 0; --j) {
const std::uint64_t numerator =
(static_cast<std::uint64_t>(normalized_lhs[j + rhs_length])
<< 32) |
normalized_lhs[j + rhs_length - 1];
std::uint64_t estimate = numerator / normalized_rhs[rhs_length - 1];
std::uint64_t estimate_remainder =
numerator % normalized_rhs[rhs_length - 1];
while (estimate == word_base() ||
estimate * normalized_rhs[rhs_length - 2] >
word_base() * estimate_remainder +
normalized_lhs[j + rhs_length - 2]) {
--estimate;
estimate_remainder += normalized_rhs[rhs_length - 1];
if (estimate_remainder >= word_base()) {
break;
}
}
quotient_words[j] = static_cast<std::uint32_t>(estimate);
if (subtract_product(normalized_lhs, j, normalized_rhs, rhs_length,
quotient_words[j])) {
--quotient_words[j];
add_back(normalized_lhs, j, normalized_rhs, rhs_length);
}
}
std::uint32_t shifted_remainder[4] = {};
shift_right_words(normalized_lhs, rhs_length, shift, shifted_remainder);
for (int i = 0; i < rhs_length; ++i) {
remainder_words[i] = shifted_remainder[i];
}
}
// 呼び出し元で lhs >= rhs > 0 が保証されていることを仮定する。
static constexpr void div_mod_32bit_words(UInt128 lhs, UInt128 rhs,
UInt128 "ient,
UInt128 &remainder) {
std::uint32_t lhs_words[4] = {};
std::uint32_t rhs_words[4] = {};
to_words(lhs, lhs_words);
to_words(rhs, rhs_words);
const int lhs_length = word_length(lhs_words);
const int rhs_length = word_length(rhs_words);
if (rhs_length == 1) {
std::uint32_t rem = 0;
quotient = div_mod_uint32(lhs, rhs_words[0], rem);
remainder = UInt128(rem);
return;
}
std::uint32_t quotient_words[4] = {};
std::uint32_t remainder_words[4] = {};
div_mod_knuth(lhs_words, lhs_length, rhs_words, rhs_length,
quotient_words, remainder_words);
quotient = from_32bit_words(quotient_words);
remainder = from_32bit_words(remainder_words);
}
// rhs ≠ 0 が呼び出し元で保証されていることを仮定する。
static constexpr void div_mod_unchecked(UInt128 lhs, UInt128 rhs,
UInt128 "ient,
UInt128 &remainder) {
UInt128 q, r;
if (rhs.high_ == 0) {
if (lhs.high_ == 0) {
q = UInt128(lhs.low_ / rhs.low_);
r = UInt128(lhs.low_ % rhs.low_);
} else if (rhs.low_ == 1) {
q = lhs;
r = UInt128{};
} else if (rhs.low_ <= std::numeric_limits<std::uint32_t>::max()) {
std::uint32_t rem = 0;
q = div_mod_uint32(lhs, static_cast<std::uint32_t>(rhs.low_),
rem);
r = UInt128(rem);
} else {
div_mod_32bit_words(lhs, rhs, q, r);
}
} else if (lhs < rhs) {
q = UInt128{};
r = lhs;
} else {
div_mod_32bit_words(lhs, rhs, q, r);
}
quotient = q;
remainder = r;
}
static constexpr UInt128 multiply_u64(std::uint64_t lhs,
std::uint64_t rhs) {
constexpr std::uint64_t mask = (std::uint64_t{1} << 32) - 1;
const std::uint64_t lhs_low = lhs & mask;
const std::uint64_t lhs_high = lhs >> 32;
const std::uint64_t rhs_low = rhs & mask;
const std::uint64_t rhs_high = rhs >> 32;
const std::uint64_t p00 = lhs_low * rhs_low;
const std::uint64_t p01 = lhs_low * rhs_high;
const std::uint64_t p10 = lhs_high * rhs_low;
const std::uint64_t p11 = lhs_high * rhs_high;
const std::uint64_t middle = (p00 >> 32) + (p01 & mask) + (p10 & mask);
return from_words(p11 + (p01 >> 32) + (p10 >> 32) + (middle >> 32),
(middle << 32) | (p00 & mask));
}
std::uint64_t high_ = 0;
std::uint64_t low_ = 0;
};
class Int128 {
public:
constexpr Int128() = default;
template <class T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
constexpr Int128(T value) : value_(value) {}
static constexpr Int128 from_words(std::uint64_t high, std::uint64_t low) {
Int128 value;
value.value_ = UInt128::from_words(high, low);
return value;
}
constexpr std::uint64_t high() const { return value_.high(); }
constexpr std::uint64_t low() const { return value_.low(); }
template <class T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
explicit constexpr operator T() const {
return static_cast<T>(value_.low());
}
explicit constexpr operator bool() const { return value_ != UInt128{}; }
constexpr bool is_negative() const { return (value_.high() >> 63) != 0; }
explicit constexpr operator long double() const {
const bool negative = is_negative();
const long double magnitude =
static_cast<long double>(negative ? -value_ : value_);
return negative ? -magnitude : magnitude;
}
friend constexpr bool operator==(Int128 lhs, Int128 rhs) {
return lhs.value_ == rhs.value_;
}
friend constexpr bool operator!=(Int128 lhs, Int128 rhs) {
return !(lhs == rhs);
}
friend constexpr bool operator<(Int128 lhs, Int128 rhs) {
const bool lhs_negative = lhs.is_negative();
const bool rhs_negative = rhs.is_negative();
if (lhs_negative != rhs_negative) {
return lhs_negative;
}
return lhs.value_ < rhs.value_;
}
friend constexpr bool operator>(Int128 lhs, Int128 rhs) {
return rhs < lhs;
}
friend constexpr bool operator<=(Int128 lhs, Int128 rhs) {
return !(rhs < lhs);
}
friend constexpr bool operator>=(Int128 lhs, Int128 rhs) {
return !(lhs < rhs);
}
constexpr Int128 operator+() const { return *this; }
constexpr Int128 operator-() const {
Int128 value;
value.value_ = -value_;
return value;
}
constexpr Int128 &operator+=(Int128 rhs) {
value_ += rhs.value_;
return *this;
}
constexpr Int128 &operator-=(Int128 rhs) {
value_ -= rhs.value_;
return *this;
}
constexpr Int128 &operator*=(Int128 rhs) {
value_ *= rhs.value_;
return *this;
}
constexpr Int128 &operator/=(Int128 rhs) {
*this = *this / rhs;
return *this;
}
constexpr Int128 &operator%=(Int128 rhs) {
*this = *this % rhs;
return *this;
}
friend constexpr Int128 operator+(Int128 lhs, Int128 rhs) {
lhs += rhs;
return lhs;
}
friend constexpr Int128 operator-(Int128 lhs, Int128 rhs) {
lhs -= rhs;
return lhs;
}
friend constexpr Int128 operator*(Int128 lhs, Int128 rhs) {
lhs *= rhs;
return lhs;
}
static constexpr void div_mod(Int128 lhs, Int128 rhs, Int128 "ient,
Int128 &remainder) {
assert(rhs != Int128{});
assert(!(lhs == min_value() && rhs == Int128(-1)));
const bool lhs_negative = lhs.is_negative();
const bool rhs_negative = rhs.is_negative();
const bool quotient_negative = lhs_negative != rhs_negative;
UInt128 quotient_abs;
UInt128 remainder_abs;
UInt128::div_mod_unchecked(abs_unsigned(lhs), abs_unsigned(rhs),
quotient_abs, remainder_abs);
const Int128 q =
from_unsigned_unchecked(quotient_abs, quotient_negative);
const Int128 r = from_unsigned_unchecked(remainder_abs, lhs_negative);
quotient = q;
remainder = r;
}
friend constexpr Int128 operator/(Int128 lhs, Int128 rhs) {
Int128 quotient;
Int128 remainder;
div_mod(lhs, rhs, quotient, remainder);
return quotient;
}
friend constexpr Int128 operator%(Int128 lhs, Int128 rhs) {
Int128 quotient;
Int128 remainder;
div_mod(lhs, rhs, quotient, remainder);
return remainder;
}
private:
static constexpr Int128 from_twos_complement(UInt128 value) {
Int128 result;
result.value_ = value;
return result;
}
static constexpr Int128 min_value() {
return from_words(std::uint64_t{1} << 63, 0);
}
static constexpr UInt128 abs_unsigned(Int128 value) {
if (!value.is_negative()) {
return value.value_;
}
return -value.value_;
}
static constexpr Int128 from_unsigned(UInt128 value, bool negative) {
if (!negative) {
assert((value.high() >> 63) == 0);
} else {
assert(value <= UInt128::from_words(std::uint64_t{1} << 63, 0));
}
return from_unsigned_unchecked(value, negative);
}
// 引数 negative で指定した通りに符号を付けた value が
// Int128 の範囲に収まることを仮定する。
static constexpr Int128 from_unsigned_unchecked(UInt128 value,
bool negative) {
return from_twos_complement(negative ? -value : value);
}
friend std::istream &operator>>(std::istream &input, Int128 &value);
friend std::ostream &operator<<(std::ostream &output, Int128 value);
UInt128 value_;
};
namespace int128_internal {
// first < text.size() が呼び出し元で保証されていることを仮定する。
inline UInt128 read_uint128_decimal(const std::string &text,
std::size_t first) {
UInt128 value = 0;
std::size_t i = first;
const std::size_t first_len = (text.size() - first) % 9;
if (first_len != 0) {
std::uint32_t chunk = 0;
for (std::size_t j = 0; j < first_len; ++j) {
assert('0' <= text[i] && text[i] <= '9');
chunk = chunk * 10 + static_cast<std::uint32_t>(text[i] - '0');
++i;
}
value = UInt128(chunk);
}
while (i < text.size()) {
std::uint32_t chunk = 0;
for (int j = 0; j < 9; ++j) {
assert('0' <= text[i] && text[i] <= '9');
chunk = chunk * 10 + static_cast<std::uint32_t>(text[i] - '0');
++i;
}
value *= UInt128(1000000000);
value += UInt128(chunk);
}
return value;
}
inline void write_padded_9(std::ostream &output, std::uint32_t value) {
char buffer[9];
for (int i = 8; i >= 0; --i) {
buffer[i] = static_cast<char>('0' + value % 10);
value /= 10;
}
output.write(buffer, 9);
}
inline void write_uint128_decimal(std::ostream &output, UInt128 value) {
if (value == UInt128(0)) {
output << '0';
return;
}
constexpr std::uint32_t base = 1000000000;
std::uint32_t parts[5] = {};
int size = 0;
while (value != UInt128(0)) {
UInt128 quotient;
UInt128 remainder;
UInt128::div_mod(value, UInt128(base), quotient, remainder);
parts[size] = static_cast<std::uint32_t>(remainder);
value = quotient;
++size;
}
output << parts[size - 1];
for (int i = size - 2; i >= 0; --i) {
write_padded_9(output, parts[i]);
}
}
} // namespace int128_internal
inline std::istream &operator>>(std::istream &input, UInt128 &value) {
std::string text;
input >> text;
if (!input) {
return input;
}
std::size_t first = 0;
if (text[first] == '+') {
++first;
}
assert(first < text.size());
assert(text[first] != '-');
value = int128_internal::read_uint128_decimal(text, first);
return input;
}
inline std::ostream &operator<<(std::ostream &output, UInt128 value) {
int128_internal::write_uint128_decimal(output, value);
return output;
}
inline std::istream &operator>>(std::istream &input, Int128 &value) {
std::string text;
input >> text;
if (!input) {
return input;
}
std::size_t first = 0;
bool negative = false;
if (text[first] == '+' || text[first] == '-') {
negative = text[first] == '-';
++first;
}
assert(first < text.size());
const UInt128 magnitude =
int128_internal::read_uint128_decimal(text, first);
value = Int128::from_unsigned(magnitude, negative);
return input;
}
inline std::ostream &operator<<(std::ostream &output, Int128 value) {
if (value.is_negative()) {
output << '-';
int128_internal::write_uint128_decimal(output,
Int128::abs_unsigned(value));
} else {
int128_internal::write_uint128_decimal(
output, UInt128::from_words(value.high(), value.low()));
}
return output;
}
} // namespace NicheLibrary
namespace std {
template <> class numeric_limits<NicheLibrary::UInt128> {
public:
static constexpr bool is_specialized = true;
static constexpr NicheLibrary::UInt128 min() noexcept { return 0; }
static constexpr NicheLibrary::UInt128 max() noexcept {
return NicheLibrary::UInt128::from_words(~std::uint64_t{},
~std::uint64_t{});
}
static constexpr NicheLibrary::UInt128 lowest() noexcept { return 0; }
static constexpr int digits = 128;
static constexpr int digits10 = 38;
static constexpr int max_digits10 = 0;
static constexpr bool is_signed = false;
static constexpr bool is_integer = true;
static constexpr bool is_exact = true;
static constexpr int radix = 2;
static constexpr NicheLibrary::UInt128 epsilon() noexcept { return 0; }
static constexpr NicheLibrary::UInt128 round_error() noexcept { return 0; }
static constexpr int min_exponent = 0;
static constexpr int min_exponent10 = 0;
static constexpr int max_exponent = 0;
static constexpr int max_exponent10 = 0;
static constexpr bool has_infinity = false;
static constexpr bool has_quiet_NaN = false;
static constexpr bool has_signaling_NaN = false;
static constexpr float_denorm_style has_denorm = denorm_absent;
static constexpr bool has_denorm_loss = false;
static constexpr NicheLibrary::UInt128 infinity() noexcept { return 0; }
static constexpr NicheLibrary::UInt128 quiet_NaN() noexcept { return 0; }
static constexpr NicheLibrary::UInt128 signaling_NaN() noexcept {
return 0;
}
static constexpr NicheLibrary::UInt128 denorm_min() noexcept { return 0; }
static constexpr bool is_iec559 = false;
static constexpr bool is_bounded = true;
static constexpr bool is_modulo = true;
static constexpr bool traps = false;
static constexpr bool tinyness_before = false;
static constexpr float_round_style round_style = round_toward_zero;
};
template <> class numeric_limits<NicheLibrary::Int128> {
public:
static constexpr bool is_specialized = true;
static constexpr NicheLibrary::Int128 min() noexcept {
return NicheLibrary::Int128::from_words(std::uint64_t{1} << 63, 0);
}
static constexpr NicheLibrary::Int128 max() noexcept {
return NicheLibrary::Int128::from_words((std::uint64_t{1} << 63) - 1,
~std::uint64_t{});
}
static constexpr NicheLibrary::Int128 lowest() noexcept { return min(); }
static constexpr int digits = 127;
static constexpr int digits10 = 38;
static constexpr int max_digits10 = 0;
static constexpr bool is_signed = true;
static constexpr bool is_integer = true;
static constexpr bool is_exact = true;
static constexpr int radix = 2;
static constexpr NicheLibrary::Int128 epsilon() noexcept { return 0; }
static constexpr NicheLibrary::Int128 round_error() noexcept { return 0; }
static constexpr int min_exponent = 0;
static constexpr int min_exponent10 = 0;
static constexpr int max_exponent = 0;
static constexpr int max_exponent10 = 0;
static constexpr bool has_infinity = false;
static constexpr bool has_quiet_NaN = false;
static constexpr bool has_signaling_NaN = false;
static constexpr float_denorm_style has_denorm = denorm_absent;
static constexpr bool has_denorm_loss = false;
static constexpr NicheLibrary::Int128 infinity() noexcept { return 0; }
static constexpr NicheLibrary::Int128 quiet_NaN() noexcept { return 0; }
static constexpr NicheLibrary::Int128 signaling_NaN() noexcept { return 0; }
static constexpr NicheLibrary::Int128 denorm_min() noexcept { return 0; }
static constexpr bool is_iec559 = false;
static constexpr bool is_bounded = true;
static constexpr bool is_modulo = false;
static constexpr bool traps = false;
static constexpr bool tinyness_before = false;
static constexpr float_round_style round_style = round_toward_zero;
};
} // namespace std
#line 25 "geometry/line-convex-polygon-intersection.hpp"
namespace line_convex_polygon_intersection_internal {
template <class T>
constexpr bool is_integer_v =
std::numeric_limits<std::remove_cv_t<T>>::is_integer;
template <class T>
constexpr bool is_signed_v =
std::numeric_limits<std::remove_cv_t<T>>::is_signed;
template <class T>
constexpr bool use_int128_by_default_v =
std::is_integral_v<std::remove_cv_t<T>> &&
!std::is_same_v<std::remove_cv_t<T>, bool> &&
std::numeric_limits<std::remove_cv_t<T>>::digits <= 64;
template <class T>
using enable_if_integer_t = std::enable_if_t<is_integer_v<T>, int>;
template <class Point, class = void> struct has_member_xy : std::false_type {};
template <class Point>
struct has_member_xy<Point,
std::void_t<decltype(std::declval<const Point &>().x),
decltype(std::declval<const Point &>().y)>>
: std::true_type {};
template <class Point> decltype(auto) get_x(const Point &point) {
if constexpr (has_member_xy<Point>::value) {
return point.x;
} else {
using std::real;
return real(point);
}
}
template <class Point> decltype(auto) get_y(const Point &point) {
if constexpr (has_member_xy<Point>::value) {
return point.y;
} else {
using std::imag;
return imag(point);
}
}
template <class Point>
using coordinate_t =
std::remove_cvref_t<decltype(get_x(std::declval<const Point &>()))>;
template <class Coord, bool = is_integer_v<Coord>> struct default_calc_type {
using type = long double;
};
template <class Coord> struct default_calc_type<Coord, true> {
using type = std::conditional_t<use_int128_by_default_v<Coord>,
NicheLibrary::Int128, Coord>;
};
template <class Point>
using default_calc_t = typename default_calc_type<coordinate_t<Point>>::type;
template <class Point, class Calc> struct resolved_calc {
using type = Calc;
};
template <class Point> struct resolved_calc<Point, void> {
using type = default_calc_t<Point>;
};
template <class Point, class Calc>
using resolved_calc_t = typename resolved_calc<Point, Calc>::type;
template <class T> int sign_value(const T &value) {
if constexpr (is_integer_v<T>) {
if (value < 0) {
return -1;
}
if (value > 0) {
return 1;
}
return 0;
} else {
constexpr long double eps = 1e-12L;
if (value < -eps) {
return -1;
}
if (value > eps) {
return 1;
}
return 0;
}
}
template <class Point, class Calc> Calc to_calc_x(const Point &point) {
return static_cast<Calc>(get_x(point));
}
template <class Point, class Calc> Calc to_calc_y(const Point &point) {
return static_cast<Calc>(get_y(point));
}
template <class Calc> struct LineParameters {
Calc base_x;
Calc base_y;
Calc direction_x;
Calc direction_y;
};
template <class Point, class Calc>
LineParameters<Calc> make_line_parameters(const Point &line_a,
const Point &line_b) {
const Calc base_x = to_calc_x<Point, Calc>(line_a);
const Calc base_y = to_calc_y<Point, Calc>(line_a);
return {base_x, base_y, to_calc_x<Point, Calc>(line_b) - base_x,
to_calc_y<Point, Calc>(line_b) - base_y};
}
template <class Point, class Calc>
Calc line_height(const Point &point, const LineParameters<Calc> &line) {
return (to_calc_x<Point, Calc>(point) - line.base_x) * line.direction_y -
(to_calc_y<Point, Calc>(point) - line.base_y) * line.direction_x;
}
template <class Point, class Calc>
Point make_point(const Calc &x, const Calc &y) {
using Coord = coordinate_t<Point>;
return Point(static_cast<Coord>(x), static_cast<Coord>(y));
}
template <class T> bool equals_value(const T &lhs, const T &rhs) {
return sign_value(lhs - rhs) == 0;
}
template <class T, enable_if_integer_t<T> = 0> T abs_value(T value) {
if constexpr (is_signed_v<T>) {
return value < 0 ? -value : value;
} else {
return value;
}
}
template <class T, enable_if_integer_t<T> = 0> T gcd_value(T a, T b) {
a = abs_value(a);
b = abs_value(b);
while (b != 0) {
T t = a % b;
a = b;
b = t;
}
return a;
}
template <class T, enable_if_integer_t<T> = 0> T gcd3_value(T a, T b, T c) {
return gcd_value(gcd_value(a, b), c);
}
template <class Func> int first_non_negative(int length, Func &&func) {
int left = 0;
int right = length;
while (left + 1 < right) {
int mid = (left + right) >> 1;
if (sign_value(func(mid)) < 0) {
left = mid;
} else {
right = mid;
}
}
return right;
}
template <class Func> int last_non_positive(int length, Func &&func) {
int left = 0;
int right = length;
while (left < right) {
int mid = (left + right + 1) >> 1;
if (sign_value(func(mid)) <= 0) {
left = mid;
} else {
right = mid - 1;
}
}
return left;
}
template <class Func> int last_non_negative(int length, Func &&func) {
int left = 0;
int right = length;
while (left < right) {
int mid = (left + right + 1) >> 1;
if (sign_value(func(mid)) >= 0) {
left = mid;
} else {
right = mid - 1;
}
}
return left;
}
inline int wrap_nearby_index(int index, int size) {
if (index < 0) {
return index + size;
}
if (index >= size) {
return index - size;
}
return index;
}
inline int distance_forward(int from, int to, int n) {
if (from <= to) {
return to - from;
}
return to + n - from;
}
template <class T> int compare_value(const T &lhs, const T &rhs) {
return sign_value(lhs - rhs);
}
template <class Func>
std::pair<int, int> find_extreme_vertices(int n, Func &&height) {
using Value = std::remove_cvref_t<decltype(height(0))>;
int minimum_index = -1;
int maximum_index = -1;
const Value first = height(0);
const Value second = height(1);
const Value last = height(n - 1);
const int last_vs_first = compare_value(last, first);
const int first_vs_second = compare_value(first, second);
if (last_vs_first == 0) {
if (first_vs_second < 0) {
minimum_index = n - 1;
} else {
maximum_index = n - 1;
}
} else if (last_vs_first > 0 && first_vs_second <= 0) {
minimum_index = 0;
} else if (last_vs_first < 0 && first_vs_second >= 0) {
maximum_index = 0;
} else {
const Value denominator = static_cast<Value>(n - 1);
auto lower_envelope_scaled = [&](int index) -> Value {
const Value left_weight = static_cast<Value>(n - 1 - index);
const Value right_weight = static_cast<Value>(index);
const Value chord_scaled =
first * left_weight + last * right_weight;
const Value height_scaled = height(index) * denominator;
if (compare_value(chord_scaled, height_scaled) < 0) {
return chord_scaled;
}
return height_scaled;
};
minimum_index = first_non_negative(n - 1, [&](int index) -> Value {
return lower_envelope_scaled(index + 1) -
lower_envelope_scaled(index);
});
}
if (maximum_index == -1) {
const int base = minimum_index - n;
const int offset = first_non_negative(n, [&](int index) -> Value {
return height(base + index) - height(base + index + 1);
});
maximum_index = wrap_nearby_index(minimum_index + offset, n);
} else {
const int base = maximum_index - n;
const int offset = first_non_negative(n, [&](int index) -> Value {
return height(base + index + 1) - height(base + index);
});
minimum_index = wrap_nearby_index(maximum_index + offset, n);
}
return {minimum_index, maximum_index};
}
template <class Real, class T> Real number_as_real(const T &value) {
return static_cast<Real>(value);
}
template <class Real> Real number_as_real(NicheLibrary::UInt128 value) {
return std::ldexp(static_cast<Real>(value.high()), 64) +
static_cast<Real>(value.low());
}
template <class Real> Real number_as_real(NicheLibrary::Int128 value) {
const NicheLibrary::UInt128 bits =
NicheLibrary::UInt128::from_words(value.high(), value.low());
if (value.is_negative()) {
return -number_as_real<Real>(-bits);
}
return number_as_real<Real>(bits);
}
} // namespace line_convex_polygon_intersection_internal
template <class T> struct LineConvexHullIntersectionPoint {
T x_numerator{};
T y_numerator{};
T denominator = T(1);
template <class Real> Real x_as() const {
assert(denominator > 0);
namespace lpi_internal = line_convex_polygon_intersection_internal;
return lpi_internal::number_as_real<Real>(x_numerator) /
lpi_internal::number_as_real<Real>(denominator);
}
template <class Real> Real y_as() const {
assert(denominator > 0);
namespace lpi_internal = line_convex_polygon_intersection_internal;
return lpi_internal::number_as_real<Real>(y_numerator) /
lpi_internal::number_as_real<Real>(denominator);
}
template <class Point> Point to_point() const {
namespace lpi_internal = line_convex_polygon_intersection_internal;
using Coord = lpi_internal::coordinate_t<Point>;
if constexpr (lpi_internal::is_integer_v<Coord>) {
assert(denominator == T(1));
return Point{static_cast<Coord>(x_numerator),
static_cast<Coord>(y_numerator)};
} else {
assert(denominator > 0);
const long double real_denominator =
lpi_internal::number_as_real<long double>(denominator);
return Point{
static_cast<Coord>(
lpi_internal::number_as_real<long double>(x_numerator) /
real_denominator),
static_cast<Coord>(
lpi_internal::number_as_real<long double>(y_numerator) /
real_denominator)};
}
}
friend constexpr bool
operator==(const LineConvexHullIntersectionPoint &,
const LineConvexHullIntersectionPoint &) = default;
};
template <class T>
using LinePolygonIntersectionPoint = LineConvexHullIntersectionPoint<T>;
namespace line_convex_polygon_intersection_internal {
template <class Point, class Calc, bool = is_integer_v<coordinate_t<Point>>>
struct result_value_type {
using type = Point;
};
template <class Point, class Calc> struct result_value_type<Point, Calc, true> {
using type = LineConvexHullIntersectionPoint<Calc>;
};
template <class Point, class Calc>
using result_value_t = typename result_value_type<Point, Calc>::type;
} // namespace line_convex_polygon_intersection_internal
template <class Point, class Calc = void>
using LineConvexHullIntersectionValue =
line_convex_polygon_intersection_internal::result_value_t<
Point, line_convex_polygon_intersection_internal::resolved_calc_t<
Point, Calc>>;
template <class Point, class Calc = void>
using LineConvexHullIntersectionResult =
std::vector<LineConvexHullIntersectionValue<Point, Calc>>;
template <class Point, class Calc = void>
using LinePolygonIntersectionValue =
LineConvexHullIntersectionValue<Point, Calc>;
template <class Point, class Calc = void>
using LinePolygonIntersectionResult =
LineConvexHullIntersectionResult<Point, Calc>;
namespace line_convex_polygon_intersection_internal {
template <class Point, class Calc>
LineConvexHullIntersectionPoint<Calc> make_integral_point(const Point &point) {
return {to_calc_x<Point, Calc>(point), to_calc_y<Point, Calc>(point),
Calc(1)};
}
template <class Point, class Calc>
LineConvexHullIntersectionPoint<Calc>
line_edge_intersection_integral(const LineParameters<Calc> &line,
const Point &segment_a,
const Point &segment_b) {
const Calc segment_a_x = to_calc_x<Point, Calc>(segment_a);
const Calc segment_a_y = to_calc_y<Point, Calc>(segment_a);
const Calc edge_x = to_calc_x<Point, Calc>(segment_b) - segment_a_x;
const Calc edge_y = to_calc_y<Point, Calc>(segment_b) - segment_a_y;
Calc denominator = edge_x * line.direction_y - edge_y * line.direction_x;
Calc numerator = (line.base_x - segment_a_x) * line.direction_y -
(line.base_y - segment_a_y) * line.direction_x;
Calc x_numerator = segment_a_x * denominator + edge_x * numerator;
Calc y_numerator = segment_a_y * denominator + edge_y * numerator;
if (denominator < 0) {
denominator = -denominator;
x_numerator = -x_numerator;
y_numerator = -y_numerator;
}
const Calc g = gcd3_value(x_numerator, y_numerator, denominator);
x_numerator /= g;
y_numerator /= g;
denominator /= g;
return {x_numerator, y_numerator, denominator};
}
template <class Point, class Calc>
Point line_edge_intersection_floating(const LineParameters<Calc> &line,
const Point &segment_a,
const Point &segment_b) {
const Calc segment_a_x = to_calc_x<Point, Calc>(segment_a);
const Calc segment_a_y = to_calc_y<Point, Calc>(segment_a);
const Calc edge_x = to_calc_x<Point, Calc>(segment_b) - segment_a_x;
const Calc edge_y = to_calc_y<Point, Calc>(segment_b) - segment_a_y;
const Calc denominator =
edge_x * line.direction_y - edge_y * line.direction_x;
const Calc numerator = (line.base_x - segment_a_x) * line.direction_y -
(line.base_y - segment_a_y) * line.direction_x;
const Calc t = numerator / denominator;
const Calc x = segment_a_x + edge_x * t;
const Calc y = segment_a_y + edge_y * t;
return make_point<Point, Calc>(x, y);
}
template <class Point, class Calc>
result_value_t<Point, Calc> make_vertex_result(const Point &point) {
if constexpr (is_integer_v<coordinate_t<Point>>) {
return make_integral_point<Point, Calc>(point);
} else {
return point;
}
}
template <class Point, class Calc>
result_value_t<Point, Calc> make_edge_result(const LineParameters<Calc> &line,
const Point &segment_a,
const Point &segment_b) {
if constexpr (is_integer_v<coordinate_t<Point>>) {
return line_edge_intersection_integral<Point, Calc>(line, segment_a,
segment_b);
} else {
return line_edge_intersection_floating<Point, Calc>(line, segment_a,
segment_b);
}
}
template <class Point, class Calc>
bool equals_point(const Point &lhs, const Point &rhs) {
return equals_value(to_calc_x<Point, Calc>(lhs),
to_calc_x<Point, Calc>(rhs)) &&
equals_value(to_calc_y<Point, Calc>(lhs),
to_calc_y<Point, Calc>(rhs));
}
template <class Point, class Calc>
std::vector<result_value_t<Point, Calc>>
line_degenerate_convex_hull_intersection(const std::vector<Point> &hull,
const LineParameters<Calc> &line) {
std::vector<result_value_t<Point, Calc>> result;
if (hull.empty()) {
return result;
}
if (hull.size() == 1) {
if (sign_value(line_height<Point, Calc>(hull[0], line)) == 0) {
result.push_back(make_vertex_result<Point, Calc>(hull[0]));
}
return result;
}
assert((!equals_point<Point, Calc>(hull[0], hull[1])));
const Calc first_height = line_height<Point, Calc>(hull[0], line);
const Calc second_height = line_height<Point, Calc>(hull[1], line);
const int first_sign = sign_value(first_height);
const int second_sign = sign_value(second_height);
if (first_sign == 0 && second_sign == 0) {
result.reserve(2);
result.push_back(make_vertex_result<Point, Calc>(hull[0]));
result.push_back(make_vertex_result<Point, Calc>(hull[1]));
return result;
}
if (first_sign == 0) {
result.push_back(make_vertex_result<Point, Calc>(hull[0]));
return result;
}
if (second_sign == 0) {
result.push_back(make_vertex_result<Point, Calc>(hull[1]));
return result;
}
if ((first_sign < 0 && second_sign > 0) ||
(first_sign > 0 && second_sign < 0)) {
result.push_back(make_edge_result<Point, Calc>(line, hull[0], hull[1]));
}
return result;
}
template <class Point, class Calc>
std::vector<result_value_t<Point, Calc>>
line_strict_convex_polygon_intersection(const std::vector<Point> &polygon,
const LineParameters<Calc> &line) {
const int n = static_cast<int>(polygon.size());
auto vertex = [&](int index) -> const Point & {
return polygon[wrap_nearby_index(index, n)];
};
auto height = [&](int index) -> Calc {
return line_height<Point, Calc>(vertex(index), line);
};
auto chain_index = [&](int start, int step, int offset) -> int {
return line_convex_polygon_intersection_internal::wrap_nearby_index(
start + step * offset, n);
};
const auto [minimum_index, maximum_index] =
find_extreme_vertices(n, height);
const Calc minimum_value = height(minimum_index);
const Calc maximum_value = height(maximum_index);
const int minimum_sign = sign_value(minimum_value);
const int maximum_sign = sign_value(maximum_value);
LinePolygonIntersectionResult<Point, Calc> result;
if (minimum_sign > 0 || maximum_sign < 0) {
return result;
}
result.reserve(2);
const int forward_length =
distance_forward(minimum_index, maximum_index, n);
const int backward_length = n - forward_length;
if (minimum_sign == 0) {
const int forward_zero =
last_non_positive(forward_length, [&](int offset) {
return height(minimum_index + offset);
});
const int backward_zero =
last_non_positive(backward_length, [&](int offset) {
return height(minimum_index - offset);
});
const int forward_index = chain_index(minimum_index, +1, forward_zero);
const int backward_index =
chain_index(minimum_index, -1, backward_zero);
result.push_back(
make_vertex_result<Point, Calc>(vertex(forward_index)));
if (forward_index != backward_index) {
result.push_back(
make_vertex_result<Point, Calc>(vertex(backward_index)));
}
return result;
}
if (maximum_sign == 0) {
const int forward_length_from_max =
distance_forward(maximum_index, minimum_index, n);
const int backward_length_from_max = n - forward_length_from_max;
const int forward_zero =
last_non_negative(forward_length_from_max, [&](int offset) {
return height(maximum_index + offset);
});
const int backward_zero =
last_non_negative(backward_length_from_max, [&](int offset) {
return height(maximum_index - offset);
});
const int forward_index = chain_index(maximum_index, +1, forward_zero);
const int backward_index =
chain_index(maximum_index, -1, backward_zero);
result.push_back(
make_vertex_result<Point, Calc>(vertex(forward_index)));
if (forward_index != backward_index) {
result.push_back(
make_vertex_result<Point, Calc>(vertex(backward_index)));
}
return result;
}
const int first_cross = first_non_negative(forward_length, [&](int offset) {
return height(minimum_index + offset);
});
const int first_cross_index = chain_index(minimum_index, +1, first_cross);
const Calc first_cross_height = height(first_cross_index);
if (sign_value(first_cross_height) == 0) {
result.push_back(
make_vertex_result<Point, Calc>(vertex(first_cross_index)));
} else {
const int prev_index = chain_index(minimum_index, +1, first_cross - 1);
result.push_back(make_edge_result<Point, Calc>(
line, vertex(prev_index), vertex(first_cross_index)));
}
const int second_cross =
first_non_negative(backward_length, [&](int offset) {
return height(minimum_index - offset);
});
const int second_cross_index = chain_index(minimum_index, -1, second_cross);
const Calc second_cross_height = height(second_cross_index);
if (sign_value(second_cross_height) == 0) {
result.push_back(
make_vertex_result<Point, Calc>(vertex(second_cross_index)));
} else {
const int prev_index = chain_index(minimum_index, -1, second_cross - 1);
result.push_back(make_edge_result<Point, Calc>(
line, vertex(prev_index), vertex(second_cross_index)));
}
return result;
}
} // namespace line_convex_polygon_intersection_internal
template <class Point, class Calc = void>
LineConvexHullIntersectionResult<Point, Calc>
line_convex_hull_intersection(const std::vector<Point> &hull,
const Point &line_a, const Point &line_b) {
namespace lpi_internal = line_convex_polygon_intersection_internal;
using Coord = lpi_internal::coordinate_t<Point>;
using Number = lpi_internal::resolved_calc_t<Point, Calc>;
static_assert(!lpi_internal::is_integer_v<Coord> ||
lpi_internal::is_signed_v<Coord>,
"integer coordinate type must be signed");
static_assert(!lpi_internal::is_integer_v<Coord> ||
lpi_internal::is_integer_v<Number>,
"integer coordinate type requires integer calculation type");
static_assert(!lpi_internal::is_integer_v<Number> ||
lpi_internal::is_signed_v<Number>,
"integer calculation type must be signed");
const auto line =
lpi_internal::make_line_parameters<Point, Number>(line_a, line_b);
assert(lpi_internal::sign_value(line.direction_x) != 0 ||
lpi_internal::sign_value(line.direction_y) != 0);
if (hull.size() <= 2) {
return lpi_internal::line_degenerate_convex_hull_intersection<Point,
Number>(
hull, line);
}
return lpi_internal::line_strict_convex_polygon_intersection<Point, Number>(
hull, line);
}
template <class Point, class Calc = void>
LinePolygonIntersectionResult<Point, Calc>
line_polygon_intersection(const std::vector<Point> &polygon,
const Point &line_a, const Point &line_b) {
return line_convex_hull_intersection<Point, Calc>(polygon, line_a, line_b);
}
#line 9 "verify/standalone-line-convex-polygon-intersection.test.cpp"
int main() {
using Point = std::complex<long long>;
using NicheLibrary::Int128;
{
std::vector<Point> hull = {Point(2, 3)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 3), Point(1, 3));
assert(result.size() == 1);
assert(result[0].x_numerator == Int128(2));
assert(result[0].y_numerator == Int128(3));
assert(result[0].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(2, 3)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 4), Point(1, 4));
assert(result.empty());
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 2), Point(1, 2));
assert(result.size() == 1);
assert(result[0].x_numerator == Int128(2));
assert(result[0].y_numerator == Int128(2));
assert(result[0].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 0), Point(1, 1));
assert(result.size() == 2);
assert(result[0].denominator == Int128(1));
assert(result[1].denominator == Int128(1));
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 5), Point(1, 5));
assert(result.empty());
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 0), Point(0, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 2), Point(1, 2));
assert(result.size() == 2);
}
{
std::vector<Point> hull = {Point(0, 0), Point(4, 0), Point(0, 4)};
const auto result =
line_convex_hull_intersection(hull, Point(0, 1), Point(3, 2));
assert(result.size() == 2);
auto is_point = [&](int index, Int128 x_numerator, Int128 y_numerator,
Int128 denominator) {
return result[index].x_numerator == x_numerator &&
result[index].y_numerator == y_numerator &&
result[index].denominator == denominator;
};
const bool expected_order =
is_point(0, Int128(0), Int128(1), Int128(1)) &&
is_point(1, Int128(9), Int128(7), Int128(4));
const bool reverse_order =
is_point(1, Int128(0), Int128(1), Int128(1)) &&
is_point(0, Int128(9), Int128(7), Int128(4));
assert(expected_order || reverse_order);
const int rational_index = result[0].denominator == Int128(4) ? 0 : 1;
const auto real_point =
result[rational_index]
.template to_point<std::complex<long double>>();
assert(std::abs(real_point.real() - 2.25L) < 1e-15L);
assert(std::abs(real_point.imag() - 1.75L) < 1e-15L);
}
return 0;
}