Convex Polygons
(geometry/convex_polygon.hpp)

View this file on GitHub
Last update: 2026-10-06 02:48:54+09:00
Include: #include "geometry/convex_polygon.hpp"

Overview

convex_polygon.hpp provides algorithms specialized for convex polygons and a ConvexPolygon<T> query object. The query object normalizes an ordered convex boundary once, then supports point containment, directional extrema, tangents, and chain-area queries efficiently. It also supports Minkowski addition with operator+ and scalar multiplication about the origin with operator*.

The free functions cover convexity testing, normalization, triangulation, diameter, half-plane cuts, intersection construction, and intersection and distance between two convex polygons. Pair queries have both linear-time vector overloads and sublinear overloads for already-normalized ConvexPolygon<T> objects. They can also return a pair of closest points, including points in edge interiors.

Polygons are represented by std::vector<Point<T>>. Their first point is not repeated at the end unless a function explicitly accepts and removes a closing copy.

Basic Functions

template <Coordinate T>
bool is_convex_polygon(
    const std::vector<Point<T>>& polygon,
    bool strict = false,
    long double eps = 1e-12L
);

template <Coordinate T>
std::vector<Point<T>> normalize_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
);

template <Coordinate T>
PointInPolygon point_in_convex_polygon(
    const std::vector<Point<T>>& polygon,
    const Point<T>& point,
    long double eps = 1e-12L
);

Function	Description	Complexity
`is_convex_polygon(polygon, strict, eps)`	Tests the turns of an ordered simple boundary. With `strict == true`, collinear consecutive edges and all-collinear polygons are rejected.	$O(N)$
`normalize_convex_polygon(polygon, eps)`	Removes a closing copy, consecutive duplicates, and redundant collinear vertices; makes the order counterclockwise; and rotates the lowest `(y, x)` vertex to index `0`.	$O(N)$
`point_in_convex_polygon(polygon, point, eps)`	Classifies a point against a strict convex boundary in clockwise or counterclockwise order.	$O(\log N)$

is_convex_polygon assumes the vertices describe a simple polygon boundary. It is not a general self-intersection test. Use is_simple_polygon from polygon.hpp when arbitrary input must be validated. In non-strict mode, an all-collinear boundary with at least three vertices is considered convex.

For predictable logarithmic containment, pass a strict boundary or use ConvexPolygon<T>, which normalizes the input first. Empty polygons, points, and segments are also classified by point_in_convex_polygon.

Query Object

template <Coordinate T>
class ConvexPolygon {
public:
    using Wide = wide_type<T>;

    explicit ConvexPolygon(
        std::vector<Point<T>> polygon,
        long double eps = 1e-12L
    );

    int size() const noexcept;
    bool empty() const noexcept;
    const std::vector<Point<T>>& vertices() const noexcept;
    const Point<T>& operator[](int index) const;
    ConvexPolygon operator+(const ConvexPolygon& other) const;

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    ConvexPolygon<std::common_type_t<T, Scalar>> operator*(Scalar scalar) const;

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    friend ConvexPolygon<std::common_type_t<T, Scalar>> operator*(
        Scalar scalar,
        const ConvexPolygon& polygon
    );

    Wide area2() const;
    Wide chain_area2(int first, int last) const;
    PointInPolygon contains(const Point<T>& point) const;
    std::pair<Wide, int> min_dot(const Point<T>& direction) const;
    std::pair<Wide, int> max_dot(const Point<T>& direction) const;
    std::pair<int, int> tangent_vertices(const Point<T>& point) const;
};

template <Coordinate T>
std::optional<Point<long double>> centroid(
    const ConvexPolygon<T>& polygon,
    long double eps = 1e-12L
);

Operation	Description	Complexity
`ConvexPolygon(polygon, eps)`	Normalizes an ordered convex boundary and builds doubled prefix areas.	$O(N)$ time and memory
`size()`, `empty()`, `vertices()`, `operator[]`	Access the normalized boundary.	$O(1)$
`ConvexPolygon operator+(const ConvexPolygon& other) const`	Returns the Minkowski sum with another polygon of the same coordinate type. Both inputs must be nonempty.	$O(N+M)$ time and memory
`ConvexPolygon<std::common_type_t<T, Scalar>> operator*(Scalar scalar) const`	Returns a polygon with each vertex multiplied by `scalar`.	$O(N)$ time and memory
`friend ConvexPolygon<std::common_type_t<T, Scalar>> operator*(Scalar scalar, const ConvexPolygon& polygon)`	Supports `scalar * polygon` with the same behavior as `polygon * scalar`.	$O(N)$ time and memory
`area2()`	Returns signed twice-area. A nondegenerate normalized polygon has positive area.	$O(1)$
`chain_area2(first, last)`	Returns signed twice-area enclosed by the counterclockwise chain from `first` through `last` and the chord back to `first`.	$O(1)$
`contains(point)`	Classifies a point as `Outside`, `Boundary`, or `Inside`.	$O(\log N)$
`min_dot(direction)`, `max_dot(direction)`	Returns the extreme dot product and one vertex attaining it.	$O(\log N)$
`tangent_vertices(point)`	Returns the two tangent-vertex indices for an external point.	$O(\log N)$
`centroid(polygon, eps)`	Returns the uniformly filled polygon’s centroid, or `nullopt` for an empty or zero-area query object.	$O(N)$

min_dot and max_dot require a nonempty polygon. tangent_vertices requires at least three vertices and a point strictly outside the polygon. Ties may return either endpoint of an extreme edge. No ordering is promised between the two tangent indices.

The arithmetic operators return new normalized query objects and do not mutate their operands. first + second represents all points a + b with a in first and b in second; points and segments are supported. It uses the larger of the two constructor tolerances for normalization and later queries. Coordinate sums and edge differences must fit T.

Multiplication scales about (0, 0) and preserves the constructor tolerance. A negative scalar also rotates the polygon by 180 degrees; the returned boundary remains counterclockwise and starts at its lowest (y, x) vertex. A zero scalar maps a nonempty polygon to the single point (0, 0). Scaling an empty polygon returns an empty polygon. The result coordinate type is std::common_type_t<T, Scalar>, as with Point<T> multiplication, so ConvexPolygon<long long> * 0.5L returns ConvexPolygon<long double>. The common type must satisfy Coordinate, and coordinate products must fit it.

centroid is a free geometry-wide overload rather than a convex-only member. The same name also supports points, segments, triangles, circles, and general simple polygon vectors. The convex overload delegates to the general polygon area-centroid calculation; preprocessing the query object does not make this particular operation constant-time.

Construction and Combination

template <Coordinate T>
std::vector<std::array<Point<T>, 3>> triangulate_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
);

template <Coordinate T>
wide_type<T> convex_diameter2(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
);

template <Coordinate T>
std::vector<Point<long double>> convex_cut(
    const std::vector<Point<T>>& polygon,
    const Line<T>& boundary,
    long double eps = 1e-12L
);

template <Coordinate T>
std::vector<Point<long double>> convex_polygon_intersection(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
bool convex_polygons_intersect(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
long double convex_polygons_distance(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
bool convex_polygons_intersect(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
long double convex_polygons_distance(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
);

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
);

Vector overloads

Function	Description	Complexity
`triangulate_convex_polygon(polygon, eps)`	Returns a counterclockwise fan triangulation after normalization.	$O(N)$
`convex_diameter2(polygon, eps)`	Returns the maximum squared distance between two vertices using rotating calipers.	$O(N)$
`convex_cut(polygon, boundary, eps)`	Intersects the polygon with the closed half-plane to the left of the directed boundary line.	$O(N)$
`convex_polygon_intersection(first, second, eps)`	Constructs the closed intersection, including point or segment degeneracies, by merging the polygons’ angle-sorted half-planes.	$O(N+M)$
`convex_polygons_intersect(first, second, eps)`	Tests whether two closed convex polygons intersect.	$O(N+M)$
`convex_polygons_distance(first, second, eps)`	Returns the minimum Euclidean distance between two closed convex polygons.	$O(N+M)$
`convex_polygons_closest_points(first, second, eps)`	Returns one pair of points attaining the minimum distance.	$O(N+M)$

The $O(N+M)$ bounds in this table apply when first and second are std::vector<Point<T>>. These overloads normalize the boundaries on every call. The intersection and distance queries materialize their Minkowski difference; the closest-points query constructs two temporary query objects. All three use $O(N+M)$ temporary memory.

Preprocessed pair-query overloads

Function	Description	Complexity
`convex_polygons_intersect(first, second, eps)`	Tests whether two closed `ConvexPolygon<T>` objects intersect.	$O(\log(N+M)\log(\min(N,M)+1))$ time, $O(1)$ extra memory
`convex_polygons_distance(first, second, eps)`	Returns the minimum Euclidean distance between two closed `ConvexPolygon<T>` objects.	$O(\log(N+M)\log(\min(N,M)+1))$ time, $O(1)$ extra memory
`convex_polygons_closest_points(first, second, eps)`	Returns one pair of points attaining the minimum distance between two closed `ConvexPolygon<T>` objects.	$O(\log(N+M)\log(\min(N,M)+1))$ time, $O(1)$ extra memory

These bounds apply only when both arguments are ConvexPolygon<T> objects. Constructing those objects still costs $O(N+M)$ total time and memory. The overloads are therefore useful when the same polygons participate in multiple queries, or when the query objects already exist for other operations.

The product of logarithms is intentional: the implementation performs $O(\log(N+M))$ searches on a virtual Minkowski-difference boundary, and one random access into that merged boundary costs $O(\log(\min(N,M)+1))$. It does not build or cache an $N+M$-vertex pairwise boundary. In particular, this complexity must not be read as $O(\log N+\log M)$.

Here $N$ and $M$ are the numbers of vertices after each query object’s normalization. Empty query objects are invalid; points and segments are supported. All pair-query overloads treat the polygons as closed, so sharing a vertex or edge counts as intersection and makes the distance zero. convex_polygons_closest_points returns ClosestPoints{first, second}, where first belongs to the first polygon and second belongs to the second. The returned points use long double because a closest point may lie inside an edge. If the polygons intersect, both returned points describe one common point; when several answers exist, any one may be returned. eps controls geometric classification during the pair query; each object’s constructor tolerance has already been applied during its normalization.

convex_cut and convex_polygon_intersection return Point<long double> because new vertices may be non-integral. A cut boundary is directed and must contain two distinct points. Intersection construction requires two nondegenerate convex inputs; the resulting intersection itself may degenerate to a point or segment.

Minkowski addition is also used internally for the vector pair queries. The three vector pair queries and the three preprocessed pair-query overloads require nonempty inputs. Coordinate negation, addition, and edge differences must fit T. Cross products, dot products, squared distances, and areas must fit wide_type<T>.

Example

#include "geometry/convex_polygon.hpp"

#include <iostream>
#include <vector>

int main() {
    using Point = m1une::geometry::Point<long long>;
    std::vector<Point> vertices;
    vertices.emplace_back(0, 0);
    vertices.emplace_back(4, 0);
    vertices.emplace_back(4, 4);
    vertices.emplace_back(0, 4);

    m1une::geometry::ConvexPolygon<long long> polygon(vertices);
    std::cout << int(polygon.contains(Point(2, 2))) << "\n";  // 2

    auto maximum = polygon.max_dot(Point(1, 0));
    std::cout << static_cast<long long>(maximum.first) << "\n";  // 4

    auto sum = polygon + polygon;
    std::cout << sum.size() << "\n";  // 4, square from (0, 0) to (8, 8)
    auto enlarged = 2 * polygon;
    auto reflected = polygon * -1;
    auto half = polygon * 0.5L;  // ConvexPolygon<long double>
    std::cout << enlarged[2].x << "\n";  // 8
    std::cout << reflected[0].x << "\n";  // -4
    std::cout << half[2].x << "\n";  // 2

    m1une::geometry::Line<long long> boundary{
        Point(2, -1),
        Point(2, 5)
    };
    auto left_half = m1une::geometry::convex_cut(vertices, boundary);
    std::cout << m1une::geometry::polygon_area(left_half) << "\n";  // 8
}

Depends on

Required by

Geometry Bundle (geometry/all.hpp)

Verified with

Code

#ifndef M1UNE_GEOMETRY_CONVEX_POLYGON_HPP
#define M1UNE_GEOMETRY_CONVEX_POLYGON_HPP 1

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <deque>
#include <limits>
#include <numbers>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>

#include "convex_hull.hpp"
#include "half_plane_intersection.hpp"
#include "minkowski_sum.hpp"
#include "polygon.hpp"

namespace m1une {
namespace geometry {

namespace convex_polygon_detail {

inline bool points_close(
    const Point<long double>& first,
    const Point<long double>& second,
    long double eps
) {
    return geometry::distance(first, second) <= eps;
}

inline std::vector<Point<long double>> clean_polygon(
    std::vector<Point<long double>> polygon,
    long double eps
) {
    if (polygon.empty()) return polygon;

    std::vector<Point<long double>> deduplicated;
    for (const Point<long double>& point : polygon) {
        if (
            deduplicated.empty() ||
            !points_close(deduplicated.back(), point, eps)
        ) {
            deduplicated.push_back(point);
        }
    }
    if (
        deduplicated.size() >= 2 &&
        points_close(
            deduplicated.front(),
            deduplicated.back(),
            eps
        )
    ) {
        deduplicated.pop_back();
    }
    if (deduplicated.size() <= 2) return deduplicated;
    std::vector<Point<long double>> cleaned;
    const std::size_t size = deduplicated.size();
    cleaned.reserve(size);
    for (std::size_t index = 0; index < size; ++index) {
        const Point<long double>& previous =
            deduplicated[(index + size - 1) % size];
        const Point<long double>& current = deduplicated[index];
        const Point<long double>& next =
            deduplicated[(index + 1) % size];
        if (
            orientation(previous, current, next, eps) != 0 ||
            dot(current - previous, next - current) < -eps
        ) {
            cleaned.push_back(current);
        }
    }
    return cleaned;
}

}  // namespace convex_polygon_detail

template <Coordinate T>
bool is_convex_polygon(
    const std::vector<Point<T>>& polygon,
    bool strict = false,
    long double eps = 1e-12L
) {
    std::size_t size = polygon.size();
    if (size >= 2 && polygon.front() == polygon.back()) size--;
    if (size < 3) return false;

    int direction = 0;
    for (std::size_t index = 0; index < size; ++index) {
        const Point<T>& current = polygon[index];
        const Point<T>& next = polygon[(index + 1) % size];
        const Point<T>& after = polygon[(index + 2) % size];
        if (current == next) return false;
        const int turn = orientation(current, next, after, eps);
        if (turn == 0) {
            if (strict) return false;
            continue;
        }
        if (direction != 0 && direction != turn) return false;
        direction = turn;
    }
    return !strict || direction != 0;
}

template <Coordinate T>
std::vector<Point<T>> normalize_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    return convex_polygon_detail::normalize_convex_boundary(
        std::move(polygon),
        eps
    );
}

template <Coordinate T>
PointInPolygon point_in_convex_polygon(
    const std::vector<Point<T>>& polygon,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    const std::size_t size = polygon.size();
    if (size == 0) return PointInPolygon::Outside;
    if (size == 1) {
        return distance(polygon[0], point) <= eps
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }
    if (size == 2) {
        return on_segment(Segment<T>{polygon[0], polygon[1]}, point, eps)
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }

    const int order = orientation(
        polygon[0],
        polygon[1],
        polygon[size - 1],
        eps
    );
    if (order == 0) return point_in_polygon(polygon, point, eps);
    auto vertex = [&](std::size_t index) -> const Point<T>& {
        if (order > 0 || index == 0) return polygon[index];
        return polygon[size - index];
    };

    const int first_side = orientation(vertex(0), vertex(1), point, eps);
    const int last_side =
        orientation(vertex(0), vertex(size - 1), point, eps);
    if (first_side < 0 || last_side > 0) {
        return PointInPolygon::Outside;
    }
    if (first_side == 0) {
        return on_segment(Segment<T>{vertex(0), vertex(1)}, point, eps)
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }
    if (last_side == 0) {
        return on_segment(
            Segment<T>{vertex(0), vertex(size - 1)},
            point,
            eps
        )
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }

    std::size_t left = 1;
    std::size_t right = size - 1;
    while (right - left >= 2) {
        const std::size_t middle = (left + right) / 2;
        if (orientation(vertex(0), vertex(middle), point, eps) >= 0) {
            left = middle;
        } else {
            right = middle;
        }
    }
    const int triangle_side =
        orientation(vertex(left), vertex(right), point, eps);
    if (triangle_side < 0) return PointInPolygon::Outside;
    if (triangle_side == 0) return PointInPolygon::Boundary;
    return PointInPolygon::Inside;
}

template <Coordinate T>
class ConvexPolygon {
   public:
    using Wide = wide_type<T>;

   private:
    std::vector<Point<T>> points;
    std::vector<Wide> area_prefix;
    long double epsilon;

    template <class Compare>
    int periodic_best(Compare better) const {
        const int size = int(points.size());
        int left = 0;
        int middle = size;
        int right = 2 * size;
        while (right - left > 2) {
            const int left_middle = (left + middle) / 2;
            const int right_middle = (middle + right + 1) / 2;
            if (better(left_middle % size, middle % size)) {
                right = middle;
                middle = left_middle;
            } else if (better(right_middle % size, middle % size)) {
                left = middle;
                middle = right_middle;
            } else {
                left = left_middle;
                right = right_middle;
            }
        }
        return middle % size;
    }

    int previous(int index) const {
        return index == 0 ? int(points.size()) - 1 : index - 1;
    }

    int next(int index) const {
        return index + 1 == int(points.size()) ? 0 : index + 1;
    }

   public:
    explicit ConvexPolygon(
        std::vector<Point<T>> polygon,
        long double eps = 1e-12L
    )
        : points(normalize_convex_polygon(std::move(polygon), eps)),
          epsilon(eps) {
        assert(
            points.size() <=
            static_cast<std::size_t>(
                std::numeric_limits<int>::max() / 2
            )
        );
        assert(
            points.size() < 3 ||
            is_convex_polygon(points, true, epsilon)
        );
        area_prefix.resize(2 * points.size() + 1, Wide(0));
        for (std::size_t index = 0; index < 2 * points.size(); ++index) {
            area_prefix[index + 1] =
                area_prefix[index] +
                cross(
                    points[index % points.size()],
                    points[(index + 1) % points.size()]
                );
        }
    }

    int size() const noexcept {
        return int(points.size());
    }

    bool empty() const noexcept {
        return points.empty();
    }

    const std::vector<Point<T>>& vertices() const noexcept {
        return points;
    }

    const Point<T>& operator[](int index) const {
        assert(0 <= index && index < size());
        return points[index];
    }

    ConvexPolygon operator+(const ConvexPolygon& other) const {
        const long double eps = std::max(epsilon, other.epsilon);
        return ConvexPolygon(
            minkowski_sum(points, other.points, eps),
            eps
        );
    }

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    ConvexPolygon<std::common_type_t<T, Scalar>> operator*(Scalar scalar) const {
        using Result = std::common_type_t<T, Scalar>;
        std::vector<Point<Result>> scaled;
        scaled.reserve(points.size());
        for (const Point<T>& point : points) {
            scaled.push_back(point * scalar);
        }
        return ConvexPolygon<Result>(std::move(scaled), epsilon);
    }

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    friend ConvexPolygon<std::common_type_t<T, Scalar>> operator*(
        Scalar scalar,
        const ConvexPolygon& polygon
    ) {
        return polygon * scalar;
    }

    Wide area2() const {
        if (points.empty()) return Wide(0);
        return area_prefix[points.size()];
    }

    Wide chain_area2(int first, int last) const {
        assert(0 <= first && first < size());
        assert(0 <= last && last < size());
        int extended_last = last;
        if (extended_last < first) extended_last += size();
        return
            area_prefix[extended_last] - area_prefix[first] +
            cross(points[last], points[first]);
    }

    PointInPolygon contains(const Point<T>& point) const {
        return point_in_convex_polygon(points, point, epsilon);
    }

    std::pair<Wide, int> min_dot(const Point<T>& direction) const {
        assert(!points.empty());
        const int index = periodic_best([&](int first, int second) {
            return dot(points[first], direction) <
                   dot(points[second], direction);
        });
        return std::pair<Wide, int>(dot(points[index], direction), index);
    }

    std::pair<Wide, int> max_dot(const Point<T>& direction) const {
        assert(!points.empty());
        const int index = periodic_best([&](int first, int second) {
            return dot(points[first], direction) >
                   dot(points[second], direction);
        });
        return std::pair<Wide, int>(dot(points[index], direction), index);
    }

    std::pair<int, int> tangent_vertices(const Point<T>& point) const {
        assert(points.size() >= 3);
        assert(contains(point) == PointInPolygon::Outside);
        int first = periodic_best([&](int left, int right) {
            return orientation(point, points[left], points[right], epsilon) < 0;
        });
        int second = periodic_best([&](int left, int right) {
            return orientation(point, points[left], points[right], epsilon) > 0;
        });
        if (
            orientation(
                point,
                points[first],
                points[previous(first)],
                epsilon
            ) == 0
        ) {
            first = previous(first);
        }
        if (
            orientation(
                point,
                points[second],
                points[next(second)],
                epsilon
            ) == 0
        ) {
            second = next(second);
        }
        return std::pair<int, int>(first, second);
    }
};

template <Coordinate T>
std::optional<Point<long double>> centroid(
    const ConvexPolygon<T>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon.vertices(), eps);
}

namespace convex_polygon_detail {

template <Coordinate T>
class MinkowskiDifferenceView {
   private:
    struct Cycle {
        const ConvexPolygon<T>* polygon;
        int start;
        bool negate;

        int edge_count() const {
            return polygon->size() >= 2 ? polygon->size() : 0;
        }

        Point<T> point(int index) const {
            const int size = polygon->size();
            const Point<T>& result = (*polygon)[(start + index) % size];
            return negate ? -result : result;
        }

        Point<T> edge(int index) const {
            return point((index + 1) % polygon->size()) - point(index);
        }
    };

    Cycle first;
    Cycle second;

    std::pair<int, int> prefixes(int rank) const {
        const int first_size = first.edge_count();
        const int second_size = second.edge_count();
        if (first_size + second_size == 0) {
            return std::pair<int, int>(0, 0);
        }

        int low = std::max(0, rank - second_size);
        int high = std::min(rank, first_size);
        while (low <= high) {
            const int first_prefix = (low + high) / 2;
            const int second_prefix = rank - first_prefix;
            if (
                first_prefix > 0 &&
                second_prefix < second_size &&
                entry_less(
                    second.edge(second_prefix),
                    1,
                    first.edge(first_prefix - 1),
                    0
                )
            ) {
                high = first_prefix - 1;
                continue;
            }
            if (
                second_prefix > 0 &&
                first_prefix < first_size &&
                entry_less(
                    first.edge(first_prefix),
                    0,
                    second.edge(second_prefix - 1),
                    1
                )
            ) {
                low = first_prefix + 1;
                continue;
            }
            return std::pair<int, int>(first_prefix, second_prefix);
        }
        assert(false);
        return std::pair<int, int>(0, 0);
    }

    static int direction_half(const Point<T>& direction) {
        return
            direction.y > 0 ||
            (direction.y == 0 && direction.x >= 0)
            ? 0
            : 1;
    }

    static bool entry_less(
        const Point<T>& left,
        int left_cycle,
        const Point<T>& right,
        int right_cycle
    ) {
        if constexpr (std::floating_point<T>) {
            long double left_angle = std::atan2(
                static_cast<long double>(left.y),
                static_cast<long double>(left.x)
            );
            long double right_angle = std::atan2(
                static_cast<long double>(right.y),
                static_cast<long double>(right.x)
            );
            if (left_angle < 0) {
                left_angle += 2 * std::numbers::pi_v<long double>;
            }
            if (right_angle < 0) {
                right_angle += 2 * std::numbers::pi_v<long double>;
            }
            if (left_angle != right_angle) return left_angle < right_angle;
            return left_cycle < right_cycle;
        }
        const int left_half = direction_half(left);
        const int right_half = direction_half(right);
        if (left_half != right_half) return left_half < right_half;
        const auto turn = cross(left, right);
        if (turn != 0) return turn > 0;
        return left_cycle < right_cycle;
    }

    static int negated_start(const ConvexPolygon<T>& polygon) {
        if (polygon.size() <= 1) return 0;
        int result = polygon.max_dot(Point<T>(0, 1)).second;
        const int previous = result == 0 ? polygon.size() - 1 : result - 1;
        const int next = result + 1 == polygon.size() ? 0 : result + 1;
        for (const int candidate : {previous, next}) {
            if (
                polygon[candidate].y == polygon[result].y &&
                polygon[candidate].x > polygon[result].x
            ) {
                result = candidate;
            }
        }
        return result;
    }

   public:
    MinkowskiDifferenceView(
        const ConvexPolygon<T>& minuend,
        const ConvexPolygon<T>& subtrahend
    )
        : first{&minuend, 0, false},
          second{&subtrahend, negated_start(subtrahend), true} {
        assert(!minuend.empty());
        assert(!subtrahend.empty());
    }

    int size() const {
        const int edge_count =
            first.edge_count() + second.edge_count();
        return edge_count == 0 ? 1 : edge_count;
    }

    Point<T> operator[](int rank) const {
        assert(0 <= rank && rank < size());
        const auto [first_prefix, second_prefix] = prefixes(rank);
        return
            first.point(first_prefix % first.polygon->size()) +
            second.point(second_prefix % second.polygon->size());
    }

    std::pair<Point<T>, Point<T>> components(int rank) const {
        assert(0 <= rank && rank < size());
        const auto [first_prefix, second_prefix] = prefixes(rank);
        return std::pair<Point<T>, Point<T>>(
            first.point(first_prefix % first.polygon->size()),
            -second.point(second_prefix % second.polygon->size())
        );
    }
};

struct OriginLocation {
    PointInPolygon location;
    int outside_edge;
    std::array<int, 3> simplex;
    int simplex_size;
};

template <Coordinate T, class Polygon>
OriginLocation locate_origin(
    const Polygon& polygon,
    long double eps
) {
    const int size = polygon.size();
    assert(size >= 3);
    const Point<T> origin;
    const Point<T> base = polygon[0];
    int first = 1;
    if (
        size >= 4 &&
        orientation(base, polygon[1], polygon[2], eps) == 0 &&
        dot(polygon[1] - base, polygon[2] - polygon[1]) > 0
    ) {
        first = 2;
    }
    const int last = size - 1;

    const int first_side = orientation(base, polygon[first], origin, eps);
    const int last_side = orientation(base, polygon[last], origin, eps);
    if (first_side < 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            0,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (last_side > 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            last,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (first_side == 0) {
        if (on_segment(Segment<T>{base, polygon[first]}, origin, eps)) {
            return OriginLocation{
                PointInPolygon::Boundary,
                -1,
                std::array<int, 3>{0, first, 0},
                2,
            };
        }
        return OriginLocation{
            PointInPolygon::Outside,
            first,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (last_side == 0) {
        if (on_segment(Segment<T>{base, polygon[last]}, origin, eps)) {
            return OriginLocation{
                PointInPolygon::Boundary,
                -1,
                std::array<int, 3>{0, last, 0},
                2,
            };
        }
        return OriginLocation{
            PointInPolygon::Outside,
            last - 1,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }

    int left = first;
    int right = last;
    while (right - left >= 2) {
        const int middle = (left + right) / 2;
        if (orientation(base, polygon[middle], origin, eps) >= 0) {
            left = middle;
        } else {
            right = middle;
        }
    }
    const int side = orientation(polygon[left], polygon[right], origin, eps);
    if (side < 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            left,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (side == 0) {
        const bool boundary = on_segment(
            Segment<T>{polygon[left], polygon[right]},
            origin,
            eps
        );
        return OriginLocation{
            boundary ? PointInPolygon::Boundary : PointInPolygon::Outside,
            boundary ? -1 : left,
            std::array<int, 3>{left, right, 0},
            boundary ? 2 : 0,
        };
    }
    return OriginLocation{
        PointInPolygon::Inside,
        -1,
        std::array<int, 3>{0, left, right},
        3,
    };
}

template <class Compare>
int periodic_best(int size, Compare better) {
    int left = 0;
    int middle = size;
    int right = 2 * size;
    while (right - left > 2) {
        const int left_middle = (left + middle) / 2;
        const int right_middle = (middle + right + 1) / 2;
        if (better(left_middle % size, middle % size)) {
            right = middle;
            middle = left_middle;
        } else if (better(right_middle % size, middle % size)) {
            left = middle;
            middle = right_middle;
        } else {
            left = left_middle;
            right = right_middle;
        }
    }
    return middle % size;
}

template <Coordinate T, class Polygon>
std::pair<int, int> tangent_vertices_from_origin(
    const Polygon& polygon,
    long double eps
) {
    const int size = polygon.size();
    const Point<T> origin;
    int first = periodic_best(size, [&](int left, int right) {
        return orientation(origin, polygon[left], polygon[right], eps) < 0;
    });
    int second = periodic_best(size, [&](int left, int right) {
        return orientation(origin, polygon[left], polygon[right], eps) > 0;
    });
    const int previous = first == 0 ? size - 1 : first - 1;
    if (orientation(origin, polygon[first], polygon[previous], eps) == 0) {
        first = previous;
    }
    const int next = second + 1 == size ? 0 : second + 1;
    if (orientation(origin, polygon[second], polygon[next], eps) == 0) {
        second = next;
    }
    return std::pair<int, int>(first, second);
}

struct ClosestBoundaryFeature {
    int first;
    int second;
    long double ratio;
    long double distance;
};

template <Coordinate T, class Polygon>
ClosestBoundaryFeature closest_boundary_feature(
    const Polygon& polygon,
    const OriginLocation& location,
    long double eps
) {
    const int size = polygon.size();
    assert(size >= 3);
    assert(location.location == PointInPolygon::Outside);
    const Point<T> origin;

    const auto tangents = tangent_vertices_from_origin<T>(polygon, eps);
    auto visible = [&](int index) {
        return orientation(
            polygon[index],
            polygon[(index + 1) % size],
            origin,
            eps
        ) < 0;
    };
    auto forward_edges = [&](int start, int finish) {
        return finish >= start ? finish - start : finish + size - start;
    };

    int witness = location.outside_edge;
    if (!visible(witness)) {
        const int previous = witness == 0 ? size - 1 : witness - 1;
        const int next = witness + 1 == size ? 0 : witness + 1;
        if (visible(previous)) {
            witness = previous;
        } else if (visible(next)) {
            witness = next;
        }
    }

    int start = tangents.first;
    int finish = tangents.second;
    if (forward_edges(start, witness) >= forward_edges(start, finish)) {
        std::swap(start, finish);
    }
    int edge_count = forward_edges(start, finish);
    if (edge_count == 0) {
        start = location.outside_edge;
        finish = (start + 1) % size;
        edge_count = 1;
    }

    auto vertex = [&](int offset) {
        return polygon[(start + offset) % size];
    };
    int left = 0;
    int right = edge_count;
    while (left < right) {
        const int middle = (left + right) / 2;
        if (norm2(vertex(middle)) <= norm2(vertex(middle + 1))) {
            right = middle;
        } else {
            left = middle + 1;
        }
    }

    ClosestBoundaryFeature result{
        (start + left) % size,
        (start + left) % size,
        0,
        norm(vertex(left)),
    };
    auto consider_edge = [&](int first_offset, int second_offset) {
        const Point<long double> first_point(vertex(first_offset));
        const Point<long double> second_point(vertex(second_offset));
        const Point<long double> direction = second_point - first_point;
        long double ratio =
            -dot(first_point, direction) / dot(direction, direction);
        ratio = std::clamp(ratio, 0.0L, 1.0L);
        const long double candidate_distance =
            norm(first_point + direction * ratio);
        if (candidate_distance < result.distance) {
            result = ClosestBoundaryFeature{
                (start + first_offset) % size,
                (start + second_offset) % size,
                ratio,
                candidate_distance,
            };
        }
    };
    if (left > 0) {
        consider_edge(left - 1, left);
    }
    if (left < edge_count) {
        consider_edge(left, left + 1);
    }
    return result;
}

template <Coordinate T, class Polygon>
long double distance_from_origin(
    const Polygon& polygon,
    long double eps
) {
    const OriginLocation location = locate_origin<T>(polygon, eps);
    if (location.location != PointInPolygon::Outside) return 0;
    return closest_boundary_feature<T>(polygon, location, eps).distance;
}

inline Point<long double> interpolate(
    const Point<long double>& first,
    const Point<long double>& second,
    long double ratio
) {
    return first + (second - first) * ratio;
}

template <Coordinate T>
ClosestPoints
closest_points_from_difference(
    const MinkowskiDifferenceView<T>& difference,
    long double eps
) {
    const OriginLocation location = locate_origin<T>(difference, eps);
    if (location.location == PointInPolygon::Outside) {
        const ClosestBoundaryFeature feature =
            closest_boundary_feature<T>(difference, location, eps);
        const auto first_components = difference.components(feature.first);
        const auto second_components = difference.components(feature.second);
        return ClosestPoints{
            interpolate(
                Point<long double>(first_components.first),
                Point<long double>(second_components.first),
                feature.ratio
            ),
            interpolate(
                Point<long double>(first_components.second),
                Point<long double>(second_components.second),
                feature.ratio
            )
        };
    }

    assert(location.simplex_size == 2 || location.simplex_size == 3);
    std::array<long double, 3> weight{0, 0, 0};
    if (location.simplex_size == 2) {
        const Point<long double> first(difference[location.simplex[0]]);
        const Point<long double> second(difference[location.simplex[1]]);
        const Point<long double> direction = second - first;
        weight[1] = -dot(first, direction) / dot(direction, direction);
        weight[1] = std::clamp(weight[1], 0.0L, 1.0L);
        weight[0] = 1 - weight[1];
    } else {
        const Point<long double> first(difference[location.simplex[0]]);
        const Point<long double> second(difference[location.simplex[1]]);
        const Point<long double> third(difference[location.simplex[2]]);
        const long double denominator = cross(
            second - first,
            third - first
        );
        weight[0] = cross(second, third) / denominator;
        weight[1] = cross(third, first) / denominator;
        weight[2] = cross(first, second) / denominator;
    }

    Point<long double> first_result;
    Point<long double> second_result;
    for (int index = 0; index < location.simplex_size; ++index) {
        const auto components = difference.components(
            location.simplex[index]
        );
        first_result += Point<long double>(components.first) * weight[index];
        second_result +=
            Point<long double>(components.second) * weight[index];
    }
    return ClosestPoints{
        first_result,
        second_result
    };
}

}  // namespace convex_polygon_detail

template <Coordinate T>
std::vector<std::array<Point<T>, 3>> triangulate_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    polygon = normalize_convex_polygon(std::move(polygon), eps);
    if (polygon.size() < 3) return {};

    std::vector<std::array<Point<T>, 3>> result;
    result.reserve(polygon.size() - 2);
    for (std::size_t index = 1; index + 1 < polygon.size(); ++index) {
        std::array<Point<T>, 3> triangle;
        triangle[0] = polygon[0];
        triangle[1] = polygon[index];
        triangle[2] = polygon[index + 1];
        result.push_back(std::move(triangle));
    }
    return result;
}

template <Coordinate T>
wide_type<T> convex_diameter2(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    polygon = normalize_convex_polygon(std::move(polygon), eps);
    const std::size_t size = polygon.size();
    if (size <= 1) return 0;
    if (size == 2) return distance2(polygon[1], polygon[0]);

    wide_type<T> result = 0;
    std::size_t opposite = 1;
    for (std::size_t index = 0; index < size; ++index) {
        const std::size_t next = (index + 1) % size;
        while (true) {
            const std::size_t candidate = (opposite + 1) % size;
            const auto current_area =
                cross(polygon[index], polygon[next], polygon[opposite]);
            const auto candidate_area =
                cross(polygon[index], polygon[next], polygon[candidate]);
            if (candidate_area <= current_area) break;
            opposite = candidate;
        }
        result = std::max(
            result,
            distance2(polygon[index], polygon[opposite])
        );
        result = std::max(
            result,
            distance2(polygon[next], polygon[opposite])
        );
    }
    return result;
}

template <Coordinate T>
std::vector<Point<long double>> convex_cut(
    const std::vector<Point<T>>& polygon,
    const Line<T>& boundary,
    long double eps = 1e-12L
) {
    assert(boundary.a != boundary.b);
    std::vector<Point<long double>> input;
    input.reserve(polygon.size());
    for (const Point<T>& point : polygon) input.emplace_back(point);
    if (input.empty()) return input;

    const Point<long double> line_start(boundary.a);
    const Point<long double> line_end(boundary.b);
    const Line<long double> line{line_start, line_end};
    std::vector<Point<long double>> result;
    Point<long double> previous = input.back();
    int previous_side = orientation(line_start, line_end, previous, eps);
    for (const Point<long double>& current : input) {
        const int current_side =
            orientation(line_start, line_end, current, eps);
        const bool previous_inside = previous_side >= 0;
        const bool current_inside = current_side >= 0;
        if (previous_inside != current_inside) {
            const Line<long double> crossing{previous, current};
            const LinearIntersection intersection =
                linear_intersection(line, crossing, eps);
            if (intersection.kind == LinearIntersectionKind::Point) {
                result.push_back(intersection.first);
            }
        }
        if (current_inside) result.push_back(current);
        previous = current;
        previous_side = current_side;
    }
    return convex_polygon_detail::clean_polygon(std::move(result), eps);
}

template <Coordinate T>
bool convex_polygons_intersect(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        if (first.size() == 1 && second.size() == 1) {
            return distance(first[0], second[0]) <= eps;
        }
        if (first.size() == 1) {
            return on_segment(
                Segment<T>{second[0], second[1]},
                first[0],
                eps
            );
        }
        if (second.size() == 1) {
            return on_segment(
                Segment<T>{first[0], first[1]},
                second[0],
                eps
            );
        }
        return intersects(
            Segment<T>{first[0], first[1]},
            Segment<T>{second[0], second[1]},
            eps
        );
    }

    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return
        convex_polygon_detail::locate_origin<T>(difference, eps).location !=
        PointInPolygon::Outside;
}

template <Coordinate T>
bool convex_polygons_intersect(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    std::vector<Point<T>> negated;
    negated.reserve(second.size());
    for (const Point<T>& point : second) negated.push_back(-point);
    const std::vector<Point<T>> difference =
        minkowski_sum(first, std::move(negated), eps);
    return
        point_in_convex_polygon(difference, Point<T>(), eps) !=
        PointInPolygon::Outside;
}

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        return closest_points(
            Segment<T>{first[0], first[first.size() - 1]},
            Segment<T>{second[0], second[second.size() - 1]},
            eps
        );
    }
    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return convex_polygon_detail::closest_points_from_difference(
        difference,
        eps
    );
}

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    const ConvexPolygon<T> first_query(first, eps);
    const ConvexPolygon<T> second_query(second, eps);
    return convex_polygons_closest_points(first_query, second_query, eps);
}

template <Coordinate T>
long double convex_polygons_distance(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        if (convex_polygons_intersect(first, second, eps)) return 0;
        if (first.size() == 1 && second.size() == 1) {
            return distance(first[0], second[0]);
        }
        if (first.size() == 1) {
            return distance(
                Segment<T>{second[0], second[1]},
                first[0]
            );
        }
        if (second.size() == 1) {
            return distance(
                Segment<T>{first[0], first[1]},
                second[0]
            );
        }
        return distance(
            Segment<T>{first[0], first[1]},
            Segment<T>{second[0], second[1]}
        );
    }

    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return convex_polygon_detail::distance_from_origin<T>(difference, eps);
}

template <Coordinate T>
std::vector<Point<long double>> convex_polygon_intersection(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    using HalfPlane = half_plane_intersection_detail::HalfPlane;
    namespace detail = half_plane_intersection_detail;

    const std::vector<Point<T>> normalized_first =
        normalize_convex_polygon(first, eps);
    const std::vector<Point<T>> normalized_second =
        normalize_convex_polygon(second, eps);
    assert(normalized_first.size() >= 3);
    assert(normalized_second.size() >= 3);
    assert(is_convex_polygon(normalized_first, true, eps));
    assert(is_convex_polygon(normalized_second, true, eps));
    if (!convex_polygons_intersect(
            normalized_first,
            normalized_second,
            eps
        )) {
        return {};
    }

    auto boundaries = [](const std::vector<Point<T>>& polygon) {
        std::vector<HalfPlane> result;
        result.reserve(polygon.size());
        for (std::size_t index = 0; index < polygon.size(); ++index) {
            const Point<long double> point(polygon[index]);
            Point<long double> direction =
                Point<long double>(polygon[(index + 1) % polygon.size()]) -
                point;
            direction = direction / norm(direction);
            result.push_back(HalfPlane{point, direction});
        }
        return result;
    };
    const std::vector<HalfPlane> first_boundaries =
        boundaries(normalized_first);
    const std::vector<HalfPlane> second_boundaries =
        boundaries(normalized_second);

    std::vector<HalfPlane> merged;
    merged.reserve(first_boundaries.size() + second_boundaries.size());
    std::size_t first_index = 0;
    std::size_t second_index = 0;
    while (
        first_index < first_boundaries.size() ||
        second_index < second_boundaries.size()
    ) {
        const bool take_first =
            second_index == second_boundaries.size() ||
            (
                first_index < first_boundaries.size() &&
                detail::direction_less(
                    first_boundaries[first_index],
                    second_boundaries[second_index]
                )
            );
        if (take_first) {
            detail::merge_same_direction(
                merged,
                first_boundaries[first_index++],
                eps
            );
        } else {
            detail::merge_same_direction(
                merged,
                second_boundaries[second_index++],
                eps
            );
        }
    }
    detail::merge_cyclic_ends(merged, eps);

    std::deque<HalfPlane> active;
    for (const HalfPlane& half_plane : merged) {
        while (active.size() >= 2) {
            const std::optional<Point<long double>> point =
                detail::intersection(
                    active[active.size() - 2],
                    active.back(),
                    eps
                );
            if (
                !point.has_value() ||
                !detail::outside(half_plane, *point, eps)
            ) {
                break;
            }
            active.pop_back();
        }
        while (active.size() >= 2) {
            const std::optional<Point<long double>> point =
                detail::intersection(active[0], active[1], eps);
            if (
                !point.has_value() ||
                !detail::outside(half_plane, *point, eps)
            ) {
                break;
            }
            active.pop_front();
        }
        active.push_back(half_plane);
    }
    while (active.size() >= 3) {
        const std::optional<Point<long double>> point =
            detail::intersection(
                active[active.size() - 2],
                active.back(),
                eps
            );
        if (
            !point.has_value() ||
            !detail::outside(active.front(), *point, eps)
        ) {
            break;
        }
        active.pop_back();
    }
    while (active.size() >= 3) {
        const std::optional<Point<long double>> point =
            detail::intersection(active[0], active[1], eps);
        if (
            !point.has_value() ||
            !detail::outside(active.back(), *point, eps)
        ) {
            break;
        }
        active.pop_front();
    }

    std::vector<Point<long double>> result;
    result.reserve(active.size());
    for (std::size_t index = 0; index < active.size(); ++index) {
        const std::optional<Point<long double>> point =
            detail::intersection(
                active[index],
                active[(index + 1) % active.size()],
                eps
            );
        if (point.has_value()) result.push_back(*point);
    }
    return convex_polygon_detail::clean_polygon(std::move(result), eps);
}

template <Coordinate T>
long double convex_polygons_distance(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    std::vector<Point<T>> negated;
    negated.reserve(second.size());
    for (const Point<T>& point : second) negated.push_back(-point);
    const std::vector<Point<T>> difference =
        minkowski_sum(first, std::move(negated), eps);
    const Point<T> origin;
    if (
        point_in_convex_polygon(difference, origin, eps) !=
        PointInPolygon::Outside
    ) {
        return 0;
    }
    if (difference.size() == 1) return distance(difference[0], origin);

    long double result = std::numeric_limits<long double>::infinity();
    for (std::size_t index = 0; index < difference.size(); ++index) {
        if (difference.size() == 2 && index == 1) break;
        result = std::min(
            result,
            distance(
                Segment<T>{
                    difference[index],
                    difference[(index + 1) % difference.size()]
                },
                origin
            )
        );
    }
    return result;
}

}  // namespace geometry
}  // namespace m1une

#endif  // M1UNE_GEOMETRY_CONVEX_POLYGON_HPP

#line 1 "geometry/convex_polygon.hpp"



#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <deque>
#include <limits>
#include <numbers>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>

#line 1 "geometry/convex_hull.hpp"



#line 8 "geometry/convex_hull.hpp"

#line 1 "geometry/point.hpp"



#line 8 "geometry/point.hpp"

#line 1 "geometry/detail/floating_predicate.hpp"



namespace m1une {
namespace geometry {
namespace predicate_detail {

template <typename T>
constexpr T absolute(T value) {
    return value < T(0) ? -value : value;
}

template <typename T>
constexpr T max_value(T first, T second) {
    return first < second ? second : first;
}

template <typename T>
constexpr T vector_scale(T x, T y) {
    return max_value(absolute(x), absolute(y));
}

template <bool Exact, typename T>
constexpr int scaled_sign(T value, T scale, long double eps) {
    if constexpr (Exact) {
        return (value > T(0)) - (value < T(0));
    } else {
        const T tolerance = T(eps) * scale;
        return (value > tolerance) - (value < -tolerance);
    }
}

template <bool Exact, typename T>
constexpr T determinant_scale(T ax, T ay, T bx, T by) {
    if constexpr (Exact) {
        return T(0);
    } else {
        return vector_scale(ax, ay) * vector_scale(bx, by);
    }
}

template <bool Exact, typename T>
constexpr int determinant_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T determinant = ax * by - ay * bx;
    return scaled_sign<Exact>(
        determinant,
        determinant_scale<Exact>(ax, ay, bx, by),
        eps
    );
}

template <bool Exact, typename T>
constexpr int orientation_sign(
    T direction_x,
    T direction_y,
    T offset_x,
    T offset_y,
    long double eps
) {
    const T determinant =
        direction_x * offset_y - direction_y * offset_x;
    T scale = T(0);
    if constexpr (!Exact) {
        const T direction_scale =
            vector_scale(direction_x, direction_y);
        scale = direction_scale * max_value(
            direction_scale,
            vector_scale(offset_x, offset_y)
        );
    }
    return scaled_sign<Exact>(determinant, scale, eps);
}

template <bool Exact, typename T>
constexpr int dot_sign(
    T ax,
    T ay,
    T bx,
    T by,
    long double eps
) {
    const T value = ax * bx + ay * by;
    T scale = T(0);
    if constexpr (!Exact) {
        scale = vector_scale(ax, ay) * vector_scale(bx, by);
    }
    return scaled_sign<Exact>(value, scale, eps);
}

}  // namespace predicate_detail
}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/point.hpp"

namespace m1une {
namespace geometry {

template <typename T>
concept Coordinate = !std::same_as<std::remove_cv_t<T>, bool> &&
    (std::is_arithmetic_v<T> ||
     (std::copyable<T> && std::totally_ordered<T> && requires(T a, T b) {
         T(0);
         T(1);
         static_cast<long double>(a);
         { +a } -> std::same_as<T>;
         { -a } -> std::same_as<T>;
         { a + b } -> std::same_as<T>;
         { a - b } -> std::same_as<T>;
         { a * b } -> std::same_as<T>;
         { a / b } -> std::same_as<T>;
         { a += b } -> std::same_as<T&>;
         { a -= b } -> std::same_as<T&>;
     }));

// Custom coordinate types keep their own exact arithmetic.
template <typename T>
concept ExactCoordinate = Coordinate<T> && !std::floating_point<T>;

template <Coordinate T>
using wide_type = std::conditional_t<std::integral<T>, __int128_t,
    std::conditional_t<std::floating_point<T>, long double, T>>;

template <Coordinate T>
struct Point {
    T x;
    T y;

    constexpr Point() : x(0), y(0) {}
    constexpr Point(T x_value, T y_value) : x(x_value), y(y_value) {}

    template <Coordinate U>
    explicit constexpr Point(const Point<U>& other)
        : x(static_cast<T>(other.x)), y(static_cast<T>(other.y)) {}

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

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

    constexpr Point operator+() const {
        return *this;
    }

    constexpr Point operator-() const {
        return Point(-x, -y);
    }

    friend constexpr Point operator+(Point left, const Point& right) {
        return left += right;
    }

    friend constexpr Point operator-(Point left, const Point& right) {
        return left -= right;
    }

    friend constexpr bool operator==(const Point&, const Point&) = default;

    friend constexpr bool operator<(const Point& left, const Point& right) {
        if (left.x != right.x) return left.x < right.x;
        return left.y < right.y;
    }
};

template <Coordinate T>
constexpr Point<long double> centroid(const Point<T>& point) {
    return Point<long double>(point);
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) * Result(scalar),
        Result(point.y) * Result(scalar)
    );
}

template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator*(Scalar scalar, const Point<T>& point) {
    return point * scalar;
}

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
constexpr auto operator/(const Point<T>& point, Scalar scalar) {
    using Result = std::common_type_t<T, Scalar>;
    return Point<Result>(
        Result(point.x) / Result(scalar),
        Result(point.y) / Result(scalar)
    );
}

template <Coordinate T>
constexpr wide_type<T> dot(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.x) + W(a.y) * W(b.y);
}

template <Coordinate T>
constexpr wide_type<T> cross(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    return W(a.x) * W(b.y) - W(a.y) * W(b.x);
}

template <Coordinate T>
constexpr wide_type<T> cross(
    const Point<T>& origin,
    const Point<T>& a,
    const Point<T>& b
) {
    using W = wide_type<T>;
    W ax = W(a.x) - W(origin.x);
    W ay = W(a.y) - W(origin.y);
    W bx = W(b.x) - W(origin.x);
    W by = W(b.y) - W(origin.y);
    return ax * by - ay * bx;
}

template <Coordinate T>
constexpr wide_type<T> norm2(const Point<T>& point) {
    return dot(point, point);
}

template <Coordinate T>
constexpr wide_type<T> distance2(const Point<T>& a, const Point<T>& b) {
    using W = wide_type<T>;
    W dx = W(a.x) - W(b.x);
    W dy = W(a.y) - W(b.y);
    return dx * dx + dy * dy;
}

template <Coordinate T>
long double norm(const Point<T>& point) {
    return std::hypot(
        static_cast<long double>(point.x),
        static_cast<long double>(point.y)
    );
}

template <Coordinate T>
long double distance(const Point<T>& a, const Point<T>& b) {
    return std::hypot(
        static_cast<long double>(a.x) - static_cast<long double>(b.x),
        static_cast<long double>(a.y) - static_cast<long double>(b.y)
    );
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> internal_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio + second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T, typename M, typename N>
requires (std::is_arithmetic_v<M> || Coordinate<M>) &&
         (std::is_arithmetic_v<N> || Coordinate<N>)
constexpr Point<long double> external_division_point(
    const Point<T>& a,
    const Point<T>& b,
    M m,
    N n
) {
    long double first_ratio = static_cast<long double>(m);
    long double second_ratio = static_cast<long double>(n);
    long double denominator = first_ratio - second_ratio;
    assert(denominator != 0);
    Point<long double> first(a);
    Point<long double> direction = Point<long double>(b) - first;
    return first + direction * (first_ratio / denominator);
}

template <Coordinate T>
constexpr int sign(wide_type<T> value, long double eps = 1e-12L) {
    return predicate_detail::scaled_sign<ExactCoordinate<T>>(
        value,
        wide_type<T>(1),
        eps
    );
}

template <Coordinate T>
constexpr int orientation(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    using W = wide_type<T>;
    const W first_x = W(b.x) - W(a.x);
    const W first_y = W(b.y) - W(a.y);
    const W second_x = W(c.x) - W(a.x);
    const W second_y = W(c.y) - W(a.y);
    return predicate_detail::orientation_sign<ExactCoordinate<T>>(
        first_x,
        first_y,
        second_x,
        second_y,
        eps
    );
}

template <Coordinate T>
constexpr bool collinear(
    const Point<T>& a,
    const Point<T>& b,
    const Point<T>& c,
    long double eps = 1e-12L
) {
    return orientation(a, b, c, eps) == 0;
}

template <Coordinate T>
Point<long double> rotate(const Point<T>& point, long double angle) {
    long double cosine = std::cos(angle);
    long double sine = std::sin(angle);
    return Point<long double>(
        static_cast<long double>(point.x) * cosine -
            static_cast<long double>(point.y) * sine,
        static_cast<long double>(point.x) * sine +
            static_cast<long double>(point.y) * cosine
    );
}

template <Coordinate T>
Point<long double> normalized(const Point<T>& point) {
    long double length = norm(point);
    assert(length != 0);
    return Point<long double>(
        static_cast<long double>(point.x) / length,
        static_cast<long double>(point.y) / length
    );
}

}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/convex_hull.hpp"

namespace m1une {
namespace geometry {

// Returns the convex hull counterclockwise from its lexicographically smallest
// point. The first point is not repeated at the end.
template <Coordinate T>
std::vector<Point<T>> convex_hull(
    std::vector<Point<T>> points,
    bool include_collinear = false
) {
    std::sort(points.begin(), points.end());
    points.erase(std::unique(points.begin(), points.end()), points.end());
    std::size_t size = points.size();
    if (size <= 1) return points;

    std::vector<Point<T>> hull;
    hull.reserve(2 * size);
    auto should_pop = [include_collinear](
        const Point<T>& first,
        const Point<T>& second,
        const Point<T>& third
    ) {
        int turn = orientation(first, second, third);
        return include_collinear ? turn < 0 : turn <= 0;
    };

    for (const Point<T>& point : points) {
        while (
            hull.size() >= 2 &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }

    std::size_t lower_size = hull.size();
    for (std::size_t index = size - 1; index-- > 0;) {
        const Point<T>& point = points[index];
        while (
            hull.size() > lower_size &&
            should_pop(hull[hull.size() - 2], hull.back(), point)
        ) {
            hull.pop_back();
        }
        hull.push_back(point);
    }
    hull.pop_back();

    if (include_collinear && hull.size() == 2 * points.size() - 2) {
        hull = std::move(points);
    }
    return hull;
}

}  // namespace geometry
}  // namespace m1une


#line 1 "geometry/half_plane_intersection.hpp"



#line 12 "geometry/half_plane_intersection.hpp"
#include <random>
#line 15 "geometry/half_plane_intersection.hpp"

#line 1 "geometry/linear.hpp"



#line 7 "geometry/linear.hpp"

#line 9 "geometry/linear.hpp"

namespace m1une {
namespace geometry {

template <Coordinate T>
struct Line {
    Point<T> a;
    Point<T> b;
};

template <Coordinate T>
struct Segment {
    Point<T> a;
    Point<T> b;
};

template <Coordinate T>
struct Ray {
    Point<T> origin;
    Point<T> through;
};

enum class LinearIntersectionKind {
    Empty,
    Point,
    Segment,
    Ray,
    Line,
};

struct LinearIntersection {
    LinearIntersectionKind kind;
    Point<long double> first;
    Point<long double> second;
};

struct ClosestPoints {
    Point<long double> first;
    Point<long double> second;
};

namespace linear_intersection_detail {

inline LinearIntersection make_empty() {
    const Point<long double> zero;
    return LinearIntersection{
        LinearIntersectionKind::Empty,
        zero,
        zero,
    };
}

template <Coordinate T>
LinearIntersection make_point(const Point<T>& point) {
    const Point<long double> converted(point);
    return LinearIntersection{
        LinearIntersectionKind::Point,
        converted,
        converted,
    };
}

template <Coordinate T>
LinearIntersection make_object(
    LinearIntersectionKind kind,
    const Point<T>& first,
    const Point<T>& second
) {
    return LinearIntersection{
        kind,
        Point<long double>(first),
        Point<long double>(second),
    };
}

}  // namespace linear_intersection_detail

template <Coordinate T>
constexpr Point<long double> centroid(const Segment<T>& segment) {
    return Point<long double>(
        (
            static_cast<long double>(segment.a.x) +
            static_cast<long double>(segment.b.x)
        ) / 2,
        (
            static_cast<long double>(segment.a.y) +
            static_cast<long double>(segment.b.y)
        ) / 2
    );
}

template <Coordinate T>
bool on_line(
    const Line<T>& line,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    assert(line.a != line.b);
    return orientation(line.a, line.b, point, eps) == 0;
}

template <Coordinate T>
bool parallel(const Line<T>& first, const Line<T>& second, long double eps = 1e-12L) {
    using W = wide_type<T>;
    W first_x = W(first.b.x) - W(first.a.x);
    W first_y = W(first.b.y) - W(first.a.y);
    W second_x = W(second.b.x) - W(second.a.x);
    W second_y = W(second.b.y) - W(second.a.y);
    return predicate_detail::determinant_sign<ExactCoordinate<T>>(
        first_x,
        first_y,
        second_x,
        second_y,
        eps
    ) == 0;
}

template <Coordinate T>
bool orthogonal(const Line<T>& first, const Line<T>& second, long double eps = 1e-12L) {
    using W = wide_type<T>;
    W first_x = W(first.b.x) - W(first.a.x);
    W first_y = W(first.b.y) - W(first.a.y);
    W second_x = W(second.b.x) - W(second.a.x);
    W second_y = W(second.b.y) - W(second.a.y);
    return predicate_detail::dot_sign<ExactCoordinate<T>>(
        first_x,
        first_y,
        second_x,
        second_y,
        eps
    ) == 0;
}

template <Coordinate T>
Point<long double> projection(const Line<T>& line, const Point<T>& point) {
    assert(line.a != line.b);
    Point<long double> a(line.a);
    Point<long double> direction(
        static_cast<long double>(line.b.x) - static_cast<long double>(line.a.x),
        static_cast<long double>(line.b.y) - static_cast<long double>(line.a.y)
    );
    Point<long double> offset(
        static_cast<long double>(point.x) - a.x,
        static_cast<long double>(point.y) - a.y
    );
    long double ratio = dot(offset, direction) / dot(direction, direction);
    return a + direction * ratio;
}

template <Coordinate T>
Point<long double> reflection(const Line<T>& line, const Point<T>& point) {
    Point<long double> projected = projection(line, point);
    return projected * 2.0L - Point<long double>(point);
}

template <Coordinate T>
bool intersects(
    const Line<T>& first,
    const Line<T>& second,
    long double eps = 1e-12L
) {
    return !parallel(first, second, eps) || on_line(first, second.a, eps);
}

template <Coordinate T>
bool on_segment(
    const Segment<T>& segment,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    if (orientation(segment.a, segment.b, point, eps) != 0) return false;
    using W = wide_type<T>;
    const W direction_x = W(segment.b.x) - W(segment.a.x);
    const W direction_y = W(segment.b.y) - W(segment.a.y);
    if (direction_x == W(0) && direction_y == W(0)) {
        if constexpr (ExactCoordinate<T>) {
            return point == segment.a;
        } else {
            return
                predicate_detail::absolute(W(point.x) - W(segment.a.x)) <= eps &&
                predicate_detail::absolute(W(point.y) - W(segment.a.y)) <= eps;
        }
    }
    const W offset_x = W(point.x) - W(segment.a.x);
    const W offset_y = W(point.y) - W(segment.a.y);
    const W projection =
        offset_x * direction_x + offset_y * direction_y;
    const W length_squared =
        direction_x * direction_x + direction_y * direction_y;
    return
        predicate_detail::scaled_sign<ExactCoordinate<T>>(
            projection,
            length_squared,
            eps
        ) >= 0 &&
        predicate_detail::scaled_sign<ExactCoordinate<T>>(
            projection - length_squared,
            length_squared,
            eps
        ) <= 0;
}

template <Coordinate T>
Point<long double> projection(
    const Segment<T>& segment,
    const Point<T>& point
) {
    const Point<long double> first(segment.a);
    const Point<long double> direction =
        Point<long double>(segment.b) - first;
    const long double length_squared = dot(direction, direction);
    if (length_squared == 0) return first;
    const long double ratio = std::clamp(
        dot(Point<long double>(point) - first, direction) / length_squared,
        0.0L,
        1.0L
    );
    return first + direction * ratio;
}

template <Coordinate T>
bool intersects(
    const Segment<T>& first,
    const Segment<T>& second,
    long double eps = 1e-12L
) {
    int abc = orientation(first.a, first.b, second.a, eps);
    int abd = orientation(first.a, first.b, second.b, eps);
    int cda = orientation(second.a, second.b, first.a, eps);
    int cdb = orientation(second.a, second.b, first.b, eps);

    if (abc == 0 && on_segment(first, second.a, eps)) return true;
    if (abd == 0 && on_segment(first, second.b, eps)) return true;
    if (cda == 0 && on_segment(second, first.a, eps)) return true;
    if (cdb == 0 && on_segment(second, first.b, eps)) return true;
    return abc * abd < 0 && cda * cdb < 0;
}

template <Coordinate T>
bool intersects(
    const Line<T>& line,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    int first_side = orientation(line.a, line.b, segment.a, eps);
    int second_side = orientation(line.a, line.b, segment.b, eps);
    return first_side == 0 || second_side == 0 || first_side != second_side;
}

template <Coordinate T>
bool intersects(
    const Segment<T>& segment,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    return intersects(line, segment, eps);
}

namespace linear_parameter_detail {

template <Coordinate T>
struct Parameters {
    wide_type<T> denominator;
    wide_type<T> denominator_scale;
    wide_type<T> first_numerator;
    wide_type<T> second_numerator;
};

template <Coordinate T>
Parameters<T> parameters(
    const Point<T>& first_origin,
    const Point<T>& first_through,
    const Point<T>& second_origin,
    const Point<T>& second_through
) {
    using W = wide_type<T>;
    W first_x = W(first_through.x) - W(first_origin.x);
    W first_y = W(first_through.y) - W(first_origin.y);
    W second_x = W(second_through.x) - W(second_origin.x);
    W second_y = W(second_through.y) - W(second_origin.y);
    W offset_x = W(second_origin.x) - W(first_origin.x);
    W offset_y = W(second_origin.y) - W(first_origin.y);
    return Parameters<T>{
        first_x * second_y - first_y * second_x,
        predicate_detail::determinant_scale<ExactCoordinate<T>>(
            first_x,
            first_y,
            second_x,
            second_y
        ),
        offset_x * second_y - offset_y * second_x,
        offset_x * first_y - offset_y * first_x
    };
}

template <Coordinate T>
int denominator_sign(const Parameters<T>& values, long double eps) {
    return predicate_detail::scaled_sign<ExactCoordinate<T>>(
        values.denominator,
        values.denominator_scale,
        eps
    );
}

template <Coordinate T>
bool ratio_nonnegative(
    wide_type<T> numerator,
    wide_type<T> denominator,
    long double eps
) {
    const int numerator_sign =
        predicate_detail::scaled_sign<ExactCoordinate<T>>(
            numerator,
            predicate_detail::absolute(denominator),
            eps
        );
    const int denominator_direction =
        (denominator > 0) - (denominator < 0);
    return
        numerator_sign == 0 ||
        numerator_sign == denominator_direction;
}

template <Coordinate T>
bool ratio_in_unit_interval(
    wide_type<T> numerator,
    wide_type<T> denominator,
    long double eps
) {
    const auto scale = predicate_detail::absolute(denominator);
    const int start_sign =
        predicate_detail::scaled_sign<ExactCoordinate<T>>(
            numerator,
            scale,
            eps
        );
    const int finish_sign =
        predicate_detail::scaled_sign<ExactCoordinate<T>>(
            numerator - denominator,
            scale,
            eps
        );
    if (denominator > 0) {
        return start_sign >= 0 && finish_sign <= 0;
    }
    return start_sign <= 0 && finish_sign >= 0;
}

}  // namespace linear_parameter_detail

template <Coordinate T>
bool on_ray(
    const Ray<T>& ray,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    assert(ray.origin != ray.through);
    if (orientation(ray.origin, ray.through, point, eps) != 0) return false;
    using W = wide_type<T>;
    W direction_x = W(ray.through.x) - W(ray.origin.x);
    W direction_y = W(ray.through.y) - W(ray.origin.y);
    W offset_x = W(point.x) - W(ray.origin.x);
    W offset_y = W(point.y) - W(ray.origin.y);
    const W projection =
        direction_x * offset_x + direction_y * offset_y;
    const W length_squared =
        direction_x * direction_x + direction_y * direction_y;
    return predicate_detail::scaled_sign<ExactCoordinate<T>>(
        projection,
        length_squared,
        eps
    ) >= 0;
}

template <Coordinate T>
Point<long double> projection(const Ray<T>& ray, const Point<T>& point) {
    assert(ray.origin != ray.through);
    Point<long double> origin(ray.origin);
    Point<long double> direction =
        Point<long double>(ray.through) - origin;
    Point<long double> offset = Point<long double>(point) - origin;
    long double ratio = dot(offset, direction) / dot(direction, direction);
    if (ratio < 0) ratio = 0;
    return origin + direction * ratio;
}

template <Coordinate T>
Ray<long double> reflection(const Line<T>& line, const Ray<T>& ray) {
    assert(ray.origin != ray.through);
    return Ray<long double>{
        reflection(line, ray.origin),
        reflection(line, ray.through)
    };
}

template <Coordinate T>
Ray<long double> reflected_ray(
    const Ray<T>& incoming,
    const Point<T>& hit,
    const Line<T>& mirror,
    long double eps = 1e-12L
) {
    assert(incoming.origin != incoming.through);
    assert(on_line(mirror, hit, eps));
    Point<T> translated = hit + (incoming.through - incoming.origin);
    return Ray<long double>{
        Point<long double>(hit),
        reflection(mirror, translated)
    };
}

template <Coordinate T>
bool intersects(
    const Ray<T>& ray,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    assert(ray.origin != ray.through);
    assert(line.a != line.b);
    linear_parameter_detail::Parameters<T> values =
        linear_parameter_detail::parameters(
        ray.origin,
        ray.through,
        line.a,
        line.b
    );
    if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
        return on_line(line, ray.origin, eps);
    }
    return linear_parameter_detail::ratio_nonnegative<T>(
        values.first_numerator,
        values.denominator,
        eps
    );
}

template <Coordinate T>
bool intersects(
    const Line<T>& line,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return intersects(ray, line, eps);
}

template <Coordinate T>
bool intersects(
    const Ray<T>& ray,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    assert(ray.origin != ray.through);
    if (segment.a == segment.b) return on_ray(ray, segment.a, eps);

    linear_parameter_detail::Parameters<T> values =
        linear_parameter_detail::parameters(
        ray.origin,
        ray.through,
        segment.a,
        segment.b
    );
    if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
        if (orientation(ray.origin, ray.through, segment.a, eps) != 0) {
            return false;
        }
        return on_ray(ray, segment.a, eps) ||
               on_ray(ray, segment.b, eps) ||
               on_segment(segment, ray.origin, eps);
    }
    return linear_parameter_detail::ratio_nonnegative<T>(
               values.first_numerator,
               values.denominator,
               eps
           ) &&
           linear_parameter_detail::ratio_in_unit_interval<T>(
               values.second_numerator,
               values.denominator,
               eps
           );
}

template <Coordinate T>
bool intersects(
    const Segment<T>& segment,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return intersects(ray, segment, eps);
}

template <Coordinate T>
bool intersects(
    const Ray<T>& first,
    const Ray<T>& second,
    long double eps = 1e-12L
) {
    assert(first.origin != first.through);
    assert(second.origin != second.through);
    linear_parameter_detail::Parameters<T> values =
        linear_parameter_detail::parameters(
        first.origin,
        first.through,
        second.origin,
        second.through
    );
    if (linear_parameter_detail::denominator_sign(values, eps) == 0) {
        if (orientation(first.origin, first.through, second.origin, eps) != 0) {
            return false;
        }
        return on_ray(first, second.origin, eps) ||
               on_ray(second, first.origin, eps);
    }
    return linear_parameter_detail::ratio_nonnegative<T>(
               values.first_numerator,
               values.denominator,
               eps
           ) &&
           linear_parameter_detail::ratio_nonnegative<T>(
               values.second_numerator,
               values.denominator,
               eps
           );
}

namespace linear_intersection_detail {

enum class Domain {
    Line,
    Segment,
    Ray,
};

template <Coordinate T>
struct ParametricObject {
    Point<T> origin;
    Point<T> through;
    Domain domain;
};

template <Coordinate T>
ParametricObject<T> parametric_object(const Line<T>& line) {
    assert(line.a != line.b);
    return ParametricObject<T>{line.a, line.b, Domain::Line};
}

template <Coordinate T>
ParametricObject<T> parametric_object(const Segment<T>& segment) {
    return ParametricObject<T>{segment.a, segment.b, Domain::Segment};
}

template <Coordinate T>
ParametricObject<T> parametric_object(const Ray<T>& ray) {
    assert(ray.origin != ray.through);
    return ParametricObject<T>{ray.origin, ray.through, Domain::Ray};
}

template <Coordinate T>
bool contains(
    const ParametricObject<T>& object,
    const Point<T>& point,
    long double eps
) {
    if (object.domain == Domain::Line) {
        return on_line(Line<T>{object.origin, object.through}, point, eps);
    }
    if (object.domain == Domain::Segment) {
        return on_segment(
            Segment<T>{object.origin, object.through},
            point,
            eps
        );
    }
    return on_ray(Ray<T>{object.origin, object.through}, point, eps);
}

template <Coordinate T>
bool accepts_parameter(
    Domain domain,
    wide_type<T> numerator,
    wide_type<T> denominator,
    long double eps
) {
    if (domain == Domain::Line) return true;
    if (domain == Domain::Ray) {
        return linear_parameter_detail::ratio_nonnegative<T>(
            numerator,
            denominator,
            eps
        );
    }
    return linear_parameter_detail::ratio_in_unit_interval<T>(
        numerator,
        denominator,
        eps
    );
}

template <Coordinate T>
Point<long double> point_at_ratio(
    const ParametricObject<T>& object,
    wide_type<T> numerator,
    wide_type<T> denominator
) {
    const long double ratio =
        static_cast<long double>(numerator) /
        static_cast<long double>(denominator);
    const Point<long double> origin(object.origin);
    const Point<long double> direction =
        Point<long double>(object.through) - origin;
    return origin + direction * ratio;
}

template <Coordinate T>
struct AxisProjection {
    bool use_x;
    bool negate;

    wide_type<T> operator()(const Point<T>& point) const {
        const wide_type<T> value = use_x
            ? wide_type<T>(point.x)
            : wide_type<T>(point.y);
        return negate ? -value : value;
    }
};

template <Coordinate T>
AxisProjection<T> axis_projection(const ParametricObject<T>& object) {
    using W = wide_type<T>;
    const W direction_x = W(object.through.x) - W(object.origin.x);
    const W direction_y = W(object.through.y) - W(object.origin.y);
    const bool use_x =
        predicate_detail::absolute(direction_x) >=
        predicate_detail::absolute(direction_y);
    const W component = use_x ? direction_x : direction_y;
    assert(component != W(0));
    return AxisProjection<T>{use_x, component < W(0)};
}

template <Coordinate T>
struct ParameterInterval {
    bool has_lower;
    bool has_upper;
    wide_type<T> lower;
    wide_type<T> upper;
};

template <Coordinate T>
ParameterInterval<T> parameter_interval(
    const ParametricObject<T>& object,
    const AxisProjection<T>& projection
) {
    using W = wide_type<T>;
    const W origin = projection(object.origin);
    const W through = projection(object.through);
    if (object.domain == Domain::Line) {
        return ParameterInterval<T>{false, false, W(0), W(0)};
    }
    if (object.domain == Domain::Segment) {
        return ParameterInterval<T>{
            true,
            true,
            std::min(origin, through),
            std::max(origin, through),
        };
    }
    if (origin < through) {
        return ParameterInterval<T>{true, false, origin, W(0)};
    }
    return ParameterInterval<T>{false, true, W(0), origin};
}

template <Coordinate T>
ParameterInterval<T> intersect_intervals(
    ParameterInterval<T> first,
    const ParameterInterval<T>& second
) {
    if (
        second.has_lower &&
        (!first.has_lower || first.lower < second.lower)
    ) {
        first.has_lower = true;
        first.lower = second.lower;
    }
    if (
        second.has_upper &&
        (!first.has_upper || second.upper < first.upper)
    ) {
        first.has_upper = true;
        first.upper = second.upper;
    }
    return first;
}

template <Coordinate T>
Point<long double> point_at_projection(
    const ParametricObject<T>& object,
    const AxisProjection<T>& projection,
    long double target
) {
    const long double origin =
        static_cast<long double>(projection(object.origin));
    const long double through =
        static_cast<long double>(projection(object.through));
    const long double ratio = (target - origin) / (through - origin);
    const Point<long double> point(object.origin);
    const Point<long double> direction =
        Point<long double>(object.through) - point;
    return point + direction * ratio;
}

template <Coordinate T>
LinearIntersection collinear_intersection(
    const ParametricObject<T>& first,
    const ParametricObject<T>& second,
    long double eps
) {
    using W = wide_type<T>;
    const AxisProjection<T> projection = axis_projection(first);
    const ParameterInterval<T> first_interval =
        parameter_interval(first, projection);
    const ParameterInterval<T> second_interval =
        parameter_interval(second, projection);
    const ParameterInterval<T> common =
        intersect_intervals(first_interval, second_interval);

    W scale = predicate_detail::absolute(
        projection(first.through) - projection(first.origin)
    );
    scale = std::max(
        scale,
        predicate_detail::absolute(
            projection(second.through) - projection(second.origin)
        )
    );

    if (common.has_lower && common.has_upper) {
        const int order = predicate_detail::scaled_sign<ExactCoordinate<T>>(
            common.lower - common.upper,
            scale,
            eps
        );
        if (order > 0) return make_empty();
        if (order == 0) {
            const long double coordinate =
                (
                    static_cast<long double>(common.lower) +
                    static_cast<long double>(common.upper)
                ) / 2.0L;
            return make_point(
                point_at_projection(first, projection, coordinate)
            );
        }
        return make_object(
            LinearIntersectionKind::Segment,
            point_at_projection(
                first,
                projection,
                static_cast<long double>(common.lower)
            ),
            point_at_projection(
                first,
                projection,
                static_cast<long double>(common.upper)
            )
        );
    }

    const Point<long double> direction =
        Point<long double>(first.through) -
        Point<long double>(first.origin);
    if (common.has_lower) {
        const Point<long double> origin = point_at_projection(
            first,
            projection,
            static_cast<long double>(common.lower)
        );
        return make_object(
            LinearIntersectionKind::Ray,
            origin,
            origin + direction
        );
    }
    if (common.has_upper) {
        const Point<long double> origin = point_at_projection(
            first,
            projection,
            static_cast<long double>(common.upper)
        );
        return make_object(
            LinearIntersectionKind::Ray,
            origin,
            origin - direction
        );
    }
    return make_object(
        LinearIntersectionKind::Line,
        first.origin,
        first.through
    );
}

template <Coordinate T>
LinearIntersection intersect(
    const ParametricObject<T>& first,
    const ParametricObject<T>& second,
    long double eps
) {
    const bool first_degenerate = first.origin == first.through;
    const bool second_degenerate = second.origin == second.through;
    if (first_degenerate) {
        assert(first.domain == Domain::Segment);
        if (contains(second, first.origin, eps)) {
            return make_point(first.origin);
        }
        return make_empty();
    }
    if (second_degenerate) {
        assert(second.domain == Domain::Segment);
        if (contains(first, second.origin, eps)) {
            return make_point(second.origin);
        }
        return make_empty();
    }

    const linear_parameter_detail::Parameters<T> values =
        linear_parameter_detail::parameters(
        first.origin,
        first.through,
        second.origin,
        second.through
    );
    if (linear_parameter_detail::denominator_sign(values, eps) != 0) {
        if (
            !accepts_parameter<T>(
                first.domain,
                values.first_numerator,
                values.denominator,
                eps
            ) ||
            !accepts_parameter<T>(
                second.domain,
                values.second_numerator,
                values.denominator,
                eps
            )
        ) {
            return make_empty();
        }
        return make_point(
            point_at_ratio(
                first,
                values.first_numerator,
                values.denominator
            )
        );
    }
    if (
        orientation(
            first.origin,
            first.through,
            second.origin,
            eps
        ) != 0
    ) {
        return make_empty();
    }
    return collinear_intersection(first, second, eps);
}

}  // namespace linear_intersection_detail

template <Coordinate T>
LinearIntersection linear_intersection(
    const Line<T>& first,
    const Line<T>& second,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(first),
        linear_intersection_detail::parametric_object(second),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Line<T>& line,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(line),
        linear_intersection_detail::parametric_object(segment),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Segment<T>& segment,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(segment),
        linear_intersection_detail::parametric_object(line),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Segment<T>& first,
    const Segment<T>& second,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(first),
        linear_intersection_detail::parametric_object(second),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Ray<T>& ray,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(ray),
        linear_intersection_detail::parametric_object(line),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Line<T>& line,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(line),
        linear_intersection_detail::parametric_object(ray),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Ray<T>& ray,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(ray),
        linear_intersection_detail::parametric_object(segment),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Segment<T>& segment,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(segment),
        linear_intersection_detail::parametric_object(ray),
        eps
    );
}

template <Coordinate T>
LinearIntersection linear_intersection(
    const Ray<T>& first,
    const Ray<T>& second,
    long double eps = 1e-12L
) {
    return linear_intersection_detail::intersect(
        linear_intersection_detail::parametric_object(first),
        linear_intersection_detail::parametric_object(second),
        eps
    );
}

namespace closest_points_detail {

inline ClosestPoints reversed(const ClosestPoints& result) {
    return ClosestPoints{result.second, result.first};
}

inline bool point_less(
    const Point<long double>& first,
    const Point<long double>& second
) {
    if (first.x != second.x) return first.x < second.x;
    return first.y < second.y;
}

inline ClosestPoints common_point(const LinearIntersection& intersection) {
    assert(intersection.kind != LinearIntersectionKind::Empty);
    Point<long double> point = intersection.first;
    if (intersection.kind == LinearIntersectionKind::Segment) {
        if (point_less(intersection.second, point)) {
            point = intersection.second;
        }
    } else if (intersection.kind == LinearIntersectionKind::Line) {
        const Line<long double> line{
            intersection.first,
            intersection.second
        };
        point = projection(line, Point<long double>(0, 0));
    }
    return ClosestPoints{point, point};
}

inline long double separation2(const ClosestPoints& result) {
    return distance2(result.first, result.second);
}

inline bool canonical_less(
    const ClosestPoints& first,
    const ClosestPoints& second
) {
    Point<long double> first_start = first.first;
    Point<long double> first_finish = first.second;
    if (point_less(first_finish, first_start)) {
        std::swap(first_start, first_finish);
    }
    Point<long double> second_start = second.first;
    Point<long double> second_finish = second.second;
    if (point_less(second_finish, second_start)) {
        std::swap(second_start, second_finish);
    }
    if (point_less(first_start, second_start)) return true;
    if (point_less(second_start, first_start)) return false;
    return point_less(first_finish, second_finish);
}

inline void consider(ClosestPoints& best, const ClosestPoints& candidate) {
    const long double best_distance = separation2(best);
    const long double candidate_distance = separation2(candidate);
    if (
        candidate_distance < best_distance ||
        (
            candidate_distance == best_distance &&
            canonical_less(candidate, best)
        )
    ) {
        best = candidate;
    }
}

}  // namespace closest_points_detail

template <Coordinate T>
ClosestPoints closest_points(
    const Point<T>& first,
    const Point<T>& second
) {
    return ClosestPoints{
        Point<long double>(first),
        Point<long double>(second),
    };
}

template <Coordinate T>
ClosestPoints closest_points(
    const Line<T>& line,
    const Point<T>& point
) {
    return ClosestPoints{
        projection(line, point),
        Point<long double>(point),
    };
}

template <Coordinate T>
ClosestPoints closest_points(
    const Point<T>& point,
    const Line<T>& line
) {
    return closest_points_detail::reversed(closest_points(line, point));
}

template <Coordinate T>
ClosestPoints closest_points(
    const Segment<T>& segment,
    const Point<T>& point
) {
    return ClosestPoints{
        projection(segment, point),
        Point<long double>(point),
    };
}

template <Coordinate T>
ClosestPoints closest_points(
    const Point<T>& point,
    const Segment<T>& segment
) {
    return closest_points_detail::reversed(closest_points(segment, point));
}

template <Coordinate T>
ClosestPoints closest_points(
    const Ray<T>& ray,
    const Point<T>& point
) {
    return ClosestPoints{
        projection(ray, point),
        Point<long double>(point),
    };
}

template <Coordinate T>
ClosestPoints closest_points(
    const Point<T>& point,
    const Ray<T>& ray
) {
    return closest_points_detail::reversed(closest_points(ray, point));
}

template <Coordinate T>
ClosestPoints closest_points(
    const Line<T>& first,
    const Line<T>& second,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(first, second, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    ClosestPoints result = closest_points(first, second.a);
    closest_points_detail::consider(
        result,
        closest_points(first.a, second)
    );
    return result;
}

template <Coordinate T>
ClosestPoints closest_points(
    const Line<T>& line,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(line, segment, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    ClosestPoints result = closest_points(line, segment.a);
    closest_points_detail::consider(
        result,
        closest_points(line, segment.b)
    );
    return result;
}

template <Coordinate T>
ClosestPoints closest_points(
    const Segment<T>& segment,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(line, segment, eps)
    );
}

template <Coordinate T>
ClosestPoints closest_points(
    const Segment<T>& first,
    const Segment<T>& second,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(first, second, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    ClosestPoints result = closest_points(first, second.a);
    closest_points_detail::consider(
        result,
        closest_points(first, second.b)
    );
    closest_points_detail::consider(
        result,
        closest_points(first.a, second)
    );
    closest_points_detail::consider(
        result,
        closest_points(first.b, second)
    );
    return result;
}

template <Coordinate T>
ClosestPoints closest_points(
    const Line<T>& line,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(line, ray, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    return closest_points(line, ray.origin);
}

template <Coordinate T>
ClosestPoints closest_points(
    const Ray<T>& ray,
    const Line<T>& line,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(closest_points(line, ray, eps));
}

template <Coordinate T>
ClosestPoints closest_points(
    const Ray<T>& ray,
    const Segment<T>& segment,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(ray, segment, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    ClosestPoints result = closest_points(ray, segment.a);
    closest_points_detail::consider(
        result,
        closest_points(ray, segment.b)
    );
    closest_points_detail::consider(
        result,
        closest_points(ray.origin, segment)
    );
    return result;
}

template <Coordinate T>
ClosestPoints closest_points(
    const Segment<T>& segment,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(ray, segment, eps)
    );
}

template <Coordinate T>
ClosestPoints closest_points(
    const Ray<T>& first,
    const Ray<T>& second,
    long double eps = 1e-12L
) {
    const LinearIntersection intersection =
        linear_intersection(first, second, eps);
    if (intersection.kind != LinearIntersectionKind::Empty) {
        return closest_points_detail::common_point(intersection);
    }
    ClosestPoints result = closest_points(first, second.origin);
    closest_points_detail::consider(
        result,
        closest_points(first.origin, second)
    );
    return result;
}

template <Coordinate T>
long double distance(const Line<T>& line, const Point<T>& point) {
    const ClosestPoints result = closest_points(line, point);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Point<T>& point, const Line<T>& line) {
    return distance(line, point);
}

template <Coordinate T>
long double distance(const Segment<T>& segment, const Point<T>& point) {
    const ClosestPoints result = closest_points(segment, point);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Point<T>& point, const Segment<T>& segment) {
    return distance(segment, point);
}

template <Coordinate T>
long double distance(const Ray<T>& ray, const Point<T>& point) {
    const ClosestPoints result = closest_points(ray, point);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Point<T>& point, const Ray<T>& ray) {
    return distance(ray, point);
}

template <Coordinate T>
long double distance(const Line<T>& first, const Line<T>& second) {
    const ClosestPoints result = closest_points(first, second);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Line<T>& line, const Segment<T>& segment) {
    const ClosestPoints result = closest_points(line, segment);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Segment<T>& segment, const Line<T>& line) {
    return distance(line, segment);
}

template <Coordinate T>
long double distance(const Segment<T>& first, const Segment<T>& second) {
    const ClosestPoints result = closest_points(first, second);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Line<T>& line, const Ray<T>& ray) {
    const ClosestPoints result = closest_points(line, ray);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Ray<T>& ray, const Line<T>& line) {
    return distance(line, ray);
}

template <Coordinate T>
long double distance(const Ray<T>& ray, const Segment<T>& segment) {
    const ClosestPoints result = closest_points(ray, segment);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T>
long double distance(const Segment<T>& segment, const Ray<T>& ray) {
    return distance(ray, segment);
}

template <Coordinate T>
long double distance(const Ray<T>& first, const Ray<T>& second) {
    const ClosestPoints result = closest_points(first, second);
    return geometry::distance(result.first, result.second);
}

}  // namespace geometry
}  // namespace m1une


#line 17 "geometry/half_plane_intersection.hpp"

namespace m1une {
namespace geometry {

enum class HalfPlaneIntersectionStatus {
    Empty,
    Unbounded,
    Degenerate,
    Bounded,
};

struct HalfPlaneIntersectionResult {
    HalfPlaneIntersectionStatus status;
    std::vector<Point<long double>> polygon;
};

namespace half_plane_intersection_detail {

struct HalfPlane {
    Point<long double> point;
    Point<long double> direction;
    long double angle;

    HalfPlane(
        const Point<long double>& point_value,
        const Point<long double>& direction_value
    ) : point(point_value), direction(direction_value) {
        angle = std::atan2(direction.y, direction.x);
        if (angle < 0) angle += 2 * std::numbers::pi_v<long double>;
    }
};

inline bool direction_less(const HalfPlane& first, const HalfPlane& second) {
    return first.angle < second.angle;
}

inline bool parallel(
    const HalfPlane& first,
    const HalfPlane& second,
    long double eps
) {
    return std::fabs(cross(first.direction, second.direction)) <= eps;
}

inline bool same_direction(
    const HalfPlane& first,
    const HalfPlane& second,
    long double eps
) {
    return parallel(first, second, eps) &&
           dot(first.direction, second.direction) > 0;
}

inline bool outside(
    const HalfPlane& half_plane,
    const Point<long double>& point,
    long double eps
) {
    return cross(half_plane.direction, point - half_plane.point) < -eps;
}

inline bool more_restrictive(
    const HalfPlane& candidate,
    const HalfPlane& current,
    long double eps
) {
    return cross(
        current.direction,
        candidate.point - current.point
    ) > eps;
}

inline std::optional<Point<long double>> intersection(
    const HalfPlane& first,
    const HalfPlane& second,
    long double eps
) {
    long double denominator = cross(first.direction, second.direction);
    if (std::fabs(denominator) <= eps) return std::nullopt;
    long double ratio = cross(
        second.point - first.point,
        second.direction
    ) / denominator;
    return first.point + first.direction * ratio;
}

inline void merge_same_direction(
    std::vector<HalfPlane>& half_planes,
    const HalfPlane& half_plane,
    long double eps
) {
    if (
        half_planes.empty() ||
        !same_direction(half_planes.back(), half_plane, eps)
    ) {
        half_planes.push_back(half_plane);
        return;
    }
    if (more_restrictive(half_plane, half_planes.back(), eps)) {
        half_planes.back() = half_plane;
    }
}

inline void merge_cyclic_ends(
    std::vector<HalfPlane>& half_planes,
    long double eps
) {
    if (
        half_planes.size() < 2 ||
        !same_direction(half_planes.front(), half_planes.back(), eps)
    ) {
        return;
    }
    if (more_restrictive(half_planes.back(), half_planes.front(), eps)) {
        half_planes.front() = half_planes.back();
    }
    half_planes.pop_back();
}

inline bool has_feasible_point(
    std::vector<HalfPlane> half_planes,
    long double eps
) {
    std::mt19937_64 generator(0x6a09e667f3bcc909ULL);
    std::shuffle(half_planes.begin(), half_planes.end(), generator);

    Point<long double> feasible(0, 0);
    for (std::size_t index = 0; index < half_planes.size(); ++index) {
        const HalfPlane& current = half_planes[index];
        if (!outside(current, feasible, eps)) continue;

        Point<long double> normal(
            -current.direction.y,
            current.direction.x
        );
        Point<long double> base = normal * dot(normal, current.point);
        long double lower = -std::numeric_limits<long double>::infinity();
        long double upper = std::numeric_limits<long double>::infinity();
        for (std::size_t previous_index = 0;
             previous_index < index;
             ++previous_index) {
            const HalfPlane& previous = half_planes[previous_index];
            long double coefficient = cross(
                previous.direction,
                current.direction
            );
            long double constant = cross(
                previous.direction,
                base - previous.point
            );
            if (std::fabs(coefficient) <= eps) {
                if (constant < -eps) return false;
                continue;
            }

            long double bound = (-eps - constant) / coefficient;
            if (coefficient > 0) {
                lower = std::max(lower, bound);
            } else {
                upper = std::min(upper, bound);
            }
            if (lower > upper) return false;
        }

        long double parameter = 0;
        if (parameter < lower) parameter = lower;
        if (parameter > upper) parameter = upper;
        feasible = base + current.direction * parameter;
    }
    return true;
}

inline bool has_bounded_recession_cone(
    const std::vector<HalfPlane>& half_planes,
    long double eps
) {
    if (half_planes.empty()) return false;

    constexpr long double pi = std::numbers::pi_v<long double>;
    long double maximum_gap =
        half_planes.front().angle + 2 * pi - half_planes.back().angle;
    for (std::size_t index = 1; index < half_planes.size(); ++index) {
        maximum_gap = std::max(
            maximum_gap,
            half_planes[index].angle - half_planes[index - 1].angle
        );
    }
    return maximum_gap < pi - eps;
}

}  // namespace half_plane_intersection_detail

// Each directed line keeps its closed left half-plane. Returns the vertices of
// a bounded intersection with positive area in counterclockwise order. Empty,
// unbounded, and bounded zero-area intersections have distinct statuses.
template <Coordinate T>
HalfPlaneIntersectionResult half_plane_intersection(
    const std::vector<Line<T>>& half_planes,
    long double eps = 1e-12L
) {
    using half_plane_intersection_detail::HalfPlane;
    namespace detail = half_plane_intersection_detail;

    assert(eps >= 0);
    std::vector<HalfPlane> sorted;
    sorted.reserve(half_planes.size());
    for (const Line<T>& line : half_planes) {
        assert(line.a != line.b);
        Point<long double> point(line.a);
        Point<long double> direction = Point<long double>(line.b) - point;
        long double length = norm(direction);
        direction = direction / length;
        sorted.push_back(HalfPlane{point, direction});
    }
    if (!detail::has_feasible_point(sorted, eps)) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Empty,
            {},
        };
    }
    std::sort(sorted.begin(), sorted.end(), detail::direction_less);
    if (!detail::has_bounded_recession_cone(sorted, eps)) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Unbounded,
            {},
        };
    }
    if (sorted.size() < 3) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Degenerate,
            {},
        };
    }

    std::vector<HalfPlane> unique;
    unique.reserve(sorted.size());
    for (const HalfPlane& half_plane : sorted) {
        detail::merge_same_direction(unique, half_plane, eps);
    }
    detail::merge_cyclic_ends(unique, eps);
    if (unique.size() < 3) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Degenerate,
            {},
        };
    }

    std::deque<HalfPlane> deque;
    for (const HalfPlane& half_plane : unique) {
        while (deque.size() >= 2) {
            auto point = detail::intersection(
                deque[deque.size() - 2],
                deque.back(),
                eps
            );
            if (!point.has_value()) {
                return HalfPlaneIntersectionResult{
                    HalfPlaneIntersectionStatus::Degenerate,
                    {},
                };
            }
            if (!detail::outside(half_plane, *point, eps)) break;
            deque.pop_back();
        }
        while (deque.size() >= 2) {
            auto point = detail::intersection(deque[0], deque[1], eps);
            if (!point.has_value()) {
                return HalfPlaneIntersectionResult{
                    HalfPlaneIntersectionStatus::Degenerate,
                    {},
                };
            }
            if (!detail::outside(half_plane, *point, eps)) break;
            deque.pop_front();
        }
        deque.push_back(half_plane);
    }

    while (deque.size() >= 3) {
        auto point = detail::intersection(
            deque[deque.size() - 2],
            deque.back(),
            eps
        );
        if (!point.has_value()) {
            return HalfPlaneIntersectionResult{
                HalfPlaneIntersectionStatus::Degenerate,
                {},
            };
        }
        if (!detail::outside(deque.front(), *point, eps)) break;
        deque.pop_back();
    }
    while (deque.size() >= 3) {
        auto point = detail::intersection(deque[0], deque[1], eps);
        if (!point.has_value()) {
            return HalfPlaneIntersectionResult{
                HalfPlaneIntersectionStatus::Degenerate,
                {},
            };
        }
        if (!detail::outside(deque.back(), *point, eps)) break;
        deque.pop_front();
    }
    if (deque.size() < 3) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Degenerate,
            {},
        };
    }

    std::vector<Point<long double>> polygon;
    polygon.reserve(deque.size());
    for (std::size_t index = 0; index < deque.size(); ++index) {
        auto point = detail::intersection(
            deque[index],
            deque[(index + 1) % deque.size()],
            eps
        );
        if (!point.has_value()) {
            return HalfPlaneIntersectionResult{
                HalfPlaneIntersectionStatus::Degenerate,
                {},
            };
        }
        if (
            polygon.empty() ||
            distance(polygon.back(), *point) > eps
        ) {
            polygon.push_back(*point);
        }
    }
    if (
        polygon.size() >= 2 &&
        distance(polygon.front(), polygon.back()) <= eps
    ) {
        polygon.pop_back();
    }
    if (polygon.size() < 3) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Degenerate,
            {},
        };
    }

    long double signed_area2 = 0;
    Point<long double> origin = polygon.front();
    for (std::size_t index = 1; index + 1 < polygon.size(); ++index) {
        signed_area2 += cross(
            polygon[index] - origin,
            polygon[index + 1] - origin
        );
    }
    if (signed_area2 <= eps) {
        return HalfPlaneIntersectionResult{
            HalfPlaneIntersectionStatus::Degenerate,
            {},
        };
    }

    auto first = std::min_element(polygon.begin(), polygon.end());
    std::rotate(polygon.begin(), first, polygon.end());
    return HalfPlaneIntersectionResult{
        HalfPlaneIntersectionStatus::Bounded,
        std::move(polygon),
    };
}

}  // namespace geometry
}  // namespace m1une


#line 1 "geometry/minkowski_sum.hpp"



#line 8 "geometry/minkowski_sum.hpp"

#line 1 "geometry/detail/convex_polygon_normalize.hpp"



#line 8 "geometry/detail/convex_polygon_normalize.hpp"

#line 10 "geometry/detail/convex_polygon_normalize.hpp"

namespace m1une {
namespace geometry {
namespace convex_polygon_detail {

template <Coordinate T>
wide_type<T> boundary_area2(const std::vector<Point<T>>& polygon) {
    wide_type<T> result = 0;
    for (std::size_t index = 0; index < polygon.size(); ++index) {
        result += cross(
            polygon[index],
            polygon[(index + 1) % polygon.size()]
        );
    }
    return result;
}

template <Coordinate T>
std::vector<Point<T>> normalize_convex_boundary(
    std::vector<Point<T>> polygon,
    long double eps
) {
    if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
        polygon.pop_back();
    }
    polygon.erase(
        std::unique(polygon.begin(), polygon.end()),
        polygon.end()
    );
    if (polygon.size() >= 2 && polygon.front() == polygon.back()) {
        polygon.pop_back();
    }
    if (polygon.size() <= 1) return polygon;
    if (
        polygon.size() >= 3 &&
        sign<T>(boundary_area2(polygon), eps) < 0
    ) {
        std::reverse(polygon.begin(), polygon.end());
    }

    const auto start = std::min_element(
        polygon.begin(),
        polygon.end(),
        [](const Point<T>& first, const Point<T>& second) {
            if (first.y != second.y) return first.y < second.y;
            return first.x < second.x;
        }
    );
    std::rotate(polygon.begin(), start, polygon.end());

    if (polygon.size() >= 3) {
        std::vector<Point<T>> cleaned;
        const std::size_t size = polygon.size();
        cleaned.reserve(size);
        for (std::size_t index = 0; index < size; ++index) {
            const Point<T>& previous = polygon[(index + size - 1) % size];
            const Point<T>& current = polygon[index];
            const Point<T>& next = polygon[(index + 1) % size];
            if (
                orientation(previous, current, next, eps) != 0 ||
                sign<T>(dot(current - previous, next - current), eps) < 0
            ) {
                cleaned.push_back(current);
            }
        }
        polygon = std::move(cleaned);
    }
    return polygon;
}

}  // namespace convex_polygon_detail
}  // namespace geometry
}  // namespace m1une


#line 10 "geometry/minkowski_sum.hpp"

namespace m1une {
namespace geometry {

// Returns the normalized boundary of the Minkowski sum of two nonempty
// ordered convex polygons.
template <Coordinate T>
std::vector<Point<T>> minkowski_sum(
    std::vector<Point<T>> first,
    std::vector<Point<T>> second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    first = convex_polygon_detail::normalize_convex_boundary(
        std::move(first),
        eps
    );
    second = convex_polygon_detail::normalize_convex_boundary(
        std::move(second),
        eps
    );

    if (first.size() == 1 || second.size() == 1) {
        if (second.size() == 1) std::swap(first, second);
        for (Point<T>& point : second) point += first[0];
        return convex_polygon_detail::normalize_convex_boundary(
            std::move(second),
            eps
        );
    }

    std::vector<Point<T>> first_edges;
    std::vector<Point<T>> second_edges;
    first_edges.reserve(first.size());
    second_edges.reserve(second.size());
    for (std::size_t index = 0; index < first.size(); ++index) {
        first_edges.push_back(
            first[(index + 1) % first.size()] - first[index]
        );
    }
    for (std::size_t index = 0; index < second.size(); ++index) {
        second_edges.push_back(
            second[(index + 1) % second.size()] - second[index]
        );
    }

    Point<T> current = first.front() + second.front();
    std::vector<Point<T>> result;
    result.reserve(first.size() + second.size());
    result.push_back(current);
    std::size_t first_index = 0;
    std::size_t second_index = 0;
    while (
        first_index < first_edges.size() ||
        second_index < second_edges.size()
    ) {
        Point<T> step;
        if (first_index == first_edges.size()) {
            step = second_edges[second_index++];
        } else if (second_index == second_edges.size()) {
            step = first_edges[first_index++];
        } else {
            const auto turn = cross(
                first_edges[first_index],
                second_edges[second_index]
            );
            if (turn > 0) {
                step = first_edges[first_index++];
            } else if (turn < 0) {
                step = second_edges[second_index++];
            } else {
                step = first_edges[first_index++] +
                       second_edges[second_index++];
            }
        }
        current += step;
        if (
            first_index < first_edges.size() ||
            second_index < second_edges.size()
        ) {
            result.push_back(current);
        }
    }
    return convex_polygon_detail::normalize_convex_boundary(
        std::move(result),
        eps
    );
}

}  // namespace geometry
}  // namespace m1une


#line 1 "geometry/polygon.hpp"



#line 14 "geometry/polygon.hpp"

#line 1 "geometry/circle.hpp"



#line 13 "geometry/circle.hpp"

#line 15 "geometry/circle.hpp"

namespace m1une {
namespace geometry {

template <Coordinate T>
struct Circle {
    Point<T> center;
    T radius;
    bool filled = true;
};

enum class PointInCircle {
    Outside = 0,
    Boundary = 1,
    Inside = 2,
};

enum class CircleRelation {
    Separate,
    ExternallyTangent,
    Intersecting,
    InternallyTangent,
    Contained,
    Coincident,
};

enum class AngularCoverageKind {
    Empty,
    Point,
    Arc,
    Full,
};

struct AngularCoverage {
    AngularCoverageKind kind = AngularCoverageKind::Empty;
    long double begin = 0.0L;
    long double end = 0.0L;
};

struct CircleContact {
    Point<long double> point;
    long double first_argument = 0.0L;
    long double second_argument = 0.0L;
};

enum class CircleContactKind {
    Empty,
    Point,
    TwoPoints,
    Coincident,
};

struct CircleCircleIntersection {
    CircleRelation relation = CircleRelation::Separate;
    CircleContactKind contact_kind = CircleContactKind::Empty;
    std::array<CircleContact, 2> contacts;
    AngularCoverage first_inside_second;
    AngularCoverage second_inside_first;

    constexpr int contact_count() const noexcept {
        if (contact_kind == CircleContactKind::Point) return 1;
        if (contact_kind == CircleContactKind::TwoPoints) return 2;
        return 0;
    }
};

struct CircleLinearContact {
    Point<long double> point;
    long double circle_argument = 0.0L;
    long double linear_parameter = 0.0L;
};

struct CircleLinearIntersection {
    int contact_count = 0;
    std::array<CircleLinearContact, 2> contacts;
};

namespace circle_detail {

inline int compare(long double first, long double second, long double eps) {
    if (first < second - eps) return -1;
    if (first > second + eps) return 1;
    return 0;
}

inline bool close(
    const Point<long double>& first,
    const Point<long double>& second,
    long double eps
) {
    return geometry::distance(first, second) <= eps;
}

inline void push_unique(
    std::vector<Point<long double>>& points,
    const Point<long double>& point,
    long double eps
) {
    for (const Point<long double>& existing : points) {
        if (close(existing, point, eps)) return;
    }
    points.push_back(point);
}

inline bool same_line(
    const Line<long double>& first,
    const Line<long double>& second,
    long double eps
) {
    Point<long double> first_direction = first.b - first.a;
    Point<long double> second_direction = second.b - second.a;
    if (std::fabs(cross(first_direction, second_direction)) > eps) {
        return false;
    }
    return std::fabs(cross(first_direction, second.a - first.a)) <= eps;
}

inline Line<long double> tangent_line(
    const Point<long double>& contact,
    Point<long double> normal,
    long double eps
) {
    Point<long double> direction(-normal.y, normal.x);
    if (
        direction.x < -eps ||
        (std::fabs(direction.x) <= eps && direction.y < 0)
    ) {
        direction = -direction;
    }
    return Line<long double>{contact, contact + direction};
}

inline long double circular_segment_angle_term(
    long double angle,
    long double sine,
    long double cosine
) {
    if (angle >= 0.01L) return angle - sine * cosine;
    const long double squared = angle * angle;
    return angle * squared * (
        2.0L / 3.0L +
        squared * (
            -2.0L / 15.0L +
            squared * (4.0L / 315.0L - squared * 2.0L / 2835.0L)
        )
    );
}

inline long double segment_disk_signed_area(
    const Point<long double>& first,
    const Point<long double>& second,
    long double radius,
    long double eps
) {
    const Point<long double> direction = second - first;
    const long double quadratic = dot(direction, direction);
    if (quadratic == 0.0L || radius == 0.0L) return 0.0L;

    std::vector<long double> cuts = {0.0L, 1.0L};
    const long double linear = 2.0L * dot(first, direction);
    const long double constant = dot(first, first) - radius * radius;
    const long double discriminant =
        linear * linear - 4.0L * quadratic * constant;
    const long double tolerance = eps * std::max({
        1.0L,
        std::fabs(linear * linear),
        std::fabs(4.0L * quadratic * constant)
    });
    if (discriminant >= -tolerance) {
        const long double root = std::sqrt(std::max(0.0L, discriminant));
        const long double first_ratio =
            (-linear - root) / (2.0L * quadratic);
        const long double second_ratio =
            (-linear + root) / (2.0L * quadratic);
        if (eps < first_ratio && first_ratio < 1.0L - eps) {
            cuts.push_back(first_ratio);
        }
        if (eps < second_ratio && second_ratio < 1.0L - eps) {
            cuts.push_back(second_ratio);
        }
    }
    std::sort(cuts.begin(), cuts.end());
    cuts.erase(
        std::unique(
            cuts.begin(),
            cuts.end(),
            [eps](long double left, long double right) {
                return std::fabs(left - right) <= eps;
            }
        ),
        cuts.end()
    );

    long double result = 0.0L;
    for (std::size_t index = 1; index < cuts.size(); ++index) {
        const long double left = cuts[index - 1];
        const long double right = cuts[index];
        const Point<long double> a = first + direction * left;
        const Point<long double> b = first + direction * right;
        const Point<long double> middle =
            first + direction * ((left + right) / 2.0L);
        if (norm(middle) <= radius + eps) {
            result += cross(a, b) / 2.0L;
        } else {
            result +=
                radius * radius * std::atan2(cross(a, b), dot(a, b)) /
                2.0L;
        }
    }
    return result;
}

}  // namespace circle_detail

template <Coordinate T>
constexpr Point<long double> centroid(const Circle<T>& circle) {
    assert(circle.radius >= 0);
    return Point<long double>(circle.center);
}

template <Coordinate T>
constexpr long double circle_circumference(const Circle<T>& circle) {
    assert(circle.radius >= 0);
    return
        2.0L * std::numbers::pi_v<long double> *
        static_cast<long double>(circle.radius);
}

template <Coordinate T>
constexpr long double circle_area(const Circle<T>& circle) {
    assert(circle.radius >= 0);
    const long double radius = static_cast<long double>(circle.radius);
    return std::numbers::pi_v<long double> * radius * radius;
}

inline long double normalize_circle_argument(long double argument) {
    const long double full = 2.0L * std::numbers::pi_v<long double>;
    argument = std::fmod(argument, full);
    if (argument < 0.0L) argument += full;
    if (argument == full) argument = 0.0L;
    return argument;
}

template <Coordinate T>
Point<long double> circle_point_at(
    const Circle<T>& circle,
    long double argument
) {
    assert(circle.radius >= 0);
    const long double radius = static_cast<long double>(circle.radius);
    return Point<long double>(circle.center) + Point<long double>(
        radius * std::cos(argument),
        radius * std::sin(argument)
    );
}

inline long double angular_measure(const AngularCoverage& coverage) {
    if (
        coverage.kind == AngularCoverageKind::Empty ||
        coverage.kind == AngularCoverageKind::Point
    ) {
        return 0.0L;
    }
    if (coverage.kind == AngularCoverageKind::Full) {
        return 2.0L * std::numbers::pi_v<long double>;
    }
    assert(coverage.kind == AngularCoverageKind::Arc);
    assert(coverage.begin <= coverage.end);
    return coverage.end - coverage.begin;
}

template <Coordinate T>
long double circle_arc_length(
    const Circle<T>& circle,
    const AngularCoverage& coverage
) {
    assert(circle.radius >= 0);
    return static_cast<long double>(circle.radius) * angular_measure(coverage);
}

template <Coordinate C, Coordinate P>
PointInCircle point_in_circle(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(eps >= 0.0L);
    if constexpr (ExactCoordinate<C> && ExactCoordinate<P>) {
        using W = std::common_type_t<wide_type<C>, wide_type<P>>;
        const W dx = W(point.x) - W(circle.center.x);
        const W dy = W(point.y) - W(circle.center.y);
        const W radius = W(circle.radius);
        const W squared_distance = dx * dx + dy * dy;
        const W squared_radius = radius * radius;
        if (squared_distance < squared_radius) return PointInCircle::Inside;
        if (squared_distance > squared_radius) return PointInCircle::Outside;
        return PointInCircle::Boundary;
    } else {
        const int relation = circle_detail::compare(
            geometry::distance(
                Point<long double>(circle.center),
                Point<long double>(point)
            ),
            static_cast<long double>(circle.radius),
            eps
        );
        if (relation < 0) return PointInCircle::Inside;
        if (relation > 0) return PointInCircle::Outside;
        return PointInCircle::Boundary;
    }
}

template <Coordinate C, Coordinate P>
bool contains(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    const PointInCircle relation = point_in_circle(circle, point, eps);
    return circle.filled
        ? relation != PointInCircle::Outside
        : relation == PointInCircle::Boundary;
}

template <Coordinate C, Coordinate P>
bool on_circle(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(eps >= 0.0L);
    if constexpr (ExactCoordinate<C> && ExactCoordinate<P>) {
        using W = std::common_type_t<wide_type<C>, wide_type<P>>;
        const W dx = W(point.x) - W(circle.center.x);
        const W dy = W(point.y) - W(circle.center.y);
        const W radius = W(circle.radius);
        return dx * dx + dy * dy == radius * radius;
    } else {
        return circle_detail::compare(
            geometry::distance(
                Point<long double>(circle.center),
                Point<long double>(point)
            ),
            static_cast<long double>(circle.radius),
            eps
        ) == 0;
    }
}

template <Coordinate C, Coordinate P>
long double circle_argument(
    const Circle<C>& circle,
    const Point<P>& point
) {
    assert(circle.radius >= 0);
    return normalize_circle_argument(std::atan2(
        static_cast<long double>(point.y) -
            static_cast<long double>(circle.center.y),
        static_cast<long double>(point.x) -
            static_cast<long double>(circle.center.x)
    ));
}

template <Coordinate C, Coordinate P>
bool intersects(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    return contains(circle, point, eps);
}

template <Coordinate P, Coordinate C>
bool intersects(
    const Point<P>& point,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return intersects(circle, point, eps);
}

template <Coordinate A, Coordinate B>
Circle<long double> circle_from_diameter(
    const Point<A>& first,
    const Point<B>& second
) {
    Point<long double> a(first);
    Point<long double> b(second);
    Point<long double> center = (a + b) / 2.0L;
    return Circle<long double>{center, geometry::distance(a, b) / 2.0L};
}

template <Coordinate T>
std::optional<Circle<long double>> incircle(
    const Point<T>& first,
    const Point<T>& second,
    const Point<T>& third,
    long double eps = 1e-12L
) {
    assert(eps >= 0.0L);
    if (orientation(first, second, third, eps) == 0) return std::nullopt;

    long double opposite_first = geometry::distance(second, third);
    long double opposite_second = geometry::distance(third, first);
    long double opposite_third = geometry::distance(first, second);
    long double perimeter =
        opposite_first + opposite_second + opposite_third;
    Point<long double> center =
        (Point<long double>(first) * opposite_first +
         Point<long double>(second) * opposite_second +
         Point<long double>(third) * opposite_third) /
        perimeter;
    long double doubled_area = std::fabs(
        static_cast<long double>(cross(first, second, third))
    );
    return Circle<long double>{center, doubled_area / perimeter};
}

template <Coordinate T>
std::optional<Circle<long double>> circumcircle(
    const Point<T>& first,
    const Point<T>& second,
    const Point<T>& third,
    long double eps = 1e-12L
) {
    assert(eps >= 0.0L);
    if (orientation(first, second, third, eps) == 0) return std::nullopt;

    Point<long double> origin(first);
    Point<long double> u = Point<long double>(second) - origin;
    Point<long double> v = Point<long double>(third) - origin;
    long double denominator = 2.0L * cross(u, v);
    long double u_norm = norm2(u);
    long double v_norm = norm2(v);
    Point<long double> offset(
        (u_norm * v.y - v_norm * u.y) / denominator,
        (u.x * v_norm - v.x * u_norm) / denominator
    );
    Point<long double> center = origin + offset;
    return Circle<long double>{center, norm(offset)};
}

template <Coordinate A, Coordinate B>
CircleRelation circle_relation(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    assert(first.radius >= 0);
    assert(second.radius >= 0);
    assert(eps >= 0.0L);
    if constexpr (ExactCoordinate<A> && ExactCoordinate<B>) {
        using W = std::common_type_t<wide_type<A>, wide_type<B>>;
        W dx = W(second.center.x) - W(first.center.x);
        W dy = W(second.center.y) - W(first.center.y);
        W squared_distance = dx * dx + dy * dy;
        W first_radius = W(first.radius);
        W second_radius = W(second.radius);
        W sum = first_radius + second_radius;
        W difference = first_radius - second_radius;
        if (difference < 0) difference = -difference;
        if (squared_distance == 0 && difference == 0) {
            return CircleRelation::Coincident;
        }
        if (squared_distance > sum * sum) return CircleRelation::Separate;
        if (squared_distance == sum * sum) {
            return CircleRelation::ExternallyTangent;
        }
        if (squared_distance < difference * difference) {
            return CircleRelation::Contained;
        }
        if (squared_distance == difference * difference) {
            return CircleRelation::InternallyTangent;
        }
        return CircleRelation::Intersecting;
    } else {
        long double center_distance = geometry::distance(
            Point<long double>(first.center),
            Point<long double>(second.center)
        );
        long double first_radius = static_cast<long double>(first.radius);
        long double second_radius = static_cast<long double>(second.radius);
        long double sum = first_radius + second_radius;
        long double difference = std::fabs(first_radius - second_radius);
        if (
            center_distance <= eps &&
            difference <= eps
        ) {
            return CircleRelation::Coincident;
        }
        int outer = circle_detail::compare(center_distance, sum, eps);
        if (outer > 0) return CircleRelation::Separate;
        if (outer == 0) return CircleRelation::ExternallyTangent;
        int inner = circle_detail::compare(center_distance, difference, eps);
        if (inner < 0) return CircleRelation::Contained;
        if (inner == 0) return CircleRelation::InternallyTangent;
        return CircleRelation::Intersecting;
    }
}

template <Coordinate C, Coordinate L>
CircleLinearIntersection circle_boundary_intersection(
    const Circle<C>& circle,
    const Line<L>& line,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(line.a != line.b);
    assert(eps >= 0.0L);

    const Point<long double> center(circle.center);
    const Point<long double> origin(line.a);
    const Point<long double> direction =
        Point<long double>(line.b) - origin;
    const long double squared_length = dot(direction, direction);
    const long double length = std::sqrt(squared_length);
    const long double foot_parameter =
        dot(center - origin, direction) / squared_length;
    const Point<long double> foot =
        origin + direction * foot_parameter;
    const long double distance_to_line = geometry::distance(center, foot);
    const long double radius = static_cast<long double>(circle.radius);
    const int relation =
        circle_detail::compare(distance_to_line, radius, eps);

    CircleLinearIntersection result;
    if (relation > 0) return result;
    if (relation == 0) {
        result.contact_count = 1;
        result.contacts[0] = CircleLinearContact{
            foot,
            circle_argument(circle, foot),
            foot_parameter
        };
        return result;
    }

    const long double offset = std::sqrt(std::max(
        0.0L,
        radius * radius - distance_to_line * distance_to_line
    ));
    const long double parameter_offset = offset / length;
    result.contact_count = 2;
    for (int index = 0; index < 2; ++index) {
        const long double parameter = foot_parameter +
            (index == 0 ? -parameter_offset : parameter_offset);
        const Point<long double> point = origin + direction * parameter;
        result.contacts[index] = CircleLinearContact{
            point,
            circle_argument(circle, point),
            parameter
        };
    }
    return result;
}

template <Coordinate L, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
    const Line<L>& line,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return circle_boundary_intersection(circle, line, eps);
}

template <Coordinate C, Coordinate R>
CircleLinearIntersection circle_boundary_intersection(
    const Circle<C>& circle,
    const Ray<R>& ray,
    long double eps = 1e-12L
) {
    assert(ray.origin != ray.through);
    const Line<R> line{ray.origin, ray.through};
    const CircleLinearIntersection line_result =
        circle_boundary_intersection(circle, line, eps);
    CircleLinearIntersection result;
    for (int index = 0; index < line_result.contact_count; ++index) {
        CircleLinearContact contact = line_result.contacts[index];
        if (contact.linear_parameter < -eps) continue;
        if (std::fabs(contact.linear_parameter) <= eps) {
            contact.linear_parameter = 0.0L;
            contact.point = Point<long double>(ray.origin);
            contact.circle_argument = circle_argument(circle, contact.point);
        }
        result.contacts[result.contact_count++] = contact;
    }
    return result;
}

template <Coordinate R, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
    const Ray<R>& ray,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return circle_boundary_intersection(circle, ray, eps);
}

template <Coordinate C, Coordinate S>
CircleLinearIntersection circle_boundary_intersection(
    const Circle<C>& circle,
    const Segment<S>& segment,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(eps >= 0.0L);
    CircleLinearIntersection result;
    if (segment.a == segment.b) {
        if (on_circle(circle, segment.a, eps)) {
            const Point<long double> point(segment.a);
            result.contact_count = 1;
            result.contacts[0] = CircleLinearContact{
                point,
                circle_argument(circle, point),
                0.0L
            };
        }
        return result;
    }

    const Line<S> line{segment.a, segment.b};
    const CircleLinearIntersection line_result =
        circle_boundary_intersection(circle, line, eps);
    for (int index = 0; index < line_result.contact_count; ++index) {
        CircleLinearContact contact = line_result.contacts[index];
        if (
            contact.linear_parameter < -eps ||
            contact.linear_parameter > 1.0L + eps
        ) {
            continue;
        }
        if (std::fabs(contact.linear_parameter) <= eps) {
            contact.linear_parameter = 0.0L;
            contact.point = Point<long double>(segment.a);
        } else if (std::fabs(contact.linear_parameter - 1.0L) <= eps) {
            contact.linear_parameter = 1.0L;
            contact.point = Point<long double>(segment.b);
        }
        contact.circle_argument = circle_argument(circle, contact.point);
        result.contacts[result.contact_count++] = contact;
    }
    return result;
}

template <Coordinate S, Coordinate C>
CircleLinearIntersection circle_boundary_intersection(
    const Segment<S>& segment,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return circle_boundary_intersection(circle, segment, eps);
}

template <Coordinate A, Coordinate B>
CircleCircleIntersection circle_boundary_intersection(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    assert(first.radius >= 0);
    assert(second.radius >= 0);
    assert(eps >= 0.0L);
    const long double full = 2.0L * std::numbers::pi_v<long double>;
    const long double first_radius = static_cast<long double>(first.radius);
    const long double second_radius = static_cast<long double>(second.radius);
    CircleCircleIntersection result;
    result.relation = circle_relation(first, second, eps);

    auto point_coverage = [](long double argument) {
        return AngularCoverage{
            AngularCoverageKind::Point,
            argument,
            argument
        };
    };
    auto full_coverage = [full]() {
        return AngularCoverage{AngularCoverageKind::Full, 0.0L, full};
    };

    if (result.relation == CircleRelation::Coincident) {
        if (first_radius == 0.0L) {
            const Point<long double> point(first.center);
            result.contact_kind = CircleContactKind::Point;
            result.contacts[0] = CircleContact{point, 0.0L, 0.0L};
            result.first_inside_second = point_coverage(0.0L);
            result.second_inside_first = point_coverage(0.0L);
        } else {
            result.contact_kind = CircleContactKind::Coincident;
            result.first_inside_second = full_coverage();
            result.second_inside_first = full_coverage();
        }
        return result;
    }

    if (result.relation == CircleRelation::Separate) return result;
    if (result.relation == CircleRelation::Contained) {
        if (first_radius < second_radius) {
            result.first_inside_second = first_radius == 0.0L
                ? point_coverage(0.0L)
                : full_coverage();
        } else {
            result.second_inside_first = second_radius == 0.0L
                ? point_coverage(0.0L)
                : full_coverage();
        }
        return result;
    }

    const Point<long double> first_center(first.center);
    const Point<long double> second_center(second.center);
    const Point<long double> center_direction = second_center - first_center;
    const long double center_distance = norm(center_direction);
    const Point<long double> unit = center_direction / center_distance;
    const long double along =
        (first_radius * first_radius - second_radius * second_radius +
         center_distance * center_distance) /
        (2.0L * center_distance);
    const Point<long double> base = first_center + unit * along;

    if (
        result.relation == CircleRelation::ExternallyTangent ||
        result.relation == CircleRelation::InternallyTangent
    ) {
        const long double first_argument =
            circle_argument(first, base);
        const long double second_argument =
            circle_argument(second, base);
        result.contact_kind = CircleContactKind::Point;
        result.contacts[0] = CircleContact{
            base,
            first_argument,
            second_argument
        };
        result.first_inside_second = point_coverage(first_argument);
        result.second_inside_first = point_coverage(second_argument);
        if (result.relation == CircleRelation::InternallyTangent) {
            if (first_radius < second_radius && first_radius > 0.0L) {
                result.first_inside_second = full_coverage();
            } else if (
                second_radius < first_radius && second_radius > 0.0L
            ) {
                result.second_inside_first = full_coverage();
            }
        }
        return result;
    }

    assert(result.relation == CircleRelation::Intersecting);
    const long double height = std::sqrt(std::max(
        0.0L,
        first_radius * first_radius - along * along
    ));
    const Point<long double> perpendicular(-unit.y, unit.x);
    const Point<long double> first_point = base - perpendicular * height;
    const Point<long double> second_point = base + perpendicular * height;
    result.contact_kind = CircleContactKind::TwoPoints;
    result.contacts[0] = CircleContact{
        first_point,
        circle_argument(first, first_point),
        circle_argument(second, first_point)
    };
    result.contacts[1] = CircleContact{
        second_point,
        circle_argument(first, second_point),
        circle_argument(second, second_point)
    };

    const long double first_begin = result.contacts[0].first_argument;
    long double first_end = result.contacts[1].first_argument;
    if (first_end <= first_begin) first_end += full;
    result.first_inside_second = AngularCoverage{
        AngularCoverageKind::Arc,
        first_begin,
        first_end
    };

    const long double second_begin = result.contacts[1].second_argument;
    long double second_end = result.contacts[0].second_argument;
    if (second_end <= second_begin) second_end += full;
    result.second_inside_first = AngularCoverage{
        AngularCoverageKind::Arc,
        second_begin,
        second_end
    };
    return result;
}

template <Coordinate C, Coordinate L>
bool intersects(
    const Circle<C>& circle,
    const Line<L>& line,
    long double eps = 1e-12L
) {
    if (circle.filled) {
        const Line<long double> converted{
            Point<long double>(line.a),
            Point<long double>(line.b)
        };
        return contains(
            circle,
            projection(converted, Point<long double>(circle.center)),
            eps
        );
    }
    return circle_boundary_intersection(circle, line, eps).contact_count > 0;
}

template <Coordinate C, Coordinate L>
bool intersects(
    const Line<L>& line,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return intersects(circle, line, eps);
}

template <Coordinate C, Coordinate R>
bool intersects(
    const Circle<C>& circle,
    const Ray<R>& ray,
    long double eps = 1e-12L
) {
    if (circle.filled) {
        const Ray<long double> converted{
            Point<long double>(ray.origin),
            Point<long double>(ray.through)
        };
        return contains(
            circle,
            projection(converted, Point<long double>(circle.center)),
            eps
        );
    }
    return circle_boundary_intersection(circle, ray, eps).contact_count > 0;
}

template <Coordinate C, Coordinate R>
bool intersects(
    const Ray<R>& ray,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return intersects(circle, ray, eps);
}

template <Coordinate C, Coordinate S>
bool intersects(
    const Circle<C>& circle,
    const Segment<S>& segment,
    long double eps = 1e-12L
) {
    if (circle.filled) {
        const Segment<long double> converted{
            Point<long double>(segment.a),
            Point<long double>(segment.b)
        };
        return contains(
            circle,
            projection(converted, Point<long double>(circle.center)),
            eps
        );
    }
    return
        circle_boundary_intersection(circle, segment, eps).contact_count > 0;
}

template <Coordinate C, Coordinate S>
bool intersects(
    const Segment<S>& segment,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return intersects(circle, segment, eps);
}

template <Coordinate A, Coordinate B>
bool intersects(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    assert(first.radius >= 0);
    assert(second.radius >= 0);
    assert(eps >= 0.0L);
    if (first.filled && second.filled) {
        if constexpr (ExactCoordinate<A> && ExactCoordinate<B>) {
            using W = std::common_type_t<wide_type<A>, wide_type<B>>;
            const W dx = W(second.center.x) - W(first.center.x);
            const W dy = W(second.center.y) - W(first.center.y);
            const W radius = W(first.radius) + W(second.radius);
            return dx * dx + dy * dy <= radius * radius;
        } else {
            const long double center_distance = geometry::distance(
                Point<long double>(first.center),
                Point<long double>(second.center)
            );
            const long double radius_sum =
                static_cast<long double>(first.radius) +
                static_cast<long double>(second.radius);
            return circle_detail::compare(
                center_distance,
                radius_sum,
                eps
            ) <= 0;
        }
    }
    if (first.filled != second.filled) {
        const long double center_distance = geometry::distance(
            Point<long double>(first.center),
            Point<long double>(second.center)
        );
        const long double boundary_radius = first.filled
            ? static_cast<long double>(second.radius)
            : static_cast<long double>(first.radius);
        const long double filled_radius = first.filled
            ? static_cast<long double>(first.radius)
            : static_cast<long double>(second.radius);
        return circle_detail::compare(
            std::fabs(center_distance - boundary_radius),
            filled_radius,
            eps
        ) <= 0;
    }
    CircleRelation relation = circle_relation(first, second, eps);
    return
        relation == CircleRelation::ExternallyTangent ||
        relation == CircleRelation::Intersecting ||
        relation == CircleRelation::InternallyTangent ||
        relation == CircleRelation::Coincident;
}

template <Coordinate R, Coordinate H, Coordinate C>
Ray<long double> reflected_ray(
    const Ray<R>& incoming,
    const Point<H>& hit,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    assert(incoming.origin != incoming.through);
    assert(eps >= 0.0L);
    assert(static_cast<long double>(circle.radius) > eps);
    assert(
        std::fabs(
            geometry::distance(
                Point<long double>(hit),
                Point<long double>(circle.center)
            ) -
            static_cast<long double>(circle.radius)
        ) <= eps
    );

    Point<long double> hit_point(hit);
    Point<long double> normal = normalized(
        hit_point - Point<long double>(circle.center)
    );
    Point<long double> incoming_direction =
        Point<long double>(incoming.through) -
        Point<long double>(incoming.origin);
    Point<long double> outgoing_direction =
        incoming_direction - normal * (2.0L * dot(incoming_direction, normal));
    return Ray<long double>{hit_point, hit_point + outgoing_direction};
}

template <Coordinate C, Coordinate P>
std::vector<Point<long double>> tangent_points(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(eps >= 0.0L);
    Point<long double> center(circle.center);
    Point<long double> external(point);
    Point<long double> direction = external - center;
    long double squared_distance = dot(direction, direction);
    long double radius = static_cast<long double>(circle.radius);
    if (radius == 0.0L) return {center};

    long double center_distance = std::sqrt(squared_distance);
    int relation = circle_detail::compare(center_distance, radius, eps);
    if (relation < 0) return {};
    if (relation == 0) {
        return {center + direction * (radius / center_distance)};
    }

    Point<long double> base =
        center + direction * (radius * radius / squared_distance);
    long double scale =
        radius * std::sqrt(std::max(
            0.0L,
            squared_distance - radius * radius
        )) /
        squared_distance;
    Point<long double> perpendicular(-direction.y, direction.x);
    Point<long double> first = base - perpendicular * scale;
    Point<long double> second = base + perpendicular * scale;
    if (second < first) std::swap(first, second);
    return {first, second};
}

template <Coordinate A, Coordinate B>
std::vector<Line<long double>> common_tangents(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    assert(first.radius >= 0);
    assert(second.radius >= 0);
    assert(eps >= 0.0L);
    Point<long double> first_center(first.center);
    Point<long double> second_center(second.center);
    Point<long double> direction = second_center - first_center;
    long double squared_distance = dot(direction, direction);
    long double center_distance = std::sqrt(squared_distance);
    if (center_distance <= eps) return {};

    long double first_radius = static_cast<long double>(first.radius);
    long double second_radius = static_cast<long double>(second.radius);
    std::vector<Line<long double>> result;
    for (int second_side : {1, -1}) {
        long double difference =
            first_radius - second_side * second_radius;
        int relation = circle_detail::compare(
            std::fabs(difference),
            center_distance,
            eps
        );
        if (relation > 0) continue;
        long double perpendicular_length = relation == 0 ? 0.0L : std::sqrt(
            std::max(0.0L, squared_distance - difference * difference)
        );
        int choices = perpendicular_length <= eps ? 1 : 2;
        for (int choice = 0; choice < choices; ++choice) {
            long double side = choice == 0 ? -1.0L : 1.0L;
            Point<long double> normal =
                direction * (difference / squared_distance) +
                Point<long double>(-direction.y, direction.x) *
                    (side * perpendicular_length / squared_distance);
            normal = normalized(normal);
            Point<long double> contact =
                first_center + normal * first_radius;
            Line<long double> tangent =
                circle_detail::tangent_line(contact, normal, eps);
            bool duplicate = false;
            for (const Line<long double>& existing : result) {
                if (circle_detail::same_line(existing, tangent, eps)) {
                    duplicate = true;
                    break;
                }
            }
            if (!duplicate) result.push_back(tangent);
        }
    }
    std::sort(
        result.begin(),
        result.end(),
        [](const Line<long double>& left, const Line<long double>& right) {
            if (left.a != right.a) return left.a < right.a;
            return left.b < right.b;
        }
    );
    return result;
}

template <Coordinate A, Coordinate B>
std::vector<Point<long double>> common_tangent_points(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    std::vector<Point<long double>> result;
    for (const Line<long double>& line : common_tangents(first, second, eps)) {
        circle_detail::push_unique(result, line.a, eps);
    }
    std::sort(result.begin(), result.end());
    return result;
}

// These area functions use the enclosed disks, independent of `filled`.
template <Coordinate A, Coordinate B>
long double circle_circle_intersection_area(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    assert(first.radius >= 0);
    assert(second.radius >= 0);
    assert(eps >= 0.0L);
    const long double first_radius = static_cast<long double>(first.radius);
    const long double second_radius = static_cast<long double>(second.radius);
    const CircleRelation relation = circle_relation(first, second, eps);
    if (
        relation == CircleRelation::Separate ||
        relation == CircleRelation::ExternallyTangent
    ) {
        return 0.0L;
    }
    if (
        relation == CircleRelation::Contained ||
        relation == CircleRelation::InternallyTangent ||
        relation == CircleRelation::Coincident
    ) {
        const long double radius = std::min(first_radius, second_radius);
        return std::numbers::pi_v<long double> * radius * radius;
    }

    const long double center_distance = geometry::distance(
        Point<long double>(first.center),
        Point<long double>(second.center)
    );
    const long double first_cosine = std::clamp(
        (
            (center_distance - second_radius) *
                (center_distance + second_radius) +
            first_radius * first_radius
        ) / (2.0L * center_distance * first_radius),
        -1.0L,
        1.0L
    );
    const long double second_cosine = std::clamp(
        (
            (center_distance - first_radius) *
                (center_distance + first_radius) +
            second_radius * second_radius
        ) / (2.0L * center_distance * second_radius),
        -1.0L,
        1.0L
    );
    const long double radicand =
        (-center_distance + first_radius + second_radius) *
        (center_distance + first_radius - second_radius) *
        (center_distance - first_radius + second_radius) *
        (center_distance + first_radius + second_radius);
    const long double height =
        std::sqrt(std::max(0.0L, radicand)) / (2.0L * center_distance);
    const long double first_sine =
        std::clamp(height / first_radius, 0.0L, 1.0L);
    const long double second_sine =
        std::clamp(height / second_radius, 0.0L, 1.0L);
    const long double first_angle = std::atan2(first_sine, first_cosine);
    const long double second_angle = std::atan2(second_sine, second_cosine);
    return
        first_radius * first_radius *
            circle_detail::circular_segment_angle_term(
                first_angle,
                first_sine,
                first_cosine
            ) +
        second_radius * second_radius *
            circle_detail::circular_segment_angle_term(
                second_angle,
                second_sine,
                second_cosine
            );
}

template <Coordinate C, Coordinate P>
long double circle_polygon_intersection_area(
    const Circle<C>& circle,
    const std::vector<Point<P>>& polygon,
    long double eps = 1e-12L
) {
    assert(circle.radius >= 0);
    assert(eps >= 0.0L);
    if (polygon.empty() || circle.radius == 0) return 0.0L;

    const Point<long double> center(circle.center);
    const long double radius = static_cast<long double>(circle.radius);
    long double result = 0.0L;
    for (std::size_t index = 0; index < polygon.size(); ++index) {
        const Point<long double> first =
            Point<long double>(polygon[index]) - center;
        const Point<long double> second =
            Point<long double>(polygon[(index + 1) % polygon.size()]) - center;
        result += circle_detail::segment_disk_signed_area(
            first,
            second,
            radius,
            eps
        );
    }
    return std::fabs(result);
}

namespace circle_detail {

template <Coordinate T>
Point<long double> point_toward(
    const Circle<T>& circle,
    const Point<long double>& target
) {
    assert(circle.radius >= 0);
    const Point<long double> center(circle.center);
    const long double radius = static_cast<long double>(circle.radius);
    const Point<long double> direction = target - center;
    const long double length = norm(direction);
    if (length == 0.0L) {
        return center + Point<long double>(-radius, 0.0L);
    }
    return center + direction * (radius / length);
}

inline void consider(
    ClosestPoints& best,
    const ClosestPoints& candidate
) {
    closest_points_detail::consider(best, candidate);
}

}  // namespace circle_detail

template <Coordinate C, Coordinate P>
ClosestPoints closest_points(
    const Circle<C>& circle,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    const Point<long double> converted(point);
    if (circle.filled && contains(circle, converted, eps)) {
        return ClosestPoints{converted, converted};
    }
    return ClosestPoints{
        circle_detail::point_toward(circle, converted),
        converted
    };
}

template <Coordinate P, Coordinate C>
ClosestPoints closest_points(
    const Point<P>& point,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(circle, point, eps)
    );
}

template <Coordinate C, Coordinate L>
ClosestPoints closest_points(
    const Circle<C>& circle,
    const Line<L>& line,
    long double eps = 1e-12L
) {
    const Line<long double> converted{
        Point<long double>(line.a),
        Point<long double>(line.b)
    };
    const Point<long double> point = projection(
        converted,
        Point<long double>(circle.center)
    );
    if (circle.filled) return closest_points(circle, point, eps);

    const CircleLinearIntersection common =
        circle_boundary_intersection(circle, line, eps);
    if (common.contact_count > 0) {
        return ClosestPoints{
            common.contacts[0].point,
            common.contacts[0].point
        };
    }
    return ClosestPoints{
        circle_detail::point_toward(circle, point),
        point
    };
}

template <Coordinate L, Coordinate C>
ClosestPoints closest_points(
    const Line<L>& line,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(closest_points(circle, line, eps));
}

template <Coordinate C, Coordinate R>
ClosestPoints closest_points(
    const Circle<C>& circle,
    const Ray<R>& ray,
    long double eps = 1e-12L
) {
    const Ray<long double> converted{
        Point<long double>(ray.origin),
        Point<long double>(ray.through)
    };
    const Point<long double> point = projection(
        converted,
        Point<long double>(circle.center)
    );
    if (circle.filled) return closest_points(circle, point, eps);

    const CircleLinearIntersection common =
        circle_boundary_intersection(circle, ray, eps);
    if (common.contact_count > 0) {
        return ClosestPoints{
            common.contacts[0].point,
            common.contacts[0].point
        };
    }
    return ClosestPoints{
        circle_detail::point_toward(circle, point),
        point
    };
}

template <Coordinate R, Coordinate C>
ClosestPoints closest_points(
    const Ray<R>& ray,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(closest_points(circle, ray, eps));
}

template <Coordinate C, Coordinate S>
ClosestPoints closest_points(
    const Circle<C>& circle,
    const Segment<S>& segment,
    long double eps = 1e-12L
) {
    const Segment<long double> converted{
        Point<long double>(segment.a),
        Point<long double>(segment.b)
    };
    const Point<long double> center(circle.center);
    const Point<long double> projected = projection(converted, center);
    if (circle.filled) return closest_points(circle, projected, eps);

    const CircleLinearIntersection common =
        circle_boundary_intersection(circle, segment, eps);
    if (common.contact_count > 0) {
        return ClosestPoints{
            common.contacts[0].point,
            common.contacts[0].point
        };
    }
    ClosestPoints result{
        circle_detail::point_toward(circle, projected),
        projected
    };
    for (const Point<long double>& point : {converted.a, converted.b}) {
        circle_detail::consider(
            result,
            ClosestPoints{circle_detail::point_toward(circle, point), point}
        );
    }
    return result;
}

template <Coordinate S, Coordinate C>
ClosestPoints closest_points(
    const Segment<S>& segment,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(circle, segment, eps)
    );
}

template <Coordinate A, Coordinate B>
ClosestPoints closest_points(
    const Circle<A>& first,
    const Circle<B>& second,
    long double eps = 1e-12L
) {
    if (first.filled && !second.filled) {
        return closest_points_detail::reversed(
            closest_points(second, first, eps)
        );
    }

    if (!first.filled && second.filled) {
        const ClosestPoints center_result =
            closest_points(first, second.center, eps);
        if (contains(second, center_result.first, eps)) {
            return ClosestPoints{center_result.first, center_result.first};
        }
        const ClosestPoints filled_result =
            closest_points(second, center_result.first, eps);
        return ClosestPoints{center_result.first, filled_result.first};
    }

    if (first.filled && second.filled) {
        assert(first.radius >= 0);
        assert(second.radius >= 0);
        const Point<long double> first_center(first.center);
        const Point<long double> second_center(second.center);
        Point<long double> direction = second_center - first_center;
        const long double center_distance = norm(direction);
        if (center_distance == 0.0L) {
            return ClosestPoints{first_center, first_center};
        }
        direction = direction / center_distance;
        const long double first_radius =
            static_cast<long double>(first.radius);
        const long double second_radius =
            static_cast<long double>(second.radius);
        if (intersects(first, second, eps)) {
            const long double left = std::max(
                -first_radius,
                center_distance - second_radius
            );
            const long double right = std::min(
                first_radius,
                center_distance + second_radius
            );
            const Point<long double> common =
                first_center + direction * ((left + right) / 2.0L);
            return ClosestPoints{common, common};
        }
        return ClosestPoints{
            first_center + direction * first_radius,
            second_center - direction * second_radius
        };
    }

    const CircleCircleIntersection common =
        circle_boundary_intersection(first, second, eps);
    if (common.contact_count() > 0) {
        return ClosestPoints{
            common.contacts[0].point,
            common.contacts[0].point
        };
    }

    const Point<long double> first_center(first.center);
    const Point<long double> second_center(second.center);
    if (circle_relation(first, second, eps) == CircleRelation::Coincident) {
        const Point<long double> point =
            circle_detail::point_toward(first, first_center);
        return ClosestPoints{point, point};
    }
    Point<long double> direction = second_center - first_center;
    const long double center_distance = norm(direction);
    if (center_distance == 0.0L) {
        const Point<long double> first_point =
            circle_detail::point_toward(first, first_center);
        const Point<long double> second_point =
            circle_detail::point_toward(second, second_center);
        return ClosestPoints{first_point, second_point};
    }
    direction = direction / center_distance;
    const long double first_radius = static_cast<long double>(first.radius);
    const long double second_radius = static_cast<long double>(second.radius);
    ClosestPoints result{
        first_center + direction * first_radius,
        second_center + direction * second_radius
    };
    for (const long double first_sign : {-1.0L, 1.0L}) {
        for (const long double second_sign : {-1.0L, 1.0L}) {
            circle_detail::consider(
                result,
                ClosestPoints{
                    first_center + direction * (first_sign * first_radius),
                    second_center + direction * (second_sign * second_radius)
                }
            );
        }
    }
    return result;
}

template <Coordinate C, Coordinate P>
long double distance(const Circle<C>& circle, const Point<P>& point) {
    const ClosestPoints result = closest_points(circle, point);
    return geometry::distance(result.first, result.second);
}

template <Coordinate P, Coordinate C>
long double distance(const Point<P>& point, const Circle<C>& circle) {
    return distance(circle, point);
}

template <Coordinate C, Coordinate L>
long double distance(const Circle<C>& circle, const Line<L>& line) {
    const ClosestPoints result = closest_points(circle, line);
    return geometry::distance(result.first, result.second);
}

template <Coordinate L, Coordinate C>
long double distance(const Line<L>& line, const Circle<C>& circle) {
    return distance(circle, line);
}

template <Coordinate C, Coordinate R>
long double distance(const Circle<C>& circle, const Ray<R>& ray) {
    const ClosestPoints result = closest_points(circle, ray);
    return geometry::distance(result.first, result.second);
}

template <Coordinate R, Coordinate C>
long double distance(const Ray<R>& ray, const Circle<C>& circle) {
    return distance(circle, ray);
}

template <Coordinate C, Coordinate S>
long double distance(const Circle<C>& circle, const Segment<S>& segment) {
    const ClosestPoints result = closest_points(circle, segment);
    return geometry::distance(result.first, result.second);
}

template <Coordinate S, Coordinate C>
long double distance(const Segment<S>& segment, const Circle<C>& circle) {
    return distance(circle, segment);
}

template <Coordinate A, Coordinate B>
long double distance(const Circle<A>& first, const Circle<B>& second) {
    const ClosestPoints result = closest_points(first, second);
    return geometry::distance(result.first, result.second);
}

}  // namespace geometry
}  // namespace m1une


#line 16 "geometry/polygon.hpp"

namespace m1une {
namespace geometry {

enum class PointInPolygon {
    Outside = 0,
    Boundary = 1,
    Inside = 2,
};

template <Coordinate T>
struct Polygon {
    std::vector<Point<T>> vertices;
    bool filled = true;
};

template <Coordinate T, typename Scalar>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
Polygon<std::common_type_t<T, Scalar>> operator*(
    const Polygon<T>& polygon,
    Scalar scalar
) {
    using Result = std::common_type_t<T, Scalar>;
    Polygon<Result> scaled;
    scaled.vertices.reserve(polygon.vertices.size());
    for (const Point<T>& point : polygon.vertices) {
        scaled.vertices.push_back(point * scalar);
    }
    scaled.filled = polygon.filled;
    return scaled;
}

template <typename Scalar, Coordinate T>
requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
Polygon<std::common_type_t<T, Scalar>> operator*(
    Scalar scalar,
    const Polygon<T>& polygon
) {
    return polygon * scalar;
}

struct ParameterInterval {
    long double begin = 0.0L;
    long double end = 0.0L;
};

template <Coordinate T>
constexpr Point<long double> centroid(
    const std::array<Point<T>, 3>& triangle
) {
    return Point<long double>(
        (
            static_cast<long double>(triangle[0].x) +
            static_cast<long double>(triangle[1].x) +
            static_cast<long double>(triangle[2].x)
        ) / 3,
        (
            static_cast<long double>(triangle[0].y) +
            static_cast<long double>(triangle[1].y) +
            static_cast<long double>(triangle[2].y)
        ) / 3
    );
}

namespace polygon_detail {

template <Coordinate T>
std::vector<Point<T>> clean_polygon_vertices(
    std::vector<Point<T>> polygon,
    long double eps
) {
    if (
        polygon.size() >= 2 &&
        polygon.front() == polygon.back()
    ) {
        polygon.pop_back();
    }

    std::vector<Point<T>> deduplicated;
    for (const Point<T>& point : polygon) {
        if (deduplicated.empty() || deduplicated.back() != point) {
            deduplicated.push_back(point);
        }
    }
    if (
        deduplicated.size() >= 2 &&
        deduplicated.front() == deduplicated.back()
    ) {
        deduplicated.pop_back();
    }

    bool changed = true;
    while (changed && deduplicated.size() >= 3) {
        changed = false;
        std::vector<Point<T>> cleaned;
        std::size_t size = deduplicated.size();
        for (std::size_t index = 0; index < size; ++index) {
            const Point<T>& previous =
                deduplicated[(index + size - 1) % size];
            const Point<T>& current = deduplicated[index];
            const Point<T>& next =
                deduplicated[(index + 1) % size];
            if (
                orientation(previous, current, next, eps) == 0 &&
                sign<T>(dot(current - previous, next - current), eps) >= 0
            ) {
                changed = true;
            } else {
                cleaned.push_back(current);
            }
        }
        deduplicated = std::move(cleaned);
    }
    return deduplicated;
}

template <Coordinate T>
bool in_ccw_triangle(
    const Point<T>& point,
    const Point<T>& first,
    const Point<T>& second,
    const Point<T>& third,
    long double eps
) {
    return
        orientation(first, second, point, eps) >= 0 &&
        orientation(second, third, point, eps) >= 0 &&
        orientation(third, first, point, eps) >= 0;
}

}  // namespace polygon_detail

template <Coordinate T>
wide_type<T> polygon_area2(const std::vector<Point<T>>& polygon) {
    wide_type<T> result = 0;
    std::size_t n = polygon.size();
    for (std::size_t i = 0; i < n; i++) {
        result += cross(polygon[i], polygon[(i + 1) % n]);
    }
    return result;
}

template <Coordinate T>
long double polygon_area(const std::vector<Point<T>>& polygon) {
    return std::fabs(static_cast<long double>(polygon_area2(polygon))) / 2;
}

template <Coordinate T>
std::optional<Point<long double>> polygon_centroid(
    const std::vector<Point<T>>& polygon,
    long double eps = 1e-12L
) {
    if (polygon.size() < 3) return std::nullopt;

    wide_type<T> signed_area2 = polygon_area2(polygon);
    if (sign<T>(signed_area2, eps) == 0) return std::nullopt;

    long double x_numerator = 0;
    long double y_numerator = 0;
    std::size_t size = polygon.size();
    for (std::size_t index = 0; index < size; ++index) {
        const Point<T>& current = polygon[index];
        const Point<T>& next = polygon[(index + 1) % size];
        long double weight = static_cast<long double>(cross(current, next));
        x_numerator +=
            (static_cast<long double>(current.x) +
             static_cast<long double>(next.x)) *
            weight;
        y_numerator +=
            (static_cast<long double>(current.y) +
             static_cast<long double>(next.y)) *
            weight;
    }
    long double denominator =
        3.0L * static_cast<long double>(signed_area2);
    return Point<long double>(
        x_numerator / denominator,
        y_numerator / denominator
    );
}

template <Coordinate T>
std::optional<Point<long double>> centroid(
    const std::vector<Point<T>>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon, eps);
}

template <Coordinate T>
std::optional<Point<long double>> polygon_center_of_gravity(
    const std::vector<Point<T>>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon, eps);
}

template <Coordinate T>
bool is_simple_polygon(
    const std::vector<Point<T>>& polygon,
    long double eps = 1e-12L
) {
    if (polygon.size() < 3) return false;
    std::size_t size = polygon.size();
    for (std::size_t index = 0; index < size; ++index) {
        const Point<T>& previous = polygon[(index + size - 1) % size];
        const Point<T>& current = polygon[index];
        const Point<T>& next = polygon[(index + 1) % size];
        if (current == next) return false;
        if (
            orientation(previous, current, next, eps) == 0 &&
            sign<T>(dot(current - previous, next - current), eps) < 0
        ) {
            return false;
        }
    }
    for (std::size_t first_index = 0; first_index < size; ++first_index) {
        Segment<T> first{
            polygon[first_index],
            polygon[(first_index + 1) % size]
        };
        for (
            std::size_t second_index = first_index + 1;
            second_index < size;
            ++second_index
        ) {
            bool adjacent =
                second_index == first_index + 1 ||
                (first_index == 0 && second_index + 1 == size);
            if (adjacent) continue;

            Segment<T> second{
                polygon[second_index],
                polygon[(second_index + 1) % size]
            };
            if (intersects(first, second, eps)) return false;
        }
    }
    return true;
}

template <Coordinate T>
std::optional<std::vector<std::array<Point<T>, 3>>> triangulate_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    polygon =
        polygon_detail::clean_polygon_vertices(std::move(polygon), eps);
    if (polygon.size() < 3) return std::nullopt;

    wide_type<T> signed_area2 = polygon_area2(polygon);
    if (sign<T>(signed_area2, eps) == 0) return std::nullopt;
    if (!is_simple_polygon(polygon, eps)) return std::nullopt;
    if (sign<T>(signed_area2, eps) < 0) {
        std::reverse(polygon.begin(), polygon.end());
    }

    std::vector<std::size_t> remaining(polygon.size());
    for (std::size_t index = 0; index < polygon.size(); ++index) {
        remaining[index] = index;
    }

    std::vector<std::array<Point<T>, 3>> result;
    result.reserve(polygon.size() - 2);
    while (remaining.size() > 3) {
        bool found_ear = false;
        std::size_t size = remaining.size();
        for (std::size_t position = 0; position < size; ++position) {
            std::size_t previous_index =
                remaining[(position + size - 1) % size];
            std::size_t current_index = remaining[position];
            std::size_t next_index =
                remaining[(position + 1) % size];
            const Point<T>& previous = polygon[previous_index];
            const Point<T>& current = polygon[current_index];
            const Point<T>& next = polygon[next_index];
            if (orientation(previous, current, next, eps) <= 0) continue;

            bool contains_vertex = false;
            for (std::size_t other_index : remaining) {
                if (
                    other_index == previous_index ||
                    other_index == current_index ||
                    other_index == next_index
                ) {
                    continue;
                }
                if (
                    polygon_detail::in_ccw_triangle(
                        polygon[other_index],
                        previous,
                        current,
                        next,
                        eps
                    )
                ) {
                    contains_vertex = true;
                    break;
                }
            }
            if (contains_vertex) continue;

            std::array<Point<T>, 3> triangle;
            triangle[0] = previous;
            triangle[1] = current;
            triangle[2] = next;
            result.push_back(std::move(triangle));
            remaining.erase(
                remaining.begin() +
                static_cast<std::ptrdiff_t>(position)
            );
            found_ear = true;
            break;
        }
        if (!found_ear) return std::nullopt;
    }

    std::array<Point<T>, 3> triangle;
    triangle[0] = polygon[remaining[0]];
    triangle[1] = polygon[remaining[1]];
    triangle[2] = polygon[remaining[2]];
    if (orientation(triangle[0], triangle[1], triangle[2], eps) <= 0) {
        return std::nullopt;
    }
    result.push_back(std::move(triangle));
    return result;
}

template <Coordinate T>
PointInPolygon point_in_polygon(
    const std::vector<Point<T>>& polygon,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    bool inside = false;
    std::size_t n = polygon.size();
    for (std::size_t i = 0; i < n; i++) {
        const Point<T>& a = polygon[i];
        const Point<T>& b = polygon[(i + 1) % n];
        if (on_segment(Segment<T>{a, b}, point, eps)) {
            return PointInPolygon::Boundary;
        }

        if (a.y <= point.y) {
            if (point.y < b.y && orientation(a, b, point, eps) > 0) {
                inside = !inside;
            }
        } else if (b.y <= point.y && orientation(a, b, point, eps) < 0) {
            inside = !inside;
        }
    }
    return inside ? PointInPolygon::Inside : PointInPolygon::Outside;
}

template <Coordinate T, Coordinate P>
PointInPolygon point_in_polygon(
    const Polygon<T>& polygon,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    assert(polygon.vertices.size() >= 3);
    if constexpr (std::is_same_v<T, P>) {
        return point_in_polygon(polygon.vertices, point, eps);
    } else {
        std::vector<Point<long double>> vertices;
        vertices.reserve(polygon.vertices.size());
        for (const Point<T>& vertex : polygon.vertices) {
            vertices.emplace_back(vertex);
        }
        return point_in_polygon(vertices, Point<long double>(point), eps);
    }
}

template <Coordinate T, Coordinate P>
bool contains(
    const Polygon<T>& polygon,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    const PointInPolygon relation = point_in_polygon(polygon, point, eps);
    return polygon.filled
        ? relation != PointInPolygon::Outside
        : relation == PointInPolygon::Boundary;
}

namespace polygon_clip_detail {

struct Event {
    long double parameter;
    bool toggle;
};

inline bool close_parameter(
    long double first,
    long double second,
    long double eps
) {
    return std::fabs(first - second) <= eps * std::max({
        1.0L,
        std::fabs(first),
        std::fabs(second)
    });
}

inline long double parameter_on_line(
    const Point<long double>& origin,
    const Point<long double>& direction,
    const Point<long double>& point
) {
    return dot(point - origin, direction) / dot(direction, direction);
}

inline std::vector<Event> grouped_events(
    std::vector<Event> events,
    long double eps
) {
    std::sort(
        events.begin(),
        events.end(),
        [](const Event& first, const Event& second) {
            return first.parameter < second.parameter;
        }
    );
    std::vector<Event> result;
    for (const Event& event : events) {
        if (
            result.empty() ||
            !close_parameter(result.back().parameter, event.parameter, eps)
        ) {
            result.push_back(event);
        } else {
            result.back().toggle = result.back().toggle != event.toggle;
        }
    }
    return result;
}

inline std::vector<ParameterInterval> merge_intervals(
    std::vector<ParameterInterval> intervals,
    long double eps
) {
    for (ParameterInterval& interval : intervals) {
        if (interval.end < interval.begin) {
            std::swap(interval.begin, interval.end);
        }
    }
    std::sort(
        intervals.begin(),
        intervals.end(),
        [](const ParameterInterval& first, const ParameterInterval& second) {
            if (first.begin != second.begin) return first.begin < second.begin;
            return first.end < second.end;
        }
    );
    std::vector<ParameterInterval> result;
    for (const ParameterInterval& interval : intervals) {
        if (
            result.empty() ||
            interval.begin > result.back().end + eps * std::max({
                1.0L,
                std::fabs(interval.begin),
                std::fabs(result.back().end)
            })
        ) {
            result.push_back(interval);
        } else if (result.back().end < interval.end) {
            result.back().end = interval.end;
        }
    }
    return result;
}

template <Coordinate T>
std::vector<ParameterInterval> clip_line(
    const Point<long double>& origin,
    const Point<long double>& direction,
    const Polygon<T>& polygon,
    long double eps
) {
    assert(polygon.vertices.size() >= 3);
    assert(direction != Point<long double>());
    std::vector<Event> events;
    std::vector<ParameterInterval> intervals;
    events.reserve(polygon.vertices.size() * 2);
    intervals.reserve(polygon.vertices.size());

    const Line<long double> line{origin, origin + direction};
    for (std::size_t index = 0; index < polygon.vertices.size(); ++index) {
        const Point<long double> first(polygon.vertices[index]);
        const Point<long double> second(
            polygon.vertices[(index + 1) % polygon.vertices.size()]
        );
        assert(first != second);
        const Segment<long double> edge{first, second};
        const LinearIntersection intersection =
            linear_intersection(line, edge, eps);
        if (intersection.kind == LinearIntersectionKind::Empty) continue;
        if (intersection.kind == LinearIntersectionKind::Segment) {
            intervals.push_back(ParameterInterval{
                parameter_on_line(origin, direction, intersection.first),
                parameter_on_line(origin, direction, intersection.second)
            });
            continue;
        }
        assert(intersection.kind == LinearIntersectionKind::Point);
        const Point<long double> point = intersection.first;
        const Point<long double> edge_direction = second - first;
        const long double edge_parameter =
            dot(point - first, edge_direction) /
            dot(edge_direction, edge_direction);
        bool toggle = false;
        if (edge_parameter <= eps) {
            toggle = orientation(line.a, line.b, second, eps) > 0;
        } else if (edge_parameter >= 1.0L - eps) {
            toggle = orientation(line.a, line.b, first, eps) > 0;
        } else {
            toggle = true;
        }
        events.push_back(Event{
            parameter_on_line(origin, direction, point),
            toggle
        });
    }

    const std::vector<Event> grouped = grouped_events(std::move(events), eps);
    if (polygon.filled) {
        bool inside = false;
        for (std::size_t index = 0; index < grouped.size(); ++index) {
            if (index > 0 && inside) {
                intervals.push_back(ParameterInterval{
                    grouped[index - 1].parameter,
                    grouped[index].parameter
                });
            }
            intervals.push_back(ParameterInterval{
                grouped[index].parameter,
                grouped[index].parameter
            });
            inside = inside != grouped[index].toggle;
        }
    } else {
        for (const Event& event : grouped) {
            intervals.push_back(ParameterInterval{
                event.parameter,
                event.parameter
            });
        }
    }
    return merge_intervals(std::move(intervals), eps);
}

inline std::vector<ParameterInterval> restrict_domain(
    const std::vector<ParameterInterval>& intervals,
    long double lower,
    long double upper,
    long double eps
) {
    std::vector<ParameterInterval> result;
    result.reserve(intervals.size());
    for (const ParameterInterval& interval : intervals) {
        long double begin = std::max(interval.begin, lower);
        long double end = std::min(interval.end, upper);
        if (
            end < begin &&
            !close_parameter(begin, end, eps)
        ) {
            continue;
        }
        if (end < begin) {
            const long double middle = (begin + end) / 2.0L;
            begin = middle;
            end = middle;
        }
        result.push_back(ParameterInterval{begin, end});
    }
    return merge_intervals(std::move(result), eps);
}

inline std::vector<Event> grouped_circle_events(
    std::vector<Event> events,
    long double eps
) {
    std::vector<Event> result = grouped_events(std::move(events), eps);
    const long double full = 2.0L * std::numbers::pi_v<long double>;
    if (
        result.size() >= 2 &&
        close_parameter(result.front().parameter + full,
                        result.back().parameter, eps)
    ) {
        result.front().parameter = 0.0L;
        result.front().toggle =
            result.front().toggle != result.back().toggle;
        result.pop_back();
    }
    return result;
}

template <Coordinate C, Coordinate T>
std::vector<AngularCoverage> clip_circle(
    const Circle<C>& circle,
    const Polygon<T>& polygon,
    long double eps
) {
    assert(circle.radius >= 0);
    assert(polygon.vertices.size() >= 3);
    if (circle.radius == 0) {
        return contains(polygon, circle.center, eps)
            ? std::vector<AngularCoverage>{AngularCoverage{
                  AngularCoverageKind::Point,
                  0.0L,
                  0.0L
              }}
            : std::vector<AngularCoverage>();
    }

    std::vector<Event> events;
    events.reserve(polygon.vertices.size() * 2);
    for (std::size_t index = 0; index < polygon.vertices.size(); ++index) {
        const Point<long double> first(polygon.vertices[index]);
        const Point<long double> second(
            polygon.vertices[(index + 1) % polygon.vertices.size()]
        );
        assert(first != second);
        const Point<long double> direction = second - first;
        const Segment<long double> edge{first, second};
        const CircleLinearIntersection intersection =
            circle_boundary_intersection(circle, edge, eps);
        for (
            int contact_index = 0;
            contact_index < intersection.contact_count;
            ++contact_index
        ) {
            const CircleLinearContact& contact =
                intersection.contacts[contact_index];
            const Point<long double> radial =
                contact.point - Point<long double>(circle.center);
            bool toggle = false;
            if (contact.linear_parameter <= eps) {
                toggle = predicate_detail::dot_sign<false>(
                    radial.x,
                    radial.y,
                    direction.x,
                    direction.y,
                    eps
                ) < 0;
            } else if (contact.linear_parameter >= 1.0L - eps) {
                toggle = predicate_detail::dot_sign<false>(
                    radial.x,
                    radial.y,
                    -direction.x,
                    -direction.y,
                    eps
                ) < 0;
            } else {
                toggle = predicate_detail::dot_sign<false>(
                    radial.x,
                    radial.y,
                    direction.x,
                    direction.y,
                    eps
                ) != 0;
            }
            events.push_back(Event{contact.circle_argument, toggle});
        }
    }

    const std::vector<Event> grouped =
        grouped_circle_events(std::move(events), eps);
    const long double full = 2.0L * std::numbers::pi_v<long double>;
    if (grouped.empty()) {
        return contains(polygon, circle_point_at(circle, 0.0L), eps)
            ? std::vector<AngularCoverage>{AngularCoverage{
                  AngularCoverageKind::Full,
                  0.0L,
                  full
              }}
            : std::vector<AngularCoverage>();
    }

    std::vector<bool> inside_gap(grouped.size(), false);
    if (polygon.filled) {
        const long double wrap_middle = normalize_circle_argument(
            (grouped.back().parameter + grouped.front().parameter + full) /
            2.0L
        );
        bool inside = point_in_polygon(
            polygon,
            circle_point_at(circle, wrap_middle),
            eps
        ) != PointInPolygon::Outside;
        for (std::size_t index = 0; index < grouped.size(); ++index) {
            inside = inside != grouped[index].toggle;
            inside_gap[index] = inside;
        }
    }

    if (
        polygon.filled &&
        std::all_of(
            inside_gap.begin(),
            inside_gap.end(),
            [](bool inside) { return inside; }
        )
    ) {
        return {AngularCoverage{AngularCoverageKind::Full, 0.0L, full}};
    }

    std::vector<AngularCoverage> result;
    for (std::size_t index = 0; index < grouped.size(); ++index) {
        const std::size_t previous =
            (index + grouped.size() - 1) % grouped.size();
        if (!inside_gap[previous] && !inside_gap[index]) {
            result.push_back(AngularCoverage{
                AngularCoverageKind::Point,
                grouped[index].parameter,
                grouped[index].parameter
            });
        }
        if (!inside_gap[previous] && inside_gap[index]) {
            std::size_t finish = index;
            while (inside_gap[finish]) {
                finish = (finish + 1) % grouped.size();
            }
            long double end = grouped[finish].parameter;
            if (end <= grouped[index].parameter) end += full;
            result.push_back(AngularCoverage{
                AngularCoverageKind::Arc,
                grouped[index].parameter,
                end
            });
        }
    }
    std::sort(
        result.begin(),
        result.end(),
        [](const AngularCoverage& first, const AngularCoverage& second) {
            return first.begin < second.begin;
        }
    );
    return result;
}

}  // namespace polygon_clip_detail

template <Coordinate L, Coordinate T>
std::vector<ParameterInterval> clip(
    const Line<L>& line,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    assert(line.a != line.b);
    assert(eps >= 0.0L);
    const Point<long double> origin(line.a);
    const Point<long double> direction =
        Point<long double>(line.b) - origin;
    return polygon_clip_detail::clip_line(
        origin,
        direction,
        polygon,
        eps
    );
}

template <Coordinate R, Coordinate T>
std::vector<ParameterInterval> clip(
    const Ray<R>& ray,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    assert(ray.origin != ray.through);
    assert(eps >= 0.0L);
    const Point<long double> origin(ray.origin);
    const Point<long double> direction =
        Point<long double>(ray.through) - origin;
    return polygon_clip_detail::restrict_domain(
        polygon_clip_detail::clip_line(origin, direction, polygon, eps),
        0.0L,
        std::numeric_limits<long double>::infinity(),
        eps
    );
}

template <Coordinate S, Coordinate T>
std::vector<ParameterInterval> clip(
    const Segment<S>& segment,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    assert(eps >= 0.0L);
    if (segment.a == segment.b) {
        return contains(polygon, segment.a, eps)
            ? std::vector<ParameterInterval>{ParameterInterval{0.0L, 0.0L}}
            : std::vector<ParameterInterval>();
    }
    const Point<long double> origin(segment.a);
    const Point<long double> direction =
        Point<long double>(segment.b) - origin;
    return polygon_clip_detail::restrict_domain(
        polygon_clip_detail::clip_line(origin, direction, polygon, eps),
        0.0L,
        1.0L,
        eps
    );
}

template <Coordinate C, Coordinate T>
std::vector<AngularCoverage> clip(
    const Circle<C>& circle,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    assert(eps >= 0.0L);
    return polygon_clip_detail::clip_circle(circle, polygon, eps);
}

template <Coordinate T>
bool intersects(
    const Ray<T>& ray,
    const std::vector<Point<T>>& polygon,
    long double eps = 1e-12L
) {
    assert(polygon.size() >= 3);
    if (point_in_polygon(polygon, ray.origin, eps) != PointInPolygon::Outside) {
        return true;
    }
    Polygon<T> region{polygon};
    return !clip(ray, region, eps).empty();
}

template <Coordinate T>
bool intersects(
    const std::vector<Point<T>>& polygon,
    const Ray<T>& ray,
    long double eps = 1e-12L
) {
    return intersects(ray, polygon, eps);
}

template <Coordinate T>
long double distance(
    const Ray<T>& ray,
    const std::vector<Point<T>>& polygon
) {
    assert(polygon.size() >= 3);
    if (intersects(ray, polygon)) return 0;
    long double result = std::numeric_limits<long double>::infinity();
    std::size_t size = polygon.size();
    for (std::size_t index = 0; index < size; ++index) {
        result = std::min(
            result,
            distance(
                ray,
                Segment<T>{
                    polygon[index],
                    polygon[(index + 1) % size]
                }
            )
        );
    }
    return result;
}

template <Coordinate T>
long double distance(
    const std::vector<Point<T>>& polygon,
    const Ray<T>& ray
) {
    return distance(ray, polygon);
}

template <Coordinate T>
bool intersects(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(first.size() >= 3);
    assert(second.size() >= 3);
    std::size_t first_size = first.size();
    std::size_t second_size = second.size();
    for (
        std::size_t first_index = 0;
        first_index < first_size;
        ++first_index
    ) {
        Segment<T> first_edge{
            first[first_index],
            first[(first_index + 1) % first_size]
        };
        for (
            std::size_t second_index = 0;
            second_index < second_size;
            ++second_index
        ) {
            Segment<T> second_edge{
                second[second_index],
                second[(second_index + 1) % second_size]
            };
            if (intersects(first_edge, second_edge, eps)) return true;
        }
    }
    return
        point_in_polygon(first, second.front(), eps) !=
            PointInPolygon::Outside ||
        point_in_polygon(second, first.front(), eps) !=
            PointInPolygon::Outside;
}

template <Coordinate T>
long double distance(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second
) {
    assert(first.size() >= 3);
    assert(second.size() >= 3);
    if (intersects(first, second)) return 0;

    long double result = std::numeric_limits<long double>::infinity();
    std::size_t first_size = first.size();
    std::size_t second_size = second.size();
    for (
        std::size_t first_index = 0;
        first_index < first_size;
        ++first_index
    ) {
        Segment<T> first_edge{
            first[first_index],
            first[(first_index + 1) % first_size]
        };
        for (
            std::size_t second_index = 0;
            second_index < second_size;
            ++second_index
        ) {
            Segment<T> second_edge{
                second[second_index],
                second[(second_index + 1) % second_size]
            };
            result = std::min(result, distance(first_edge, second_edge));
        }
    }
    return result;
}

template <Coordinate T>
wide_type<T> polygon_area2(const Polygon<T>& polygon) {
    return polygon_area2(polygon.vertices);
}

template <Coordinate T>
long double polygon_area(const Polygon<T>& polygon) {
    return polygon_area(polygon.vertices);
}

template <Coordinate T>
std::optional<Point<long double>> polygon_centroid(
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon.vertices, eps);
}

template <Coordinate T>
std::optional<Point<long double>> centroid(
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon.vertices, eps);
}

namespace polygon_detail {

template <Coordinate T>
Segment<long double> edge(const Polygon<T>& polygon, std::size_t index) {
    return Segment<long double>{
        Point<long double>(polygon.vertices[index]),
        Point<long double>(
            polygon.vertices[(index + 1) % polygon.vertices.size()]
        )
    };
}

template <Coordinate T>
ClosestPoints closest_boundary_point(
    const Polygon<T>& polygon,
    const Point<long double>& point
) {
    assert(polygon.vertices.size() >= 3);
    ClosestPoints result = closest_points(edge(polygon, 0), point);
    for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
        closest_points_detail::consider(
            result,
            closest_points(edge(polygon, index), point)
        );
    }
    return result;
}

template <Coordinate T, class Object>
ClosestPoints closest_boundary_object(
    const Polygon<T>& polygon,
    const Object& object
) {
    assert(polygon.vertices.size() >= 3);
    ClosestPoints result = closest_points(edge(polygon, 0), object);
    for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
        closest_points_detail::consider(
            result,
            closest_points(edge(polygon, index), object)
        );
    }
    return result;
}

template <Coordinate A, Coordinate B>
ClosestPoints closest_boundaries(
    const Polygon<A>& first,
    const Polygon<B>& second
) {
    assert(first.vertices.size() >= 3);
    assert(second.vertices.size() >= 3);
    ClosestPoints result = closest_points(edge(first, 0), edge(second, 0));
    for (
        std::size_t first_index = 0;
        first_index < first.vertices.size();
        ++first_index
    ) {
        for (
            std::size_t second_index = 0;
            second_index < second.vertices.size();
            ++second_index
        ) {
            closest_points_detail::consider(
                result,
                closest_points(
                    edge(first, first_index),
                    edge(second, second_index)
                )
            );
        }
    }
    return result;
}

}  // namespace polygon_detail

template <Coordinate T, Coordinate P>
ClosestPoints closest_points(
    const Polygon<T>& polygon,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    assert(polygon.vertices.size() >= 3);
    const Point<long double> converted(point);
    if (polygon.filled && contains(polygon, point, eps)) {
        return ClosestPoints{converted, converted};
    }
    return polygon_detail::closest_boundary_point(polygon, converted);
}

template <Coordinate P, Coordinate T>
ClosestPoints closest_points(
    const Point<P>& point,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(polygon, point, eps)
    );
}

template <Coordinate T, Coordinate S>
ClosestPoints closest_points(
    const Polygon<T>& polygon,
    const Segment<S>& segment,
    long double eps = 1e-12L
) {
    assert(polygon.vertices.size() >= 3);
    const Segment<long double> converted{
        Point<long double>(segment.a),
        Point<long double>(segment.b)
    };
    if (polygon.filled) {
        if (contains(polygon, segment.a, eps)) {
            const Point<long double> point(segment.a);
            return ClosestPoints{point, point};
        }
        if (contains(polygon, segment.b, eps)) {
            const Point<long double> point(segment.b);
            return ClosestPoints{point, point};
        }
    }
    return polygon_detail::closest_boundary_object(polygon, converted);
}

template <Coordinate S, Coordinate T>
ClosestPoints closest_points(
    const Segment<S>& segment,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(polygon, segment, eps)
    );
}

template <Coordinate T, Coordinate R>
ClosestPoints closest_points(
    const Polygon<T>& polygon,
    const Ray<R>& ray,
    long double eps = 1e-12L
) {
    assert(polygon.vertices.size() >= 3);
    const Ray<long double> converted{
        Point<long double>(ray.origin),
        Point<long double>(ray.through)
    };
    if (polygon.filled && contains(polygon, ray.origin, eps)) {
        const Point<long double> point(ray.origin);
        return ClosestPoints{point, point};
    }
    return polygon_detail::closest_boundary_object(polygon, converted);
}

template <Coordinate R, Coordinate T>
ClosestPoints closest_points(
    const Ray<R>& ray,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(polygon, ray, eps)
    );
}

template <Coordinate A, Coordinate B>
ClosestPoints closest_points(
    const Polygon<A>& first,
    const Polygon<B>& second,
    long double eps = 1e-12L
) {
    assert(first.vertices.size() >= 3);
    assert(second.vertices.size() >= 3);
    ClosestPoints result = polygon_detail::closest_boundaries(first, second);
    if (geometry::distance(result.first, result.second) <= eps) return result;

    if (first.filled) {
        for (const Point<B>& vertex : second.vertices) {
            if (contains(first, vertex, eps)) {
                const Point<long double> point(vertex);
                return ClosestPoints{point, point};
            }
        }
    }
    if (second.filled) {
        for (const Point<A>& vertex : first.vertices) {
            if (contains(second, vertex, eps)) {
                const Point<long double> point(vertex);
                return ClosestPoints{point, point};
            }
        }
    }
    return result;
}

template <Coordinate C, Coordinate T>
ClosestPoints closest_points(
    const Circle<C>& circle,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    assert(polygon.vertices.size() >= 3);
    ClosestPoints result = closest_points(
        circle,
        polygon_detail::edge(polygon, 0),
        eps
    );
    for (std::size_t index = 1; index < polygon.vertices.size(); ++index) {
        closest_points_detail::consider(
            result,
            closest_points(circle, polygon_detail::edge(polygon, index), eps)
        );
    }
    if (geometry::distance(result.first, result.second) <= eps) return result;

    if (polygon.filled) {
        Point<long double> member(circle.center);
        if (!circle.filled) {
            member = circle_detail::point_toward(circle, member);
        }
        if (contains(polygon, member, eps)) {
            return ClosestPoints{member, member};
        }
    }
    return result;
}

template <Coordinate T, Coordinate C>
ClosestPoints closest_points(
    const Polygon<T>& polygon,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return closest_points_detail::reversed(
        closest_points(circle, polygon, eps)
    );
}

template <Coordinate T, Coordinate P>
bool intersects(
    const Polygon<T>& polygon,
    const Point<P>& point,
    long double eps = 1e-12L
) {
    return contains(polygon, point, eps);
}

template <Coordinate P, Coordinate T>
bool intersects(
    const Point<P>& point,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return intersects(polygon, point, eps);
}

template <Coordinate T, Coordinate S>
bool intersects(
    const Polygon<T>& polygon,
    const Segment<S>& segment,
    long double eps = 1e-12L
) {
    const ClosestPoints result = closest_points(polygon, segment, eps);
    return geometry::distance(result.first, result.second) <= eps;
}

template <Coordinate S, Coordinate T>
bool intersects(
    const Segment<S>& segment,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return intersects(polygon, segment, eps);
}

template <Coordinate T, Coordinate R>
bool intersects(
    const Polygon<T>& polygon,
    const Ray<R>& ray,
    long double eps = 1e-12L
) {
    const ClosestPoints result = closest_points(polygon, ray, eps);
    return geometry::distance(result.first, result.second) <= eps;
}

template <Coordinate R, Coordinate T>
bool intersects(
    const Ray<R>& ray,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return intersects(polygon, ray, eps);
}

template <Coordinate A, Coordinate B>
bool intersects(
    const Polygon<A>& first,
    const Polygon<B>& second,
    long double eps = 1e-12L
) {
    const ClosestPoints result = closest_points(first, second, eps);
    return geometry::distance(result.first, result.second) <= eps;
}

template <Coordinate C, Coordinate T>
bool intersects(
    const Circle<C>& circle,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    const ClosestPoints result = closest_points(circle, polygon, eps);
    return geometry::distance(result.first, result.second) <= eps;
}

template <Coordinate T, Coordinate C>
bool intersects(
    const Polygon<T>& polygon,
    const Circle<C>& circle,
    long double eps = 1e-12L
) {
    return intersects(circle, polygon, eps);
}

template <Coordinate A, Coordinate B>
long double distance(
    const Polygon<A>& first,
    const Polygon<B>& second
) {
    const ClosestPoints result = closest_points(first, second);
    return geometry::distance(result.first, result.second);
}

template <Coordinate C, Coordinate T>
long double distance(
    const Circle<C>& circle,
    const Polygon<T>& polygon
) {
    const ClosestPoints result = closest_points(circle, polygon);
    return geometry::distance(result.first, result.second);
}

template <Coordinate T, Coordinate C>
long double distance(
    const Polygon<T>& polygon,
    const Circle<C>& circle
) {
    return distance(circle, polygon);
}

template <Coordinate C, Coordinate T>
long double circle_polygon_intersection_area(
    const Circle<C>& circle,
    const Polygon<T>& polygon,
    long double eps = 1e-12L
) {
    return circle_polygon_intersection_area(circle, polygon.vertices, eps);
}

template <Coordinate T, Coordinate P>
long double distance(
    const Polygon<T>& polygon,
    const Point<P>& point
) {
    const ClosestPoints result = closest_points(polygon, point);
    return geometry::distance(result.first, result.second);
}

template <Coordinate P, Coordinate T>
long double distance(
    const Point<P>& point,
    const Polygon<T>& polygon
) {
    return distance(polygon, point);
}

template <Coordinate T, Coordinate S>
long double distance(
    const Polygon<T>& polygon,
    const Segment<S>& segment
) {
    const ClosestPoints result = closest_points(polygon, segment);
    return geometry::distance(result.first, result.second);
}

template <Coordinate S, Coordinate T>
long double distance(
    const Segment<S>& segment,
    const Polygon<T>& polygon
) {
    return distance(polygon, segment);
}

template <Coordinate T, Coordinate R>
long double distance(
    const Polygon<T>& polygon,
    const Ray<R>& ray
) {
    const ClosestPoints result = closest_points(polygon, ray);
    return geometry::distance(result.first, result.second);
}

template <Coordinate R, Coordinate T>
long double distance(
    const Ray<R>& ray,
    const Polygon<T>& polygon
) {
    return distance(polygon, ray);
}

}  // namespace geometry
}  // namespace m1une


#line 22 "geometry/convex_polygon.hpp"

namespace m1une {
namespace geometry {

namespace convex_polygon_detail {

inline bool points_close(
    const Point<long double>& first,
    const Point<long double>& second,
    long double eps
) {
    return geometry::distance(first, second) <= eps;
}

inline std::vector<Point<long double>> clean_polygon(
    std::vector<Point<long double>> polygon,
    long double eps
) {
    if (polygon.empty()) return polygon;

    std::vector<Point<long double>> deduplicated;
    for (const Point<long double>& point : polygon) {
        if (
            deduplicated.empty() ||
            !points_close(deduplicated.back(), point, eps)
        ) {
            deduplicated.push_back(point);
        }
    }
    if (
        deduplicated.size() >= 2 &&
        points_close(
            deduplicated.front(),
            deduplicated.back(),
            eps
        )
    ) {
        deduplicated.pop_back();
    }
    if (deduplicated.size() <= 2) return deduplicated;
    std::vector<Point<long double>> cleaned;
    const std::size_t size = deduplicated.size();
    cleaned.reserve(size);
    for (std::size_t index = 0; index < size; ++index) {
        const Point<long double>& previous =
            deduplicated[(index + size - 1) % size];
        const Point<long double>& current = deduplicated[index];
        const Point<long double>& next =
            deduplicated[(index + 1) % size];
        if (
            orientation(previous, current, next, eps) != 0 ||
            dot(current - previous, next - current) < -eps
        ) {
            cleaned.push_back(current);
        }
    }
    return cleaned;
}

}  // namespace convex_polygon_detail

template <Coordinate T>
bool is_convex_polygon(
    const std::vector<Point<T>>& polygon,
    bool strict = false,
    long double eps = 1e-12L
) {
    std::size_t size = polygon.size();
    if (size >= 2 && polygon.front() == polygon.back()) size--;
    if (size < 3) return false;

    int direction = 0;
    for (std::size_t index = 0; index < size; ++index) {
        const Point<T>& current = polygon[index];
        const Point<T>& next = polygon[(index + 1) % size];
        const Point<T>& after = polygon[(index + 2) % size];
        if (current == next) return false;
        const int turn = orientation(current, next, after, eps);
        if (turn == 0) {
            if (strict) return false;
            continue;
        }
        if (direction != 0 && direction != turn) return false;
        direction = turn;
    }
    return !strict || direction != 0;
}

template <Coordinate T>
std::vector<Point<T>> normalize_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    return convex_polygon_detail::normalize_convex_boundary(
        std::move(polygon),
        eps
    );
}

template <Coordinate T>
PointInPolygon point_in_convex_polygon(
    const std::vector<Point<T>>& polygon,
    const Point<T>& point,
    long double eps = 1e-12L
) {
    const std::size_t size = polygon.size();
    if (size == 0) return PointInPolygon::Outside;
    if (size == 1) {
        return distance(polygon[0], point) <= eps
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }
    if (size == 2) {
        return on_segment(Segment<T>{polygon[0], polygon[1]}, point, eps)
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }

    const int order = orientation(
        polygon[0],
        polygon[1],
        polygon[size - 1],
        eps
    );
    if (order == 0) return point_in_polygon(polygon, point, eps);
    auto vertex = [&](std::size_t index) -> const Point<T>& {
        if (order > 0 || index == 0) return polygon[index];
        return polygon[size - index];
    };

    const int first_side = orientation(vertex(0), vertex(1), point, eps);
    const int last_side =
        orientation(vertex(0), vertex(size - 1), point, eps);
    if (first_side < 0 || last_side > 0) {
        return PointInPolygon::Outside;
    }
    if (first_side == 0) {
        return on_segment(Segment<T>{vertex(0), vertex(1)}, point, eps)
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }
    if (last_side == 0) {
        return on_segment(
            Segment<T>{vertex(0), vertex(size - 1)},
            point,
            eps
        )
            ? PointInPolygon::Boundary
            : PointInPolygon::Outside;
    }

    std::size_t left = 1;
    std::size_t right = size - 1;
    while (right - left >= 2) {
        const std::size_t middle = (left + right) / 2;
        if (orientation(vertex(0), vertex(middle), point, eps) >= 0) {
            left = middle;
        } else {
            right = middle;
        }
    }
    const int triangle_side =
        orientation(vertex(left), vertex(right), point, eps);
    if (triangle_side < 0) return PointInPolygon::Outside;
    if (triangle_side == 0) return PointInPolygon::Boundary;
    return PointInPolygon::Inside;
}

template <Coordinate T>
class ConvexPolygon {
   public:
    using Wide = wide_type<T>;

   private:
    std::vector<Point<T>> points;
    std::vector<Wide> area_prefix;
    long double epsilon;

    template <class Compare>
    int periodic_best(Compare better) const {
        const int size = int(points.size());
        int left = 0;
        int middle = size;
        int right = 2 * size;
        while (right - left > 2) {
            const int left_middle = (left + middle) / 2;
            const int right_middle = (middle + right + 1) / 2;
            if (better(left_middle % size, middle % size)) {
                right = middle;
                middle = left_middle;
            } else if (better(right_middle % size, middle % size)) {
                left = middle;
                middle = right_middle;
            } else {
                left = left_middle;
                right = right_middle;
            }
        }
        return middle % size;
    }

    int previous(int index) const {
        return index == 0 ? int(points.size()) - 1 : index - 1;
    }

    int next(int index) const {
        return index + 1 == int(points.size()) ? 0 : index + 1;
    }

   public:
    explicit ConvexPolygon(
        std::vector<Point<T>> polygon,
        long double eps = 1e-12L
    )
        : points(normalize_convex_polygon(std::move(polygon), eps)),
          epsilon(eps) {
        assert(
            points.size() <=
            static_cast<std::size_t>(
                std::numeric_limits<int>::max() / 2
            )
        );
        assert(
            points.size() < 3 ||
            is_convex_polygon(points, true, epsilon)
        );
        area_prefix.resize(2 * points.size() + 1, Wide(0));
        for (std::size_t index = 0; index < 2 * points.size(); ++index) {
            area_prefix[index + 1] =
                area_prefix[index] +
                cross(
                    points[index % points.size()],
                    points[(index + 1) % points.size()]
                );
        }
    }

    int size() const noexcept {
        return int(points.size());
    }

    bool empty() const noexcept {
        return points.empty();
    }

    const std::vector<Point<T>>& vertices() const noexcept {
        return points;
    }

    const Point<T>& operator[](int index) const {
        assert(0 <= index && index < size());
        return points[index];
    }

    ConvexPolygon operator+(const ConvexPolygon& other) const {
        const long double eps = std::max(epsilon, other.epsilon);
        return ConvexPolygon(
            minkowski_sum(points, other.points, eps),
            eps
        );
    }

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    ConvexPolygon<std::common_type_t<T, Scalar>> operator*(Scalar scalar) const {
        using Result = std::common_type_t<T, Scalar>;
        std::vector<Point<Result>> scaled;
        scaled.reserve(points.size());
        for (const Point<T>& point : points) {
            scaled.push_back(point * scalar);
        }
        return ConvexPolygon<Result>(std::move(scaled), epsilon);
    }

    template <typename Scalar>
    requires (std::is_arithmetic_v<Scalar> || Coordinate<Scalar>)
    friend ConvexPolygon<std::common_type_t<T, Scalar>> operator*(
        Scalar scalar,
        const ConvexPolygon& polygon
    ) {
        return polygon * scalar;
    }

    Wide area2() const {
        if (points.empty()) return Wide(0);
        return area_prefix[points.size()];
    }

    Wide chain_area2(int first, int last) const {
        assert(0 <= first && first < size());
        assert(0 <= last && last < size());
        int extended_last = last;
        if (extended_last < first) extended_last += size();
        return
            area_prefix[extended_last] - area_prefix[first] +
            cross(points[last], points[first]);
    }

    PointInPolygon contains(const Point<T>& point) const {
        return point_in_convex_polygon(points, point, epsilon);
    }

    std::pair<Wide, int> min_dot(const Point<T>& direction) const {
        assert(!points.empty());
        const int index = periodic_best([&](int first, int second) {
            return dot(points[first], direction) <
                   dot(points[second], direction);
        });
        return std::pair<Wide, int>(dot(points[index], direction), index);
    }

    std::pair<Wide, int> max_dot(const Point<T>& direction) const {
        assert(!points.empty());
        const int index = periodic_best([&](int first, int second) {
            return dot(points[first], direction) >
                   dot(points[second], direction);
        });
        return std::pair<Wide, int>(dot(points[index], direction), index);
    }

    std::pair<int, int> tangent_vertices(const Point<T>& point) const {
        assert(points.size() >= 3);
        assert(contains(point) == PointInPolygon::Outside);
        int first = periodic_best([&](int left, int right) {
            return orientation(point, points[left], points[right], epsilon) < 0;
        });
        int second = periodic_best([&](int left, int right) {
            return orientation(point, points[left], points[right], epsilon) > 0;
        });
        if (
            orientation(
                point,
                points[first],
                points[previous(first)],
                epsilon
            ) == 0
        ) {
            first = previous(first);
        }
        if (
            orientation(
                point,
                points[second],
                points[next(second)],
                epsilon
            ) == 0
        ) {
            second = next(second);
        }
        return std::pair<int, int>(first, second);
    }
};

template <Coordinate T>
std::optional<Point<long double>> centroid(
    const ConvexPolygon<T>& polygon,
    long double eps = 1e-12L
) {
    return polygon_centroid(polygon.vertices(), eps);
}

namespace convex_polygon_detail {

template <Coordinate T>
class MinkowskiDifferenceView {
   private:
    struct Cycle {
        const ConvexPolygon<T>* polygon;
        int start;
        bool negate;

        int edge_count() const {
            return polygon->size() >= 2 ? polygon->size() : 0;
        }

        Point<T> point(int index) const {
            const int size = polygon->size();
            const Point<T>& result = (*polygon)[(start + index) % size];
            return negate ? -result : result;
        }

        Point<T> edge(int index) const {
            return point((index + 1) % polygon->size()) - point(index);
        }
    };

    Cycle first;
    Cycle second;

    std::pair<int, int> prefixes(int rank) const {
        const int first_size = first.edge_count();
        const int second_size = second.edge_count();
        if (first_size + second_size == 0) {
            return std::pair<int, int>(0, 0);
        }

        int low = std::max(0, rank - second_size);
        int high = std::min(rank, first_size);
        while (low <= high) {
            const int first_prefix = (low + high) / 2;
            const int second_prefix = rank - first_prefix;
            if (
                first_prefix > 0 &&
                second_prefix < second_size &&
                entry_less(
                    second.edge(second_prefix),
                    1,
                    first.edge(first_prefix - 1),
                    0
                )
            ) {
                high = first_prefix - 1;
                continue;
            }
            if (
                second_prefix > 0 &&
                first_prefix < first_size &&
                entry_less(
                    first.edge(first_prefix),
                    0,
                    second.edge(second_prefix - 1),
                    1
                )
            ) {
                low = first_prefix + 1;
                continue;
            }
            return std::pair<int, int>(first_prefix, second_prefix);
        }
        assert(false);
        return std::pair<int, int>(0, 0);
    }

    static int direction_half(const Point<T>& direction) {
        return
            direction.y > 0 ||
            (direction.y == 0 && direction.x >= 0)
            ? 0
            : 1;
    }

    static bool entry_less(
        const Point<T>& left,
        int left_cycle,
        const Point<T>& right,
        int right_cycle
    ) {
        if constexpr (std::floating_point<T>) {
            long double left_angle = std::atan2(
                static_cast<long double>(left.y),
                static_cast<long double>(left.x)
            );
            long double right_angle = std::atan2(
                static_cast<long double>(right.y),
                static_cast<long double>(right.x)
            );
            if (left_angle < 0) {
                left_angle += 2 * std::numbers::pi_v<long double>;
            }
            if (right_angle < 0) {
                right_angle += 2 * std::numbers::pi_v<long double>;
            }
            if (left_angle != right_angle) return left_angle < right_angle;
            return left_cycle < right_cycle;
        }
        const int left_half = direction_half(left);
        const int right_half = direction_half(right);
        if (left_half != right_half) return left_half < right_half;
        const auto turn = cross(left, right);
        if (turn != 0) return turn > 0;
        return left_cycle < right_cycle;
    }

    static int negated_start(const ConvexPolygon<T>& polygon) {
        if (polygon.size() <= 1) return 0;
        int result = polygon.max_dot(Point<T>(0, 1)).second;
        const int previous = result == 0 ? polygon.size() - 1 : result - 1;
        const int next = result + 1 == polygon.size() ? 0 : result + 1;
        for (const int candidate : {previous, next}) {
            if (
                polygon[candidate].y == polygon[result].y &&
                polygon[candidate].x > polygon[result].x
            ) {
                result = candidate;
            }
        }
        return result;
    }

   public:
    MinkowskiDifferenceView(
        const ConvexPolygon<T>& minuend,
        const ConvexPolygon<T>& subtrahend
    )
        : first{&minuend, 0, false},
          second{&subtrahend, negated_start(subtrahend), true} {
        assert(!minuend.empty());
        assert(!subtrahend.empty());
    }

    int size() const {
        const int edge_count =
            first.edge_count() + second.edge_count();
        return edge_count == 0 ? 1 : edge_count;
    }

    Point<T> operator[](int rank) const {
        assert(0 <= rank && rank < size());
        const auto [first_prefix, second_prefix] = prefixes(rank);
        return
            first.point(first_prefix % first.polygon->size()) +
            second.point(second_prefix % second.polygon->size());
    }

    std::pair<Point<T>, Point<T>> components(int rank) const {
        assert(0 <= rank && rank < size());
        const auto [first_prefix, second_prefix] = prefixes(rank);
        return std::pair<Point<T>, Point<T>>(
            first.point(first_prefix % first.polygon->size()),
            -second.point(second_prefix % second.polygon->size())
        );
    }
};

struct OriginLocation {
    PointInPolygon location;
    int outside_edge;
    std::array<int, 3> simplex;
    int simplex_size;
};

template <Coordinate T, class Polygon>
OriginLocation locate_origin(
    const Polygon& polygon,
    long double eps
) {
    const int size = polygon.size();
    assert(size >= 3);
    const Point<T> origin;
    const Point<T> base = polygon[0];
    int first = 1;
    if (
        size >= 4 &&
        orientation(base, polygon[1], polygon[2], eps) == 0 &&
        dot(polygon[1] - base, polygon[2] - polygon[1]) > 0
    ) {
        first = 2;
    }
    const int last = size - 1;

    const int first_side = orientation(base, polygon[first], origin, eps);
    const int last_side = orientation(base, polygon[last], origin, eps);
    if (first_side < 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            0,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (last_side > 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            last,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (first_side == 0) {
        if (on_segment(Segment<T>{base, polygon[first]}, origin, eps)) {
            return OriginLocation{
                PointInPolygon::Boundary,
                -1,
                std::array<int, 3>{0, first, 0},
                2,
            };
        }
        return OriginLocation{
            PointInPolygon::Outside,
            first,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (last_side == 0) {
        if (on_segment(Segment<T>{base, polygon[last]}, origin, eps)) {
            return OriginLocation{
                PointInPolygon::Boundary,
                -1,
                std::array<int, 3>{0, last, 0},
                2,
            };
        }
        return OriginLocation{
            PointInPolygon::Outside,
            last - 1,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }

    int left = first;
    int right = last;
    while (right - left >= 2) {
        const int middle = (left + right) / 2;
        if (orientation(base, polygon[middle], origin, eps) >= 0) {
            left = middle;
        } else {
            right = middle;
        }
    }
    const int side = orientation(polygon[left], polygon[right], origin, eps);
    if (side < 0) {
        return OriginLocation{
            PointInPolygon::Outside,
            left,
            std::array<int, 3>{0, 0, 0},
            0,
        };
    }
    if (side == 0) {
        const bool boundary = on_segment(
            Segment<T>{polygon[left], polygon[right]},
            origin,
            eps
        );
        return OriginLocation{
            boundary ? PointInPolygon::Boundary : PointInPolygon::Outside,
            boundary ? -1 : left,
            std::array<int, 3>{left, right, 0},
            boundary ? 2 : 0,
        };
    }
    return OriginLocation{
        PointInPolygon::Inside,
        -1,
        std::array<int, 3>{0, left, right},
        3,
    };
}

template <class Compare>
int periodic_best(int size, Compare better) {
    int left = 0;
    int middle = size;
    int right = 2 * size;
    while (right - left > 2) {
        const int left_middle = (left + middle) / 2;
        const int right_middle = (middle + right + 1) / 2;
        if (better(left_middle % size, middle % size)) {
            right = middle;
            middle = left_middle;
        } else if (better(right_middle % size, middle % size)) {
            left = middle;
            middle = right_middle;
        } else {
            left = left_middle;
            right = right_middle;
        }
    }
    return middle % size;
}

template <Coordinate T, class Polygon>
std::pair<int, int> tangent_vertices_from_origin(
    const Polygon& polygon,
    long double eps
) {
    const int size = polygon.size();
    const Point<T> origin;
    int first = periodic_best(size, [&](int left, int right) {
        return orientation(origin, polygon[left], polygon[right], eps) < 0;
    });
    int second = periodic_best(size, [&](int left, int right) {
        return orientation(origin, polygon[left], polygon[right], eps) > 0;
    });
    const int previous = first == 0 ? size - 1 : first - 1;
    if (orientation(origin, polygon[first], polygon[previous], eps) == 0) {
        first = previous;
    }
    const int next = second + 1 == size ? 0 : second + 1;
    if (orientation(origin, polygon[second], polygon[next], eps) == 0) {
        second = next;
    }
    return std::pair<int, int>(first, second);
}

struct ClosestBoundaryFeature {
    int first;
    int second;
    long double ratio;
    long double distance;
};

template <Coordinate T, class Polygon>
ClosestBoundaryFeature closest_boundary_feature(
    const Polygon& polygon,
    const OriginLocation& location,
    long double eps
) {
    const int size = polygon.size();
    assert(size >= 3);
    assert(location.location == PointInPolygon::Outside);
    const Point<T> origin;

    const auto tangents = tangent_vertices_from_origin<T>(polygon, eps);
    auto visible = [&](int index) {
        return orientation(
            polygon[index],
            polygon[(index + 1) % size],
            origin,
            eps
        ) < 0;
    };
    auto forward_edges = [&](int start, int finish) {
        return finish >= start ? finish - start : finish + size - start;
    };

    int witness = location.outside_edge;
    if (!visible(witness)) {
        const int previous = witness == 0 ? size - 1 : witness - 1;
        const int next = witness + 1 == size ? 0 : witness + 1;
        if (visible(previous)) {
            witness = previous;
        } else if (visible(next)) {
            witness = next;
        }
    }

    int start = tangents.first;
    int finish = tangents.second;
    if (forward_edges(start, witness) >= forward_edges(start, finish)) {
        std::swap(start, finish);
    }
    int edge_count = forward_edges(start, finish);
    if (edge_count == 0) {
        start = location.outside_edge;
        finish = (start + 1) % size;
        edge_count = 1;
    }

    auto vertex = [&](int offset) {
        return polygon[(start + offset) % size];
    };
    int left = 0;
    int right = edge_count;
    while (left < right) {
        const int middle = (left + right) / 2;
        if (norm2(vertex(middle)) <= norm2(vertex(middle + 1))) {
            right = middle;
        } else {
            left = middle + 1;
        }
    }

    ClosestBoundaryFeature result{
        (start + left) % size,
        (start + left) % size,
        0,
        norm(vertex(left)),
    };
    auto consider_edge = [&](int first_offset, int second_offset) {
        const Point<long double> first_point(vertex(first_offset));
        const Point<long double> second_point(vertex(second_offset));
        const Point<long double> direction = second_point - first_point;
        long double ratio =
            -dot(first_point, direction) / dot(direction, direction);
        ratio = std::clamp(ratio, 0.0L, 1.0L);
        const long double candidate_distance =
            norm(first_point + direction * ratio);
        if (candidate_distance < result.distance) {
            result = ClosestBoundaryFeature{
                (start + first_offset) % size,
                (start + second_offset) % size,
                ratio,
                candidate_distance,
            };
        }
    };
    if (left > 0) {
        consider_edge(left - 1, left);
    }
    if (left < edge_count) {
        consider_edge(left, left + 1);
    }
    return result;
}

template <Coordinate T, class Polygon>
long double distance_from_origin(
    const Polygon& polygon,
    long double eps
) {
    const OriginLocation location = locate_origin<T>(polygon, eps);
    if (location.location != PointInPolygon::Outside) return 0;
    return closest_boundary_feature<T>(polygon, location, eps).distance;
}

inline Point<long double> interpolate(
    const Point<long double>& first,
    const Point<long double>& second,
    long double ratio
) {
    return first + (second - first) * ratio;
}

template <Coordinate T>
ClosestPoints
closest_points_from_difference(
    const MinkowskiDifferenceView<T>& difference,
    long double eps
) {
    const OriginLocation location = locate_origin<T>(difference, eps);
    if (location.location == PointInPolygon::Outside) {
        const ClosestBoundaryFeature feature =
            closest_boundary_feature<T>(difference, location, eps);
        const auto first_components = difference.components(feature.first);
        const auto second_components = difference.components(feature.second);
        return ClosestPoints{
            interpolate(
                Point<long double>(first_components.first),
                Point<long double>(second_components.first),
                feature.ratio
            ),
            interpolate(
                Point<long double>(first_components.second),
                Point<long double>(second_components.second),
                feature.ratio
            )
        };
    }

    assert(location.simplex_size == 2 || location.simplex_size == 3);
    std::array<long double, 3> weight{0, 0, 0};
    if (location.simplex_size == 2) {
        const Point<long double> first(difference[location.simplex[0]]);
        const Point<long double> second(difference[location.simplex[1]]);
        const Point<long double> direction = second - first;
        weight[1] = -dot(first, direction) / dot(direction, direction);
        weight[1] = std::clamp(weight[1], 0.0L, 1.0L);
        weight[0] = 1 - weight[1];
    } else {
        const Point<long double> first(difference[location.simplex[0]]);
        const Point<long double> second(difference[location.simplex[1]]);
        const Point<long double> third(difference[location.simplex[2]]);
        const long double denominator = cross(
            second - first,
            third - first
        );
        weight[0] = cross(second, third) / denominator;
        weight[1] = cross(third, first) / denominator;
        weight[2] = cross(first, second) / denominator;
    }

    Point<long double> first_result;
    Point<long double> second_result;
    for (int index = 0; index < location.simplex_size; ++index) {
        const auto components = difference.components(
            location.simplex[index]
        );
        first_result += Point<long double>(components.first) * weight[index];
        second_result +=
            Point<long double>(components.second) * weight[index];
    }
    return ClosestPoints{
        first_result,
        second_result
    };
}

}  // namespace convex_polygon_detail

template <Coordinate T>
std::vector<std::array<Point<T>, 3>> triangulate_convex_polygon(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    polygon = normalize_convex_polygon(std::move(polygon), eps);
    if (polygon.size() < 3) return {};

    std::vector<std::array<Point<T>, 3>> result;
    result.reserve(polygon.size() - 2);
    for (std::size_t index = 1; index + 1 < polygon.size(); ++index) {
        std::array<Point<T>, 3> triangle;
        triangle[0] = polygon[0];
        triangle[1] = polygon[index];
        triangle[2] = polygon[index + 1];
        result.push_back(std::move(triangle));
    }
    return result;
}

template <Coordinate T>
wide_type<T> convex_diameter2(
    std::vector<Point<T>> polygon,
    long double eps = 1e-12L
) {
    polygon = normalize_convex_polygon(std::move(polygon), eps);
    const std::size_t size = polygon.size();
    if (size <= 1) return 0;
    if (size == 2) return distance2(polygon[1], polygon[0]);

    wide_type<T> result = 0;
    std::size_t opposite = 1;
    for (std::size_t index = 0; index < size; ++index) {
        const std::size_t next = (index + 1) % size;
        while (true) {
            const std::size_t candidate = (opposite + 1) % size;
            const auto current_area =
                cross(polygon[index], polygon[next], polygon[opposite]);
            const auto candidate_area =
                cross(polygon[index], polygon[next], polygon[candidate]);
            if (candidate_area <= current_area) break;
            opposite = candidate;
        }
        result = std::max(
            result,
            distance2(polygon[index], polygon[opposite])
        );
        result = std::max(
            result,
            distance2(polygon[next], polygon[opposite])
        );
    }
    return result;
}

template <Coordinate T>
std::vector<Point<long double>> convex_cut(
    const std::vector<Point<T>>& polygon,
    const Line<T>& boundary,
    long double eps = 1e-12L
) {
    assert(boundary.a != boundary.b);
    std::vector<Point<long double>> input;
    input.reserve(polygon.size());
    for (const Point<T>& point : polygon) input.emplace_back(point);
    if (input.empty()) return input;

    const Point<long double> line_start(boundary.a);
    const Point<long double> line_end(boundary.b);
    const Line<long double> line{line_start, line_end};
    std::vector<Point<long double>> result;
    Point<long double> previous = input.back();
    int previous_side = orientation(line_start, line_end, previous, eps);
    for (const Point<long double>& current : input) {
        const int current_side =
            orientation(line_start, line_end, current, eps);
        const bool previous_inside = previous_side >= 0;
        const bool current_inside = current_side >= 0;
        if (previous_inside != current_inside) {
            const Line<long double> crossing{previous, current};
            const LinearIntersection intersection =
                linear_intersection(line, crossing, eps);
            if (intersection.kind == LinearIntersectionKind::Point) {
                result.push_back(intersection.first);
            }
        }
        if (current_inside) result.push_back(current);
        previous = current;
        previous_side = current_side;
    }
    return convex_polygon_detail::clean_polygon(std::move(result), eps);
}

template <Coordinate T>
bool convex_polygons_intersect(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        if (first.size() == 1 && second.size() == 1) {
            return distance(first[0], second[0]) <= eps;
        }
        if (first.size() == 1) {
            return on_segment(
                Segment<T>{second[0], second[1]},
                first[0],
                eps
            );
        }
        if (second.size() == 1) {
            return on_segment(
                Segment<T>{first[0], first[1]},
                second[0],
                eps
            );
        }
        return intersects(
            Segment<T>{first[0], first[1]},
            Segment<T>{second[0], second[1]},
            eps
        );
    }

    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return
        convex_polygon_detail::locate_origin<T>(difference, eps).location !=
        PointInPolygon::Outside;
}

template <Coordinate T>
bool convex_polygons_intersect(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    std::vector<Point<T>> negated;
    negated.reserve(second.size());
    for (const Point<T>& point : second) negated.push_back(-point);
    const std::vector<Point<T>> difference =
        minkowski_sum(first, std::move(negated), eps);
    return
        point_in_convex_polygon(difference, Point<T>(), eps) !=
        PointInPolygon::Outside;
}

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        return closest_points(
            Segment<T>{first[0], first[first.size() - 1]},
            Segment<T>{second[0], second[second.size() - 1]},
            eps
        );
    }
    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return convex_polygon_detail::closest_points_from_difference(
        difference,
        eps
    );
}

template <Coordinate T>
ClosestPoints
convex_polygons_closest_points(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    const ConvexPolygon<T> first_query(first, eps);
    const ConvexPolygon<T> second_query(second, eps);
    return convex_polygons_closest_points(first_query, second_query, eps);
}

template <Coordinate T>
long double convex_polygons_distance(
    const ConvexPolygon<T>& first,
    const ConvexPolygon<T>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    if (first.size() <= 2 && second.size() <= 2) {
        if (convex_polygons_intersect(first, second, eps)) return 0;
        if (first.size() == 1 && second.size() == 1) {
            return distance(first[0], second[0]);
        }
        if (first.size() == 1) {
            return distance(
                Segment<T>{second[0], second[1]},
                first[0]
            );
        }
        if (second.size() == 1) {
            return distance(
                Segment<T>{first[0], first[1]},
                second[0]
            );
        }
        return distance(
            Segment<T>{first[0], first[1]},
            Segment<T>{second[0], second[1]}
        );
    }

    const convex_polygon_detail::MinkowskiDifferenceView<T> difference(
        first,
        second
    );
    return convex_polygon_detail::distance_from_origin<T>(difference, eps);
}

template <Coordinate T>
std::vector<Point<long double>> convex_polygon_intersection(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    using HalfPlane = half_plane_intersection_detail::HalfPlane;
    namespace detail = half_plane_intersection_detail;

    const std::vector<Point<T>> normalized_first =
        normalize_convex_polygon(first, eps);
    const std::vector<Point<T>> normalized_second =
        normalize_convex_polygon(second, eps);
    assert(normalized_first.size() >= 3);
    assert(normalized_second.size() >= 3);
    assert(is_convex_polygon(normalized_first, true, eps));
    assert(is_convex_polygon(normalized_second, true, eps));
    if (!convex_polygons_intersect(
            normalized_first,
            normalized_second,
            eps
        )) {
        return {};
    }

    auto boundaries = [](const std::vector<Point<T>>& polygon) {
        std::vector<HalfPlane> result;
        result.reserve(polygon.size());
        for (std::size_t index = 0; index < polygon.size(); ++index) {
            const Point<long double> point(polygon[index]);
            Point<long double> direction =
                Point<long double>(polygon[(index + 1) % polygon.size()]) -
                point;
            direction = direction / norm(direction);
            result.push_back(HalfPlane{point, direction});
        }
        return result;
    };
    const std::vector<HalfPlane> first_boundaries =
        boundaries(normalized_first);
    const std::vector<HalfPlane> second_boundaries =
        boundaries(normalized_second);

    std::vector<HalfPlane> merged;
    merged.reserve(first_boundaries.size() + second_boundaries.size());
    std::size_t first_index = 0;
    std::size_t second_index = 0;
    while (
        first_index < first_boundaries.size() ||
        second_index < second_boundaries.size()
    ) {
        const bool take_first =
            second_index == second_boundaries.size() ||
            (
                first_index < first_boundaries.size() &&
                detail::direction_less(
                    first_boundaries[first_index],
                    second_boundaries[second_index]
                )
            );
        if (take_first) {
            detail::merge_same_direction(
                merged,
                first_boundaries[first_index++],
                eps
            );
        } else {
            detail::merge_same_direction(
                merged,
                second_boundaries[second_index++],
                eps
            );
        }
    }
    detail::merge_cyclic_ends(merged, eps);

    std::deque<HalfPlane> active;
    for (const HalfPlane& half_plane : merged) {
        while (active.size() >= 2) {
            const std::optional<Point<long double>> point =
                detail::intersection(
                    active[active.size() - 2],
                    active.back(),
                    eps
                );
            if (
                !point.has_value() ||
                !detail::outside(half_plane, *point, eps)
            ) {
                break;
            }
            active.pop_back();
        }
        while (active.size() >= 2) {
            const std::optional<Point<long double>> point =
                detail::intersection(active[0], active[1], eps);
            if (
                !point.has_value() ||
                !detail::outside(half_plane, *point, eps)
            ) {
                break;
            }
            active.pop_front();
        }
        active.push_back(half_plane);
    }
    while (active.size() >= 3) {
        const std::optional<Point<long double>> point =
            detail::intersection(
                active[active.size() - 2],
                active.back(),
                eps
            );
        if (
            !point.has_value() ||
            !detail::outside(active.front(), *point, eps)
        ) {
            break;
        }
        active.pop_back();
    }
    while (active.size() >= 3) {
        const std::optional<Point<long double>> point =
            detail::intersection(active[0], active[1], eps);
        if (
            !point.has_value() ||
            !detail::outside(active.back(), *point, eps)
        ) {
            break;
        }
        active.pop_front();
    }

    std::vector<Point<long double>> result;
    result.reserve(active.size());
    for (std::size_t index = 0; index < active.size(); ++index) {
        const std::optional<Point<long double>> point =
            detail::intersection(
                active[index],
                active[(index + 1) % active.size()],
                eps
            );
        if (point.has_value()) result.push_back(*point);
    }
    return convex_polygon_detail::clean_polygon(std::move(result), eps);
}

template <Coordinate T>
long double convex_polygons_distance(
    const std::vector<Point<T>>& first,
    const std::vector<Point<T>>& second,
    long double eps = 1e-12L
) {
    assert(!first.empty());
    assert(!second.empty());
    std::vector<Point<T>> negated;
    negated.reserve(second.size());
    for (const Point<T>& point : second) negated.push_back(-point);
    const std::vector<Point<T>> difference =
        minkowski_sum(first, std::move(negated), eps);
    const Point<T> origin;
    if (
        point_in_convex_polygon(difference, origin, eps) !=
        PointInPolygon::Outside
    ) {
        return 0;
    }
    if (difference.size() == 1) return distance(difference[0], origin);

    long double result = std::numeric_limits<long double>::infinity();
    for (std::size_t index = 0; index < difference.size(); ++index) {
        if (difference.size() == 2 && index == 1) break;
        result = std::min(
            result,
            distance(
                Segment<T>{
                    difference[index],
                    difference[(index + 1) % difference.size()]
                },
                origin
            )
        );
    }
    return result;
}

}  // namespace geometry
}  // namespace m1une

m1une's library

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

Convex Polygons
(geometry/convex_polygon.hpp)

Overview

Basic Functions

Query Object

Construction and Combination

Vector overloads

Preprocessed pair-query overloads

Example

Depends on

Required by

Verified with

Code

Convex Polygons (geometry/convex_polygon.hpp)

Overview

Basic Functions

Query Object

Construction and Combination

Vector overloads

Preprocessed pair-query overloads

Example

Depends on

Required by

Verified with

Code

Convex Polygons
(geometry/convex_polygon.hpp)