geometry/all.hpp

View this file on GitHub
Last update: 2024-03-24 14:28:32+09:00
Include: #include "geometry/all.hpp"

Depends on

Code

#pragma once

#include "geometry/point.hpp"
#include "geometry/line.hpp"
#include "geometry/segment.hpp"
#include "geometry/circle.hpp"
#include "geometry/polygon.hpp"

#include "geometry/projection.hpp"
#include "geometry/reflection.hpp"
#include "geometry/ccw.hpp"
#include "geometry/is_orthogonal.hpp"
#include "geometry/is_parallel.hpp"

#include "geometry/is_intersect_ll.hpp"
#include "geometry/is_intersect_lp.hpp"
#include "geometry/is_intersect_ss.hpp"
#include "geometry/is_intersect_sp.hpp"
#include "geometry/tangent_number_cc.hpp"
#include "geometry/is_intersect_cc.hpp"
#include "geometry/is_intersect_cp.hpp"
#include "geometry/is_intersect_cl.hpp"

#include "geometry/cross_point_ll.hpp"
#include "geometry/cross_point_ss.hpp"
#include "geometry/cross_point_cl.hpp"
#include "geometry/cross_point_cc.hpp"

#include "geometry/distance_lp.hpp"
#include "geometry/distance_sp.hpp"
#include "geometry/distance_ss.hpp"

#include "geometry/tangent_point_cp.hpp"
#include "geometry/incircle.hpp"
#include "geometry/circumscribed_circle.hpp"

#include "geometry/polygon_area.hpp"
#include "geometry/polygon_is_convex.hpp"
#include "geometry/polygon_contain.hpp"
#include "geometry/monotone_chain.hpp"
#include "geometry/convex_polygon_diameter.hpp"
#include "geometry/convex_polygon_cut.hpp"

#include "geometry/closest_pair.hpp"
#include "geometry/farthest_pair.hpp"

#line 2 "geometry/all.hpp"

#line 2 "geometry/point.hpp"

// point
template <typename T> struct Point {
    static T EPS;
    static constexpr T PI = 3.1415926535'8979323846'2643383279L;
    static void set_eps(const T &e) { EPS = e; }
    T x, y;
    Point(const T x = T(0), const T y = T(0)) : x(x), y(y) {}
    Point &operator+=(const Point &p) {
        x += p.x;
        y += p.y;
        return *this;
    }
    Point &operator-=(const Point &p) {
        x -= p.x;
        y -= p.y;
        return *this;
    }
    Point &operator*=(const Point &p) { return *this = Point(x * p.x - y * p.y, x * p.y + y * p.x); }
    Point &operator*=(const T &k) {
        x *= k;
        y *= k;
        return *this;
    }
    Point &operator/=(const Point &p) { return *this = Point(x * p.x + y * p.y, -x * p.y + y * p.x) / (p.x * p.x + p.y * p.y); }
    Point &operator/=(const T &k) {
        x /= k;
        y /= k;
        return *this;
    }

    Point operator+() const { return *this; }
    Point operator-() const { return Point(-x, -y); }

    friend Point operator+(const Point &a, const Point &b) { return Point(a) += b; }
    friend Point operator-(const Point &a, const Point &b) { return Point(a) -= b; }
    friend Point operator*(const Point &a, const Point &b) { return Point(a) *= b; }
    friend Point operator*(const Point &p, const T &k) { return Point(p) *= k; }
    friend Point operator/(const Point &a, const Point &b) { return Point(a) /= b; }
    friend Point operator/(const Point &p, const T &k) { return Point(p) /= k; }
    // for std::set, std::map, compare_arg, ...
    friend bool operator<(const Point &a, const Point &b) { return a.x == b.x ? a.y < b.y : a.x < b.x; }
    friend bool operator>(const Point &a, const Point &b) { return a.x == b.x ? a.y > b.y : a.x > b.x; }
    // I/O
    friend std::istream &operator>>(std::istream &is, Point &p) { return is >> p.x >> p.y; }
    friend std::ostream &operator<<(std::ostream &os, const Point &p) { return os << '(' << p.x << ' ' << p.y << ')'; }
};

// template
template <typename T> inline int sign(const T &x) { return x < -Point<T>::EPS ? -1 : (x > Point<T>::EPS ? 1 : 0); }
template <typename T> inline bool equal(const T &a, const T &b) { return sign(a - b) == 0; }
template <typename T> inline T radian_to_degree(const T &r) { return r * 180.0 / Point<T>::PI; }
template <typename T> inline T degree_to_radian(const T &d) { return d * Point<T>::PI / 180.0; }

// contain enum
constexpr int IN = 2;
constexpr int ON = 1;
constexpr int OUT = 0;

// equal (point and point)
template <typename T> inline bool equal(const Point<T> &a, const Point<T> &b) { return equal(a.x, b.x) and equal(a.y, b.y); }
// inner product
template <typename T> inline T dot(const Point<T> &a, const Point<T> &b) { return a.x * b.x + a.y * b.y; }
// outer product
template <typename T> inline T cross(const Point<T> &a, const Point<T> &b) { return a.x * b.y - a.y * b.x; }
// rotate Point p counterclockwise by theta radian
template <typename T> inline Point<T> rotate(const Point<T> &p, const T &theta) { return p * Point<T>(std::cos(theta), std::sin(theta)); }
// compare (x, y)
template <typename T> inline bool compare_x(const Point<T> &a, const Point<T> &b) { return equal(a.x, b.x) ? sign(a.y - b.y) < 0 : sign(a.x - b.x) < 0; }
// compare (y, x)
template <typename T> inline bool compare_y(const Point<T> &a, const Point<T> &b) { return equal(a.y, b.y) ? sign(a.x - b.x) < 0 : sign(a.y - b.y) < 0; }
// compare by (arg(p), norm(p)) [0, 360)
template <typename T> inline bool compare_arg(const Point<T> &a, const Point<T> &b) {
    // https://ngtkana.hatenablog.com/entry/2021/11/13/202103
    assert(!equal(a, Point<T>(0, 0)));
    assert(!equal(b, Point<T>(0, 0)));
    if ((Point<T>(0, 0) < Point<T>(a.y, a.x)) == (Point<T>(0, 0) < Point<T>(b.y, b.x))) {
        return (a.x * b.y == a.y * b.x) ? norm(a) < norm(b) : a.x * b.y > a.y * b.x;
    } else {
        return Point<T>(a.y, a.x) > Point<T>(b.y, b.x);
    }
}
// |p| ^ 2
template <typename T> inline T norm(const Point<T> &p) { return p.x * p.x + p.y * p.y; }
// |p|
template <typename T> inline T abs(const Point<T> &p) { return std::sqrt(norm(p)); }
// arg
template <typename T> inline T arg(const Point<T> &p) { return std::atan2(p.y, p.x); }
// polar
template <typename T> inline Point<T> polar(const T &rho, const T &theta = T(0)) { return rotate(Point<T>(rho, 0), theta); }
// EPS
template <> double Point<double>::EPS = 1e-9;
template <> long double Point<long double>::EPS = 1e-12;
template <> long long Point<long long>::EPS = 0;
template <> __int128_t Point<__int128_t>::EPS = 0;
// change EPS
// using Double = double;
// using Pt = Point<Double>;
// Point<Double>::set_eps(new_eps);
#line 2 "geometry/line.hpp"

#line 4 "geometry/line.hpp"

// line
template <typename T> struct Line {
    Point<T> a, b;

    Line() = default;

    Line(const Point<T> &a, const Point<T> &b) : a(a), b(b) {}

    // Ax + By = C
    Line(const T &A, const T &B, const T &C) {
        assert(!(equal(A, T(0)) and equal(B, T(0))));
        if (equal(A, T(0))) {
            a = Point<T>(T(0), C / B), b = Point<T>(T(1), C / B);
        } else if (equal(B, T(0))) {
            a = Point<T>(C / A, T(0)), b = Point<T>(C / A, T(1));
        } else if (equal(C, T(0))) {
            a = Point<T>(T(0), T(0)), b = Point<T>(T(1), B / A);
        } else {
            a = Point<T>(T(0), C / B), b = Point<T>(C / A, T(0));
        }
    }

    friend std::istream &operator>>(std::istream &is, Line &p) { return is >> p.a >> p.b; }
    friend std::ostream &operator<<(std::ostream &os, const Line &p) { return os << p.a << "->" << p.b; }
};
#line 2 "geometry/segment.hpp"

#line 4 "geometry/segment.hpp"

// segment
template <typename T> struct Segment : Line<T> {
    Segment() = default;

    Segment(const Point<T> &a, const Point<T> &b) : Line<T>(a, b) {}
};
#line 2 "geometry/circle.hpp"

#line 4 "geometry/circle.hpp"

// circle
template <typename T> struct Circle {
    Point<T> o;
    T r;

    Circle() = default;

    Circle(const Point<T> &o, const T &r) : o(o), r(r) {}

    friend std::istream &operator>>(std::istream &is, Circle &c) { return is >> c.o >> c.r; }  // format : x y r
    friend std::ostream &operator<<(std::ostream &os, const Circle &c) { return os << c.o << ' ' << c.r; }
};
#line 2 "geometry/polygon.hpp"

#line 4 "geometry/polygon.hpp"

// polygon
template <typename T> using Polygon = std::vector<Point<T>>;
template <typename T> std::istream &operator>>(std::istream &is, Polygon<T> &p) {
    for (auto &&pi : p) is >> pi;
    return is;
}
template <typename T> std::ostream &operator<<(std::ostream &os, const Polygon<T> &p) {
    for (auto &&pi : p) os << pi << " -> ";
    return os;
}
#line 8 "geometry/all.hpp"

#line 2 "geometry/projection.hpp"

#line 5 "geometry/projection.hpp"

// projection
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_A
template <typename T> Point<T> projection(const Line<T> &l, const Point<T> &p) {
    T t = dot(p - l.a, l.b - l.a) / norm(l.b - l.a);
    return l.a + t * (l.b - l.a);
}
#line 2 "geometry/reflection.hpp"

#line 6 "geometry/reflection.hpp"

// reflection
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_B
template <typename T> Point<T> reflection(const Line<T> &l, const Point<T> &p) { return p + (projection(l, p) - p) * T(2); }
#line 2 "geometry/ccw.hpp"

#line 4 "geometry/ccw.hpp"

// counter clockwise
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C
constexpr int COUNTER_CLOCKWISE = 1;  // a-b-c counter clockwise
constexpr int CLOCKWISE = -1;         // a-b-c clockwise
constexpr int ONLINE_BACK = 2;        // c-a-b line
constexpr int ONLINE_FRONT = -2;      // a-b-c line
constexpr int ON_SEGMENT = 0;         // a-c-b line
template <typename T> int ccw(const Point<T> &a, Point<T> b, Point<T> c) {
    b = b - a, c = c - a;
    if (sign(cross(b, c)) == 1) return COUNTER_CLOCKWISE;
    if (sign(cross(b, c)) == -1) return CLOCKWISE;
    if (sign(dot(b, c)) == -1) return ONLINE_BACK;
    if (norm(b) < norm(c)) return ONLINE_FRONT;
    return ON_SEGMENT;
}
#line 2 "geometry/is_orthogonal.hpp"

#line 4 "geometry/is_orthogonal.hpp"

// orthogonal
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <typename T> inline bool is_orthogonal(const Line<T> &l1, const Line<T> &l2) { return sign(dot(l1.b - l1.a, l2.b - l2.a)) == 0; }
#line 2 "geometry/is_parallel.hpp"

#line 4 "geometry/is_parallel.hpp"

// parallel
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <typename T> inline bool is_parallel(const Line<T> &l1, const Line<T> &l2) { return sign(cross(l1.b - l1.a, l2.b - l2.a)) == 0; }
#line 14 "geometry/all.hpp"

#line 2 "geometry/is_intersect_ll.hpp"

#line 4 "geometry/is_intersect_ll.hpp"

// intersection (line and line)
template <typename T> bool is_intersect_ll(const Line<T> &l1, const Line<T> &l2) {
    Point<T> base = l1.b - l1.a;
    T d12 = cross(base, l2.b - l2.a);
    T d1 = cross(base, l1.b - l2.a);
    if (sign(d12) == 0) {
        // parallel
        if (sign(d1) == 0) {
            // cross
            return true;
        } else {
            // not cross
            return false;
        }
    }
    return true;
}
#line 2 "geometry/is_intersect_lp.hpp"

#line 5 "geometry/is_intersect_lp.hpp"

// intersection (line and point)
// ccw(a, b, c) == ON_SEGMENT or ONLINE_BACK or ONLINE_FRONT
template <typename T> inline bool is_intersect_lp(const Line<T> &l, const Point<T> &p) {
    int res = ccw(l.a, l.b, p);
    return (res == ONLINE_BACK or res == ONLINE_FRONT or res == ON_SEGMENT);
}
#line 2 "geometry/is_intersect_ss.hpp"

#line 5 "geometry/is_intersect_ss.hpp"

// intersection (segment and segment)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_B
template <typename T> inline bool is_intersect_ss(const Segment<T> &s1, const Segment<T> &s2) { return (ccw(s1.a, s1.b, s2.a) * ccw(s1.a, s1.b, s2.b) <= 0 and ccw(s2.a, s2.b, s1.a) * ccw(s2.a, s2.b, s1.b) <= 0); }
#line 2 "geometry/is_intersect_sp.hpp"

#line 5 "geometry/is_intersect_sp.hpp"

// intersection (segment and point)
// ccw(a, b, c) == ON_SEGMENT -> a - c - b
template <typename T> inline bool is_intersect_sp(const Segment<T> &s, const Point<T> &p) { return ccw(s.a, s.b, p) == ON_SEGMENT or sign(norm(s.a - p)) == 0 or sign(norm(s.b - p)) == 0; }
#line 3 "geometry/tangent_number_cc.hpp"

// return the number of tangent
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_A
template <typename T> int tangent_number_cc(Circle<T> c1, Circle<T> c2) {
    if (c1.r < c2.r) std::swap(c1, c2);
    const T d2 = norm(c1.o - c2.o);
    if (sign(d2 - (c1.r + c2.r) * (c1.r + c2.r)) == 1) return 4;  // d > c1.r + c2.r and c1.r + c2.r >= 0 <=> d * d > (c1.r + c2.r) * (c1.r + c2.r)
    if (sign(d2 - (c1.r + c2.r) * (c1.r + c2.r)) == 0) return 3;  // d = c1.r + c2.r and c1.r + c2.r >= 0 <=> d * d = (c1.r + c2.r) * (c1.r + c2.r)
    if (sign(d2 - (c1.r - c2.r) * (c1.r - c2.r)) == 1) return 2;  // d > c1.r - c2.r and c1.r - c2.r >= 0 <=> d * d > (c1.r - c2.r) * (c1.r - c2.r)
    if (sign(d2 - (c1.r - c2.r) * (c1.r - c2.r)) == 0) return 1;  // d = c1.r - c2.r and c1.r - c2.r >= 0 <=> d * d = (c1.r - c2.r) * (c1.r - c2.r)
    return 0;
}
#line 2 "geometry/is_intersect_cc.hpp"

#line 5 "geometry/is_intersect_cc.hpp"

// intersection (circle and circle)
// intersect = number of tangent is 1, 2, 3
template <typename T> inline bool is_intersect_cc(const Circle<T> &c1, const Circle<T> &c2) {
    int num = tangent_number_cc(c1, c2);
    return 1 <= num and num <= 3;
}
#line 2 "geometry/is_intersect_cp.hpp"

#line 5 "geometry/is_intersect_cp.hpp"
// intersection (circle and point)
template <typename T> inline bool is_intersect_cp(const Circle<T> &c, const Point<T> &p) { return equal(norm(p - c.o), c.r * c.r); }
#line 2 "geometry/is_intersect_cl.hpp"

#line 2 "geometry/distance_lp.hpp"

#line 6 "geometry/distance_lp.hpp"

// distance (line and point) (Double = double or long)
template <typename T> T distance_lp(const Line<T> &l, const Point<T> &p) { return abs(p - projection(l, p)); }
template <typename T> T distance2_lp(const Line<T> &l, const Point<T> &p) { return norm(p - projection(l, p)); }
#line 5 "geometry/is_intersect_cl.hpp"

// intersection (circle and line)
template <typename T> inline bool is_intersect_cl(const Circle<T> &c, const Line<T> &l) { return sign(c.r * c.r - distance2_lp(l, c.o)) >= 0; }
#line 23 "geometry/all.hpp"

#line 2 "geometry/cross_point_ll.hpp"

#line 4 "geometry/cross_point_ll.hpp"

// cross point (line and line)
template <typename T> Point<T> cross_point_ll(const Line<T> &l1, const Line<T> &l2) {
    Point<T> base = l1.b - l1.a;
    T d12 = cross(base, l2.b - l2.a);
    T d1 = cross(base, l1.b - l2.a);
    if (sign(d12) == 0) {
        // parallel
        if (sign(d1) == 0) {
            // cross
            return l2.a;
        } else {
            // not cross
            assert(false);
        }
    }
    return l2.a + (l2.b - l2.a) * (d1 / d12);
}
#line 2 "geometry/cross_point_ss.hpp"

#line 6 "geometry/cross_point_ss.hpp"

// cross point (segment and segment)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_C
template <typename T> Point<T> cross_point_ss(const Segment<T> &s1, const Segment<T> &s2) {
    // check intersection s1 and s2
    assert(is_intersect_ss(s1, s2));
    Point<T> base = s1.b - s1.a;
    T d12 = cross(base, s2.b - s2.a);
    T d1 = cross(base, s1.b - s2.a);
    if (sign(d12) == 0) {
        // parallel
        if (sign(d1) == 0) {
            // equal
            if (is_intersect_sp(s1, s2.a)) return s2.a;
            if (is_intersect_sp(s1, s2.b)) return s2.b;
            if (is_intersect_sp(s2, s1.a)) return s1.a;
            assert(is_intersect_sp(s2, s1.b));
            return s1.b;
        } else {
            // excepted by is_intersect_ss(s1, s2)
            assert(0);
        }
    }
    return s2.a + (s2.b - s2.a) * (d1 / d12);
}
#line 2 "geometry/cross_point_cl.hpp"

#line 5 "geometry/cross_point_cl.hpp"

// cross point (circle and line)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_D
template <typename T> std::vector<Point<T>> cross_point_cl(const Circle<T> &c, const Line<T> &l) {
    assert(is_intersect_cl(c, l));
    auto pr = projection(l, c.o);
    if (equal(norm(pr - c.o), c.r * c.r)) return {pr};
    Point<T> e = (l.b - l.a) * (T(1) / abs(l.b - l.a));
    auto k = sqrt(c.r * c.r - norm(pr - c.o));
    return {pr - e * k, pr + e * k};
}
#line 2 "geometry/cross_point_cc.hpp"

#line 4 "geometry/cross_point_cc.hpp"

// cross point (circle and circle)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_E
template <typename T> std::vector<Point<T>> cross_point_cc(const Circle<T> &c1, const Circle<T> &c2) {
    if (!is_intersect_cc(c1, c2)) return {};
    T d = abs(c1.o - c2.o);
    T a = std::acos((c1.r * c1.r - c2.r * c2.r + d * d) / (T(2) * c1.r * d));
    T t = arg(c2.o - c1.o);
    Point<T> p = c1.o + polar(c1.r, t + a);
    Point<T> q = c1.o + polar(c1.r, t - a);
    if (equal(p.x, q.x) and equal(p.y, q.y)) return {p};
    return {p, q};
}
#line 28 "geometry/all.hpp"

#line 2 "geometry/distance_sp.hpp"

#line 7 "geometry/distance_sp.hpp"

// distance (segment and point)
template <typename T> T distance_sp(const Segment<T> &s, const Point<T> &p) {
    Point<T> r = projection(s, p);
    if (is_intersect_sp(s, r)) {
        return abs(r - p);
    }
    return std::min(abs(s.a - p), abs(s.b - p));
}
#line 2 "geometry/distance_ss.hpp"

#line 6 "geometry/distance_ss.hpp"

// distance (segment and segment)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_D
template <typename T> T distance_ss(const Segment<T> &s1, const Segment<T> &s2) {
    if (is_intersect_ss(s1, s2)) return T(0);
    return std::min({distance_sp(s1, s2.a), distance_sp(s1, s2.b), distance_sp(s2, s1.a), distance_sp(s2, s1.b)});
}
#line 32 "geometry/all.hpp"

#line 2 "geometry/tangent_point_cp.hpp"

#line 4 "geometry/tangent_point_cp.hpp"

// tangent point (circle and point)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_F
template <typename T> std::pair<Point<T>, Point<T>> tangent_point_cp(const Circle<T> &c, const Point<T> &p) {
    assert(sign(abs(c.o - p) - c.r) == 1);
    auto res = cross_point_cc(c, Circle(p, sqrt(norm(c.o - p) - c.r * c.r)));
    return {res[0], res[1]};
}
#line 2 "geometry/incircle.hpp"

#line 6 "geometry/incircle.hpp"

// incircle of a triangle
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_B
// https://drken1215.hatenablog.com/entry/2020/10/16/073700
template <typename T> Circle<T> incircle(const Point<T> &a, const Point<T> &b, const Point<T> &c) {
    T A = arg((c - a) / (b - a)), B = arg((a - b) / (c - b));
    Line l1(a, a + rotate(b - a, A / 2)), l2(b, b + rotate(c - b, B / 2));
    auto o = cross_point_ll(l1, l2);
    auto r = distance_lp(Line(a, b), o);
    return Circle(o, r);
}
#line 2 "geometry/circumscribed_circle.hpp"

#line 5 "geometry/circumscribed_circle.hpp"

// circumscribed circle of a triangle
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_C
// https://drken1215.hatenablog.com/entry/2020/10/16/074400
template <typename T> Circle<T> circumscribed_circle(const Point<T> &a, const Point<T> &b, const Point<T> &c) {
    Line l1((a + b) / 2, (a + b) / 2 + rotate(b - a, Point<T>::PI / 2)), l2((b + c) / 2, (b + c) / 2 + rotate(c - b, Point<T>::PI / 2));
    auto o = cross_point_ll(l1, l2);
    auto r = abs(o - a);
    return Circle(o, r);
}
#line 36 "geometry/all.hpp"

#line 2 "geometry/polygon_area.hpp"

#line 4 "geometry/polygon_area.hpp"

// area of polygon
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_A
// return area * 2
template <typename T> T polygon_area2(const Polygon<T> &p) {
    int n = (int)p.size();
    assert(n >= 2);
    T ret = T(0);
    for (int i = 0; i < n - 1; i++) {
        ret += cross(p[i], p[i + 1]);
    }
    ret += cross(p[n - 1], p[0]);
    // counter clockwise: ret > 0
    // clockwise: ret < 0
    return ret;
}
template <typename T> T polygon_area(const Polygon<T> &p) { return polygon_area2(p) / T(2); }
#line 2 "geometry/polygon_is_convex.hpp"

#line 5 "geometry/polygon_is_convex.hpp"

// check polygon is convex (not strictly, 0 <= angle <= 180 degrees)
// angle = 180 degrees -> ON_SEGMENT
// angle = 0 degrees -> ONLINE_FRONT or ONLINE_BACK
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_B
template <typename T> bool polygon_is_convex(const Polygon<T> &p) {
    int n = (int)p.size();
    assert(n >= 3);
    bool okccw = true, okcw = true;
    for (int i = 0; i < n - 2; i++) {
        int res = ccw(p[i], p[i + 1], p[i + 2]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    {
        int res = ccw(p[n - 2], p[n - 1], p[0]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    {
        int res = ccw(p[n - 1], p[0], p[1]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    return true;
}
#line 2 "geometry/polygon_contain.hpp"

#line 5 "geometry/polygon_contain.hpp"

// polygon contain point -> 2 (IN)
// polygon cross point -> 1 (ON)
// otherwise -> 0 (OUT)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_C
template <typename T> int polygon_contain(const Polygon<T> &q, const Point<T> &p) {
    bool x = false;
    int n = (int)q.size();
    for (int i = 0; i < n - 1; i++) {
        if (is_intersect_sp(Segment(q[i], q[i + 1]), p)) return ON;
        Point a = q[i] - p, b = q[i + 1] - p;
        if (a.y > b.y) std::swap(a, b);
        // a.y < b.y
        // check each point's y is 0 at most 1 times
        if (sign(a.y) <= 0 and sign(b.y) > 0 and sign(cross(a, b)) > 0) x = !x;
    }
    {
        if (is_intersect_sp(Segment(q[n - 1], q[0]), p)) return ON;
        Point a = q[n - 1] - p, b = q[0] - p;
        if (a.y > b.y) std::swap(a, b);
        if (sign(a.y) <= 0 and sign(b.y) > 0 and sign(cross(a, b)) > 0) x = !x;
    }
    return (x ? IN : OUT);
}
#line 2 "geometry/monotone_chain.hpp"

#line 5 "geometry/monotone_chain.hpp"

// convex hull (Andrew's monotone chain convex hull algorithm)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_A
// sort (x, y) by lexicographical order, use stack, calculate upper convex hull and lower convex hull
// counter clockwise order
// assume the return value of ccw is not ONLINE_BACK or ONLINE_FRONT (lexicographical order)
// strict is true : points on the edges of the convex hull are not included (the number of points is minimized)
// complexity: O(n \log n) (n: the number of points)
template <typename T> Polygon<T> monotone_chain(std::vector<Point<T>> &p, bool strict = true) {
    int n = (int)p.size();
    if (n <= 2) return p;
    std::sort(p.begin(), p.end(), compare_x<T>);
    Polygon<T> r;
    r.reserve(n * 2);
    if (strict) {
        for (int i = 0; i < n; i++) {
            while (r.size() >= 2 and ccw(r[r.size() - 2], r[r.size() - 1], p[i]) != CLOCKWISE) {
                r.pop_back();
            }
            r.push_back(p[i]);
        }
        int t = r.size() + 1;
        for (int i = n - 2; i >= 0; i--) {
            while (r.size() >= t and ccw(r[r.size() - 2], r[r.size() - 1], p[i]) != CLOCKWISE) {
                r.pop_back();
            }
            r.push_back(p[i]);
        }
    } else {
        for (int i = 0; i < n; i++) {
            while (r.size() >= 2 and ccw(r[r.size() - 2], r[r.size() - 1], p[i]) == COUNTER_CLOCKWISE) {
                r.pop_back();
            }
            r.push_back(p[i]);
        }
        int t = r.size() + 1;
        for (int i = n - 2; i >= 0; i--) {
            while (r.size() >= t and ccw(r[r.size() - 2], r[r.size() - 1], p[i]) == COUNTER_CLOCKWISE) {
                r.pop_back();
            }
            r.push_back(p[i]);
        }
    }
    r.pop_back();
    std::reverse(r.begin(), r.end());
    return r;
}
#line 2 "geometry/convex_polygon_diameter.hpp"

#line 5 "geometry/convex_polygon_diameter.hpp"

// convex polygon diameter
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_B
// return {index1, index2, diameter}
// using the method of rotating calipers (https://en.wikipedia.org/wiki/Rotating_calipers)
// complexity: O(n)
template <typename T> std::tuple<int, int, T> convex_polygon_diameter(const Polygon<T> &p) {
    assert(polygon_is_convex(p));
    int n = (int)p.size();
    assert(n >= 2);
    if (n == 2) {
        return {0, 1, abs(p[0] - p[1])};
    }
    auto [it_min, it_max] = std::minmax_element(p.begin(), p.end(), compare_x<T>);
    int idx_min = it_min - p.begin();
    int idx_max = it_max - p.begin();

    T maxdis = norm(p[idx_max] - p[idx_min]);
    int maxi = idx_min, i = idx_min, maxj = idx_max, j = idx_max;
    do {
        int ni = (i + 1 == n ? 0 : i + 1), nj = (j + 1 == n ? 0 : j + 1);
        if (sign(cross(p[ni] - p[i], p[nj] - p[j])) < 0) {
            i = ni;
        } else {
            j = nj;
        }
        if (norm(p[i] - p[j]) > maxdis) {
            maxdis = norm(p[i] - p[j]);
            maxi = i;
            maxj = j;
        }
    } while (i != idx_min || j != idx_max);
    return {maxi, maxj, abs(p[maxi] - p[maxj])};
}
#line 2 "geometry/convex_polygon_cut.hpp"

#line 5 "geometry/convex_polygon_cut.hpp"

// cut convex polygon p by line l
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_C
// return {left polygon, right polygon}
// whether each point is included is determined by the sign of the outer product of the two vectors to the endpoints of the line
template <typename T> std::pair<Polygon<T>, Polygon<T>> convex_polygon_cut(const Polygon<T> &p, const Line<T> &l) {
    int n = (int)p.size();
    assert(n >= 3);
    Polygon<T> pl, pr;
    for (int i = 0; i < n; i++) {
        int s1 = sign(cross(l.a - p[i], l.b - p[i]));
        int s2 = sign(cross(l.a - p[i + 1 == n ? 0 : i + 1], l.b - p[i + 1 == n ? 0 : i + 1]));
        if (s1 >= 0) {
            pl.push_back(p[i]);
        }
        if (s1 <= 0) {
            pr.push_back(p[i]);
        }
        if (s1 * s2 < 0) {
            // don't use "<=", use "<" to exclude endpoints
            auto pc = cross_point_ll(Line(p[i], p[i + 1 == n ? 0 : i + 1]), l);
            pl.push_back(pc);
            pr.push_back(pc);
        }
    }
    return {pl, pr};
}
#line 43 "geometry/all.hpp"

#line 2 "geometry/closest_pair.hpp"

#line 4 "geometry/closest_pair.hpp"

// closest pair
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_5_A
// return {index1, index2, distance}
// using divide-and-conquer algorithm
// complexity: O(n \log n) (n: the number of points)
template <typename T> std::tuple<int, int, T> closest_pair(const std::vector<Point<T>> &p) {
    int n = int(p.size());
    assert(n >= 2);
    if (n == 2) {
        return {0, 1, abs(p[0] - p[1])};
    }
    // may not be efficient due to indirect references ...
    std::vector<int> ind(n);
    std::iota(ind.begin(), ind.end(), 0);
    std::sort(ind.begin(), ind.end(), [&](int i, int j) { return compare_x(p[i], p[j]); });
    auto divide_and_conquer = [&](auto f, int l, int r) -> std::tuple<int, int, T> {
        if (r - l <= 1) return {-1, -1, std::numeric_limits<T>::max()};
        int md = (l + r) / 2;
        T x = p[ind[md]].x;
        // divide and conquer
        auto [i1l, i2l, dl] = f(f, l, md);
        auto [i1r, i2r, dr] = f(f, md, r);
        int i1, i2;
        T d;
        if (dl < dr) {
            d = dl, i1 = i1l, i2 = i2l;
        } else {
            d = dr, i1 = i1r, i2 = i2r;
        }
        std::inplace_merge(ind.begin() + l, ind.begin() + md, ind.begin() + r, [&](int i, int j) { return compare_y(p[i], p[j]); });
        // ind are sorted by y
        std::vector<int> near_x;  // index of vertices whose distance from the line x is less than d
        for (int i = l; i < r; i++) {
            if (sign(std::abs(p[ind[i]].x - x) - d) >= 0) continue;  // std::abs(p[ind[i]].x - x) >= d
            int sz = int(near_x.size());
            // iterate from the end until the distance in y-coordinates is greater than or equal to d
            for (int j = sz - 1; j >= 0; j--) {
                Point cp = p[ind[i]] - p[near_x[j]];
                if (sign(cp.y - d) >= 0) break;  // cp.y >= d
                T cd = abs(cp);
                if (cd < d) {
                    d = cd, i1 = ind[i], i2 = near_x[j];
                }
            }
            near_x.push_back(ind[i]);
        }
        return {i1, i2, d};
    };
    return divide_and_conquer(divide_and_conquer, 0, n);
}
#line 2 "geometry/farthest_pair.hpp"

#line 5 "geometry/farthest_pair.hpp"

// farthest pair
// return {index1, index2, distance}
// using monotone chain (convex hull) and convex polygon diameter
// complexity: O(n \log n) (n: the number of points)
template <typename T> std::tuple<int, int, T> farthest_pair(const std::vector<Point<T>> &p) {
    int n = int(p.size());
    assert(n >= 2);
    if (n == 2) {
        return {0, 1, abs(p[0] - p[1])};
    }
    auto q = p;
    auto ch = monotone_chain(q);                   // O(n \log n)
    auto [i, j, d] = convex_polygon_diameter(ch);  // O(|ch|)
    int resi, resj;
    for (int k = 0; k < n; k++) {
        if (p[k] == ch[i]) {
            resi = k;
        }
        if (p[k] == ch[j]) {
            resj = k;
        }
    }
    return {resi, resj, d};
}
#line 46 "geometry/all.hpp"