library

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

View the Project on GitHub gyouzasushi/library

:heavy_check_mark: test/unit/geometry.test.cpp

Depends on

Code

#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <random>
#include <set>
#include <vector>

#include "geometry/circle.hpp"
#include "geometry/line.hpp"
#include "geometry/point.hpp"
#include "geometry/polygon.hpp"
using P = geo::point<long long>;
using D = geo::point<double>;
using i128 = __int128;
// independent segment intersection via exact parametric solve
bool brute_intersect(P a, P b, P c, P d) {
    auto cr = [](P p, P q) { return (i128)p.x * q.y - (i128)p.y * q.x; };
    auto on = [&](P p, P q, P x) {
        return cr(q - p, x - p) == 0 && std::min(p.x, q.x) <= x.x && x.x <= std::max(p.x, q.x) &&
               std::min(p.y, q.y) <= x.y && x.y <= std::max(p.y, q.y);
    };
    if (a == b) return on(c, d, a);
    if (c == d) return on(a, b, c);
    P r = b - a, s = d - c;
    i128 den = cr(r, s);
    if (den != 0) {
        i128 tn = cr(c - a, s), un = cr(c - a, r);
        if (den < 0) den = -den, tn = -tn, un = -un;
        return 0 <= tn && tn <= den && 0 <= un && un <= den;
    }
    if (cr(c - a, r) != 0) return false;  // parallel, not collinear
    return on(a, b, c) || on(a, b, d) || on(c, d, a) || on(c, d, b);
}
bool in_triangle_closed(P a, P b, P c, P p) {
    int s1 = geo::ccw(a, b, p), s2 = geo::ccw(b, c, p), s3 = geo::ccw(c, a, p);
    if (geo::ccw(a, b, c) == 0) return false;
    return (s1 >= 0 && s2 >= 0 && s3 >= 0) || (s1 <= 0 && s2 <= 0 && s3 <= 0);
}
int main() {
    std::mt19937 rng(2024);
    auto rp = [&](int w) { return P(int(rng() % (2 * w + 1)) - w, int(rng() % (2 * w + 1)) - w); };
    for (int it = 0; it < 3000; it++) {
        // ccw sign == sign of cross, and matches double version
        P a = rp(5), b = rp(5), c = rp(5);
        D da(a.x, a.y), db(b.x, b.y), dc(c.x, c.y);
        assert(geo::ccw(a, b, c) == geo::sgn(geo::cross(b - a, c - a)));
        assert(geo::ccw(a, b, c) == geo::ccw(da, db, dc));
        // segment intersection vs independent implementation
        P d = rp(5);
        geo::segment<long long> s(a, b), t(c, d);
        assert(geo::intersect(s, t) == brute_intersect(a, b, c, d));
        assert(geo::intersect(s, t) == geo::intersect(geo::segment<double>(da, db), geo::segment<double>(dc, D(d.x, d.y))));
        // on_segment
        assert(geo::on_segment(s, c) == brute_intersect(a, b, c, c));
    }
    for (int it = 0; it < 2000; it++) {
        int n = 1 + rng() % 9;
        std::vector<P> ps(n);
        for (auto& p : ps) p = rp(4);
        // strict hull: vertices are exactly the points not inside/on any triangle or segment of others
        geo::polygon<long long> ch = geo::convex_hull(ps);
        std::set<std::pair<long long, long long>> hull_set;
        for (auto& p : ch) hull_set.insert({p.x, p.y});
        std::set<std::pair<long long, long long>> uniq;
        for (auto& p : ps) uniq.insert({p.x, p.y});
        std::vector<P> us;
        for (auto [x, y] : uniq) us.emplace_back(x, y);
        int m = int(us.size());
        for (int i = 0; i < m; i++) {
            bool vertex = true;
            for (int j = 0; j < m && vertex; j++) {
                if (j == i) continue;
                for (int k = j + 1; k < m && vertex; k++) {
                    if (k == i) continue;
                    if (geo::on_segment(geo::segment<long long>(us[j], us[k]), us[i])) vertex = false;
                    for (int l = k + 1; l < m && vertex; l++) {
                        if (l == i) continue;
                        if (in_triangle_closed(us[j], us[k], us[l], us[i])) vertex = false;
                    }
                }
            }
            assert(vertex == hull_set.count({us[i].x, us[i].y}));
        }
        // ccw order and convexity
        int h = int(ch.size());
        if (h >= 3) {
            assert(geo::is_convex(ch));
            for (int i = 0; i < h; i++) assert(geo::ccw(ch[i], ch[(i + 1) % h], ch[(i + 2) % h]) > 0);
            assert(geo::area2(ch) > 0);
            // contains on convex polygon vs brute
            for (int q = 0; q < 20; q++) {
                P p = rp(5);
                int expect = 2;
                for (int i = 0; i < h; i++) {
                    int s = geo::ccw(ch[i], ch[(i + 1) % h], p);
                    if (s < 0) expect = 0;
                    if (s == 0 && expect != 0 && geo::on_segment(geo::segment<long long>(ch[i], ch[(i + 1) % h]), p)) expect = 1;
                }
                assert(geo::contains(ch, p) == expect);
            }
            // area2 vs double shoelace
            geo::polygon<double> dch;
            for (auto& p : ch) dch.emplace_back(p.x, p.y);
            assert(std::abs(geo::area(dch) * 2 - (double)geo::area2(ch)) < 1e-9);
        }
        // non-strict hull: superset of strict, all on boundary, contains everything
        geo::polygon<long long> ch2 = geo::convex_hull(ps, false);
        assert(ch2.size() >= ch.size());
        if (h >= 3) {
            for (auto& p : ch2) assert(geo::contains(ch, p) == 1);
            for (auto& p : ps) assert(geo::contains(ch2, p) >= 1);
            int h2 = int(ch2.size());
            for (int i = 0; i < h2; i++) assert(geo::ccw(ch2[i], ch2[(i + 1) % h2], ch2[(i + 2) % h2]) >= 0);
        }
        // closest / farthest pair vs brute
        if (n >= 2) {
            auto [i, j] = geo::closest_pair(ps);
            auto [k, l] = geo::farthest_pair(ps);
            assert(i != j && k != l);
            long long best = geo::norm(ps[i] - ps[j]), worst = geo::norm(ps[k] - ps[l]);
            for (int x = 0; x < n; x++) {
                for (int y = x + 1; y < n; y++) {
                    assert(geo::norm(ps[x] - ps[y]) >= best);
                    assert(geo::norm(ps[x] - ps[y]) <= worst);
                }
            }
        }
        // arg sort vs atan2 (-pi, pi], origin treated as 0
        std::vector<P> qs = ps;
        std::sort(qs.begin(), qs.end(), geo::arg_less{});
        auto ang = [](P p) { return p.x == 0 && p.y == 0 ? 0.0L : atan2l((long double)p.y, (long double)p.x); };
        for (int i = 0; i + 1 < n; i++) assert(ang(qs[i]) <= ang(qs[i + 1]) + 1e-15L);
    }
    // non-convex polygon containment (handmade)
    geo::polygon<long long> u = {{0, 0}, {4, 0}, {4, 4}, {3, 4}, {3, 1}, {1, 1}, {1, 4}, {0, 4}};
    assert(geo::contains(u, {2, 2}) == 0);
    assert(geo::contains(u, {2, 0}) == 1);
    assert(geo::contains(u, {3, 2}) == 1);
    assert(geo::contains(u, {2, 1}) == 1);
    assert(geo::contains(u, {3, 1}) == 1);
    assert(geo::contains(u, {1, 2}) == 1);
    assert(geo::contains(u, {0, 2}) == 1);
    assert(geo::contains(u, {2, 3}) == 0);
    assert(geo::contains(u, {3, 5}) == 0);
    for (long long x = -1; x <= 5; x++) {
        for (long long y = -1; y <= 5; y++) {
            P p(x, y);
            bool inside = (0 < x && x < 4 && 0 < y && y < 4) && !(1 <= x && x <= 3 && 1 <= y);
            bool boundary = !inside && ((0 <= x && x <= 4 && (y == 0 || y == 4) && !(y == 4 && 1 < x && x < 3)) ||
                                        ((x == 0 || x == 4) && 0 <= y && y <= 4) ||
                                        ((x == 1 || x == 3) && 1 <= y && y <= 4) || (y == 1 && 1 <= x && x <= 3));
            assert(geo::contains(u, p) == (inside ? 2 : boundary ? 1 : 0));
        }
    }
    // double-only: circle sanity
    geo::circle<double> o(D(0, 0), 1);
    assert(geo::contains(o, D(0.5, 0)) == 2 && geo::contains(o, D(1, 0)) == 1 && geo::contains(o, D(2, 0)) == 0);
    auto cp = geo::cross_points(o, geo::line<double>(D(-2, 0), D(2, 0)));
    assert(cp.size() == 2 && cp[0] == D(-1, 0) && cp[1] == D(1, 0));
    assert(std::abs(geo::area_of_intersection(o, geo::circle<double>(D(0, 0), 1)) - M_PI) < 1e-9);
    assert(std::abs(geo::area_of_intersection(o, geo::circle<double>(D(2, 0), 1))) < 1e-9);
    std::cout << "Hello World" << std::endl;
}
#line 1 "test/unit/geometry.test.cpp"
#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iostream>
#include <random>
#include <set>
#include <vector>

#line 5 "geometry/circle.hpp"
#include <concepts>
#include <numbers>
#line 9 "geometry/circle.hpp"

#line 7 "geometry/line.hpp"

#line 5 "geometry/point.hpp"
#include <type_traits>
namespace geo {
inline long double eps = 1e-9;
template <class T> int sgn(T a) {
    if constexpr (std::is_integral_v<T>) {
        return (a > 0) - (a < 0);
    } else {
        return a < -eps ? -1 : a > eps ? 1 : 0;
    }
}
template <class T> struct point {
    T x, y;
    point() : x(0), y(0) {
    }
    point(T x, T y) : 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*=(T k) {
        x *= k, y *= k;
        return *this;
    }
    point& operator/=(T k) requires std::floating_point<T> {
        x /= k, y /= k;
        return *this;
    }
    point operator-() const {
        return point(-x, -y);
    }
    friend point operator+(point a, const point& b) {
        return a += b;
    }
    friend point operator-(point a, const point& b) {
        return a -= b;
    }
    friend point operator*(point a, T k) {
        return a *= k;
    }
    friend point operator*(T k, point a) {
        return a *= k;
    }
    friend point operator/(point a, T k) requires std::floating_point<T> {
        return a /= k;
    }
    friend bool operator==(const point& a, const point& b) {
        return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
    }
    friend bool operator!=(const point& a, const point& b) {
        return !(a == b);
    }
    friend bool operator<(const point& a, const point& b) {
        if (sgn(a.x - b.x) != 0) return a.x < b.x;
        return sgn(a.y - b.y) < 0;
    }
    friend bool operator>(const point& a, const point& b) {
        return b < a;
    }
    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 <class T> T dot(const point<T>& a, const point<T>& b) {
    return a.x * b.x + a.y * b.y;
}
template <class T> T cross(const point<T>& a, const point<T>& b) {
    return a.x * b.y - a.y * b.x;
}
template <class T> T norm(const point<T>& p) {
    return p.x * p.x + p.y * p.y;
}
template <std::floating_point T> T abs(const point<T>& p) {
    return std::sqrt(norm(p));
}
template <std::floating_point T> T arg(const point<T>& p) {
    return std::atan2(p.y, p.x);
}
template <class T> point<T> normal(const point<T>& p) {
    return point<T>(-p.y, p.x);
}
template <std::floating_point T> point<T> unit(const point<T>& p) {
    return p / abs(p);
}
template <std::floating_point T> point<T> rotate(const point<T>& p, T theta) {
    T c = std::cos(theta), s = std::sin(theta);
    return point<T>(p.x * c - p.y * s, p.x * s + p.y * c);
}
template <class T> int ccw(const point<T>& a, const point<T>& b, const point<T>& c) {
    return sgn(cross(b - a, c - a));
}
struct arg_less {
    template <class T> static int group(const point<T>& p) {
        if (sgn(p.y) < 0) return 0;
        if (sgn(p.y) == 0 && sgn(p.x) >= 0) return 1;
        return 2;
    }
    template <class T> bool operator()(const point<T>& a, const point<T>& b) const {
        if (group(a) != group(b)) return group(a) < group(b);
        return sgn(cross(a, b)) > 0;
    }
};
}  // namespace geo
#line 9 "geometry/line.hpp"
namespace geo {
template <class T> struct line {
    point<T> a, b;
    line() {
    }
    line(const point<T>& a, const point<T>& b) : a(a), b(b) {
    }
    friend std::istream& operator>>(std::istream& is, line& l) {
        return is >> l.a >> l.b;
    }
};
template <class T> struct segment {
    point<T> a, b;
    segment() {
    }
    segment(const point<T>& a, const point<T>& b) : a(a), b(b) {
    }
    friend std::istream& operator>>(std::istream& is, segment& s) {
        return is >> s.a >> s.b;
    }
};
template <class T> line<T> to_line(const segment<T>& s) {
    return line<T>(s.a, s.b);
}
template <class T> bool parallel(const line<T>& l, const line<T>& m) {
    return sgn(cross(l.b - l.a, m.b - m.a)) == 0;
}
template <class T> bool orthogonal(const line<T>& l, const line<T>& m) {
    return sgn(dot(l.b - l.a, m.b - m.a)) == 0;
}
template <class T> bool on_line(const line<T>& l, const point<T>& p) {
    return ccw(l.a, l.b, p) == 0;
}
template <class T> bool on_segment(const segment<T>& s, const point<T>& p) {
    return ccw(s.a, s.b, p) == 0 && sgn(dot(s.a - p, s.b - p)) <= 0;
}
template <class T> bool intersect(const line<T>& l, const line<T>& m) {
    return !parallel(l, m) || on_line(l, m.a);
}
template <class T> bool intersect(const line<T>& l, const segment<T>& s) {
    return ccw(l.a, l.b, s.a) * ccw(l.a, l.b, s.b) <= 0;
}
template <class T> bool intersect(const segment<T>& s, const line<T>& l) {
    return intersect(l, s);
}
template <class T> bool intersect(const segment<T>& s, const segment<T>& t) {
    if (ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) > 0) return false;
    if (ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) > 0) return false;
    if (ccw(s.a, s.b, t.a) != 0 || ccw(s.a, s.b, t.b) != 0) return true;
    return on_segment(s, t.a) || on_segment(s, t.b) || on_segment(t, s.a) || on_segment(t, s.b);
}
template <std::floating_point T> point<T> projection(const line<T>& l, const point<T>& p) {
    point<T> d = l.b - l.a;
    return l.a + d * (dot(d, p - l.a) / norm(d));
}
template <std::floating_point T> point<T> reflection(const line<T>& l, const point<T>& p) {
    return projection(l, p) * T(2) - p;
}
template <std::floating_point T> line<T> perpendicular_bisector(const point<T>& a, const point<T>& b) {
    point<T> m = (a + b) / T(2);
    return line<T>(m, m + normal(b - a));
}
template <std::floating_point T> T dist(const point<T>& a, const point<T>& b) {
    return abs(a - b);
}
template <std::floating_point T> T dist(const line<T>& l, const point<T>& p) {
    return std::abs(cross(l.b - l.a, p - l.a)) / abs(l.b - l.a);
}
template <std::floating_point T> T dist(const point<T>& p, const line<T>& l) {
    return dist(l, p);
}
template <std::floating_point T> T dist(const segment<T>& s, const point<T>& p) {
    if (sgn(dot(s.b - s.a, p - s.a)) < 0) return abs(p - s.a);
    if (sgn(dot(s.a - s.b, p - s.b)) < 0) return abs(p - s.b);
    return dist(to_line(s), p);
}
template <std::floating_point T> T dist(const point<T>& p, const segment<T>& s) {
    return dist(s, p);
}
template <std::floating_point T> T dist(const line<T>& l, const line<T>& m) {
    if (intersect(l, m)) return 0;
    return dist(l, m.a);
}
template <std::floating_point T> T dist(const line<T>& l, const segment<T>& s) {
    if (intersect(l, s)) return 0;
    return std::min(dist(l, s.a), dist(l, s.b));
}
template <std::floating_point T> T dist(const segment<T>& s, const line<T>& l) {
    return dist(l, s);
}
template <std::floating_point T> T dist(const segment<T>& s, const segment<T>& t) {
    if (intersect(s, t)) return 0;
    return std::min({dist(s, t.a), dist(s, t.b), dist(t, s.a), dist(t, s.b)});
}
template <std::floating_point T> point<T> cross_point(const line<T>& l, const line<T>& m) {
    assert(!parallel(l, m));
    point<T> d = l.b - l.a;
    return l.a + d * (cross(m.b - m.a, m.a - l.a) / cross(m.b - m.a, d));
}
template <std::floating_point T> point<T> cross_point(const line<T>& l, const segment<T>& s) {
    assert(intersect(l, s));
    return cross_point(l, to_line(s));
}
template <std::floating_point T> point<T> cross_point(const segment<T>& s, const line<T>& l) {
    return cross_point(l, s);
}
template <std::floating_point T> point<T> cross_point(const segment<T>& s, const segment<T>& t) {
    assert(intersect(s, t));
    return cross_point(to_line(s), to_line(t));
}
}  // namespace geo
#line 6 "geometry/polygon.hpp"
#include <utility>
#line 8 "geometry/polygon.hpp"

#line 11 "geometry/polygon.hpp"
namespace geo {
template <class T> using polygon = std::vector<point<T>>;
template <class T> T area2(const polygon<T>& poly) {
    int n = int(poly.size());
    if constexpr (std::is_integral_v<T>) {
        __int128 ret = 0;
        for (int i = 0; i < n; i++) ret += (__int128)cross(poly[i], poly[(i + 1) % n]);
        return T(ret);
    } else {
        T ret = 0;
        for (int i = 0; i < n; i++) ret += cross(poly[i], poly[(i + 1) % n]);
        return ret;
    }
}
template <std::floating_point T> T area(const polygon<T>& poly) {
    return area2(poly) / T(2);
}
template <class T> bool is_convex(const polygon<T>& poly) {
    int n = int(poly.size());
    for (int i = 0; i < n; i++) {
        if (ccw(poly[i], poly[(i + 1) % n], poly[(i + 2) % n]) < 0) return false;
    }
    return true;
}
template <class T> int contains(const polygon<T>& poly, const point<T>& p) {
    int n = int(poly.size());
    bool in = false;
    for (int i = 0; i < n; i++) {
        point<T> a = poly[i] - p, b = poly[(i + 1) % n] - p;
        if (ccw(a, b, point<T>()) == 0 && sgn(dot(a, b)) <= 0) return 1;
        if (sgn(a.y) > sgn(b.y)) std::swap(a, b);
        if (sgn(a.y) <= 0 && 0 < sgn(b.y) && sgn(cross(a, b)) > 0) in = !in;
    }
    return in ? 2 : 0;
}
template <class T> polygon<T> convex_hull(std::vector<point<T>> ps, bool strict = true) {
    int n = int(ps.size());
    std::sort(ps.begin(), ps.end());
    ps.erase(std::unique(ps.begin(), ps.end()), ps.end());
    n = int(ps.size());
    if (n <= 2) return ps;
    polygon<T> ret(2 * n);
    int k = 0;
    for (int i = 0; i < n; ret[k++] = ps[i++]) {
        while (k >= 2 && ccw(ret[k - 2], ret[k - 1], ps[i]) < (strict ? 1 : 0)) k--;
    }
    for (int i = n - 2, t = k + 1; i >= 0; ret[k++] = ps[i--]) {
        while (k >= t && ccw(ret[k - 2], ret[k - 1], ps[i]) < (strict ? 1 : 0)) k--;
    }
    ret.resize(k - 1);
    return ret;
}
template <class T> std::pair<int, int> farthest_pair(const std::vector<point<T>>& ps) {
    int n = int(ps.size());
    assert(n >= 2);
    polygon<T> ch = convex_hull(ps);
    int m = int(ch.size());
    auto index_of = [&](const point<T>& p) {
        return int(std::find(ps.begin(), ps.end(), p) - ps.begin());
    };
    if (m == 1) return {0, 1};
    if (m == 2) return {index_of(ch[0]), index_of(ch[1])};
    std::pair<int, int> ret = {0, 1};
    T best = norm(ch[0] - ch[1]);
    for (int i = 0, j = 1; i < m; i++) {
        while (norm(ch[i] - ch[(j + 1) % m]) > norm(ch[i] - ch[j])) j = (j + 1) % m;
        if (norm(ch[i] - ch[j]) > best) {
            best = norm(ch[i] - ch[j]);
            ret = {i, j};
        }
    }
    return {index_of(ch[ret.first]), index_of(ch[ret.second])};
}
template <class T> std::pair<int, int> closest_pair(const std::vector<point<T>>& ps) {
    int n = int(ps.size());
    assert(n >= 2);
    std::vector<int> id(n);
    for (int i = 0; i < n; i++) id[i] = i;
    std::sort(id.begin(), id.end(), [&](int i, int j) { return ps[i] < ps[j]; });
    std::pair<int, int> ret = {id[0], id[1]};
    T best = norm(ps[id[0]] - ps[id[1]]);
    std::vector<int> buf(n);
    auto rec = [&](auto rec, int l, int r) -> void {
        if (r - l <= 1) return;
        int m = (l + r) / 2;
        T x = ps[id[m]].x;
        rec(rec, l, m);
        rec(rec, m, r);
        std::inplace_merge(id.begin() + l, id.begin() + m, id.begin() + r,
                           [&](int i, int j) { return ps[i].y < ps[j].y; });
        int k = 0;
        for (int t = l; t < r; t++) {
            int i = id[t];
            T dx = ps[i].x - x;
            if (dx * dx >= best) continue;
            for (int u = k - 1; u >= 0; u--) {
                int j = buf[u];
                T dy = ps[i].y - ps[j].y;
                if (dy * dy >= best) break;
                T d = norm(ps[i] - ps[j]);
                if (d < best) {
                    best = d;
                    ret = {j, i};
                }
            }
            buf[k++] = i;
        }
    };
    rec(rec, 0, n);
    return ret;
}
template <std::floating_point T> polygon<T> convex_cut(const polygon<T>& poly, const line<T>& l) {
    int n = int(poly.size());
    polygon<T> ret;
    for (int i = 0; i < n; i++) {
        const point<T>& p = poly[i];
        const point<T>& q = poly[(i + 1) % n];
        int cp = ccw(l.a, l.b, p), cq = ccw(l.a, l.b, q);
        if (cp >= 0) ret.push_back(p);
        if (cp * cq < 0) ret.push_back(cross_point(l, line<T>(p, q)));
    }
    return ret;
}
}  // namespace geo
#line 13 "geometry/circle.hpp"
namespace geo {
template <std::floating_point T> struct circle {
    point<T> c;
    T r;
    circle() : r(0) {
    }
    circle(const point<T>& c, T r) : c(c), r(r) {
    }
    friend std::istream& operator>>(std::istream& is, circle& o) {
        return is >> o.c >> o.r;
    }
};
template <std::floating_point T> T area(const circle<T>& o) {
    return std::numbers::pi_v<T> * o.r * o.r;
}
template <std::floating_point T> int contains(const circle<T>& o, const point<T>& p) {
    int s = sgn(abs(p - o.c) - o.r);
    return s > 0 ? 0 : s == 0 ? 1 : 2;
}
template <std::floating_point T> bool intersect(const line<T>& l, const circle<T>& o) {
    return sgn(dist(l, o.c) - o.r) <= 0;
}
template <std::floating_point T> bool intersect(const circle<T>& o, const line<T>& l) {
    return intersect(l, o);
}
template <std::floating_point T> bool intersect(const segment<T>& s, const circle<T>& o) {
    if (sgn(dist(s, o.c) - o.r) > 0) return false;
    return sgn(abs(s.a - o.c) - o.r) >= 0 || sgn(abs(s.b - o.c) - o.r) >= 0;
}
template <std::floating_point T> bool intersect(const circle<T>& o, const segment<T>& s) {
    return intersect(s, o);
}
template <std::floating_point T> int count_common_tangents(const circle<T>& a, const circle<T>& b) {
    T d = abs(a.c - b.c);
    if (sgn(d - (a.r + b.r)) > 0) return 4;
    if (sgn(d - (a.r + b.r)) == 0) return 3;
    if (sgn(d - std::abs(a.r - b.r)) > 0) return 2;
    if (sgn(d - std::abs(a.r - b.r)) == 0) return 1;
    return 0;
}
template <std::floating_point T> bool intersect(const circle<T>& a, const circle<T>& b) {
    int k = count_common_tangents(a, b);
    return 1 <= k && k <= 3;
}
template <std::floating_point T> std::vector<point<T>> cross_points(const line<T>& l, const circle<T>& o) {
    if (!intersect(l, o)) return {};
    point<T> p = projection(l, o.c);
    T h = std::sqrt(std::max<T>(0, o.r * o.r - norm(p - o.c)));
    point<T> v = unit(l.b - l.a) * h;
    std::vector<point<T>> ret = {p - v, p + v};
    if (ret[1] < ret[0]) std::swap(ret[0], ret[1]);
    return ret;
}
template <std::floating_point T> std::vector<point<T>> cross_points(const circle<T>& o, const line<T>& l) {
    return cross_points(l, o);
}
template <std::floating_point T> std::vector<point<T>> cross_points(const segment<T>& s, const circle<T>& o) {
    std::vector<point<T>> ret = cross_points(to_line(s), o);
    ret.erase(std::remove_if(ret.begin(), ret.end(), [&](const point<T>& p) { return !on_segment(s, p); }),
              ret.end());
    return ret;
}
template <std::floating_point T> std::vector<point<T>> cross_points(const circle<T>& o, const segment<T>& s) {
    return cross_points(s, o);
}
template <std::floating_point T> std::vector<point<T>> cross_points(const circle<T>& a, const circle<T>& b) {
    if (!intersect(a, b)) return {};
    T d = abs(a.c - b.c);
    T x = (d * d + a.r * a.r - b.r * b.r) / (2 * d);
    T h = std::sqrt(std::max<T>(0, a.r * a.r - x * x));
    point<T> e = unit(b.c - a.c);
    point<T> m = a.c + e * x, v = normal(e) * h;
    std::vector<point<T>> ret = {m - v, m + v};
    if (ret[1] < ret[0]) std::swap(ret[0], ret[1]);
    return ret;
}
template <std::floating_point T> std::vector<point<T>> tangent_points(const circle<T>& o, const point<T>& p) {
    assert(sgn(abs(p - o.c) - o.r) >= 0);
    T r = std::sqrt(std::max<T>(0, norm(p - o.c) - o.r * o.r));
    return cross_points(o, circle<T>(p, r));
}
template <std::floating_point T> std::vector<point<T>> common_tangents(const circle<T>& a, const circle<T>& b) {
    std::vector<point<T>> ret;
    T d = abs(a.c - b.c);
    if (sgn(d - std::abs(a.r - b.r)) >= 0) {
        T r = std::sqrt(std::max<T>(0, d * d - (a.r - b.r) * (a.r - b.r) + b.r * b.r));
        std::vector<point<T>> ps = cross_points(a, circle<T>(b.c, r));
        ret.insert(ret.end(), ps.begin(), ps.end());
    }
    if (sgn(d - (a.r + b.r)) >= 0) {
        point<T> p = a.c + (b.c - a.c) * (a.r / (a.r + b.r));
        std::vector<point<T>> ps = tangent_points(a, p);
        ret.insert(ret.end(), ps.begin(), ps.end());
    }
    std::sort(ret.begin(), ret.end());
    ret.erase(std::unique(ret.begin(), ret.end()), ret.end());
    return ret;
}
template <std::floating_point T> circle<T> incircle(const point<T>& a, const point<T>& b, const point<T>& c) {
    T x = abs(b - c), y = abs(c - a), z = abs(a - b);
    point<T> p = (a * x + b * y + c * z) / (x + y + z);
    return circle<T>(p, dist(line<T>(a, b), p));
}
template <std::floating_point T> circle<T> circumcircle(const point<T>& a, const point<T>& b, const point<T>& c) {
    point<T> p = cross_point(perpendicular_bisector(a, b), perpendicular_bisector(a, c));
    return circle<T>(p, abs(a - p));
}
template <std::floating_point T> circle<T> smallest_enclosing_circle(std::vector<point<T>> ps) {
    int n = int(ps.size());
    assert(n >= 1);
    std::shuffle(ps.begin(), ps.end(), std::mt19937(123456789));
    circle<T> ret(ps[0], 0);
    for (int i = 1; i < n; i++) {
        if (contains(ret, ps[i])) continue;
        ret = circle<T>(ps[i], 0);
        for (int j = 0; j < i; j++) {
            if (contains(ret, ps[j])) continue;
            ret = circle<T>((ps[i] + ps[j]) / T(2), abs(ps[i] - ps[j]) / T(2));
            for (int k = 0; k < j; k++) {
                if (contains(ret, ps[k])) continue;
                ret = circumcircle(ps[i], ps[j], ps[k]);
            }
        }
    }
    return ret;
}
template <std::floating_point T> T area_of_intersection(const circle<T>& o, const polygon<T>& poly) {
    int n = int(poly.size());
    auto tri = [](const point<T>& a, const point<T>& b) { return cross(a, b); };
    auto sector = [&](const point<T>& a, const point<T>& b) { return o.r * o.r * arg(rotate(b, -arg(a))); };
    auto inside = [&](const point<T>& a) { return sgn(abs(a) - o.r) < 0; };
    T ret = 0;
    for (int i = 0; i < n; i++) {
        point<T> a = poly[i] - o.c, b = poly[(i + 1) % n] - o.c;
        segment<T> s(a, b);
        circle<T> c(point<T>(), o.r);
        if (!intersect(s, c)) {
            ret += inside(a) ? tri(a, b) : sector(a, b);
            continue;
        }
        std::vector<point<T>> ps = cross_points(s, c);
        point<T> p = ps[0], q = ps[ps.size() - 1];
        if (sgn(dot(b - a, q - p)) < 0) std::swap(p, q);
        ret += inside(a) ? tri(a, p) : sector(a, p);
        ret += tri(p, q);
        ret += inside(b) ? tri(q, b) : sector(q, b);
    }
    return ret / 2;
}
template <std::floating_point T> T area_of_intersection(const polygon<T>& poly, const circle<T>& o) {
    return area_of_intersection(o, poly);
}
template <std::floating_point T> T area_of_intersection(const circle<T>& a, const circle<T>& b) {
    T d = abs(a.c - b.c);
    if (sgn(a.r + b.r - d) <= 0) return 0;
    if (sgn(std::abs(a.r - b.r) - d) >= 0) return std::min(area(a), area(b));
    auto part = [&](const circle<T>& x, const circle<T>& y) {
        T cs = (d * d + x.r * x.r - y.r * y.r) / (2 * d * x.r);
        T theta = std::acos(std::clamp<T>(cs, -1, 1));
        return x.r * x.r * (theta - std::sin(2 * theta) / 2);
    };
    return part(a, b) + part(b, a);
}
template <std::floating_point T> circle<T> circle_of_apollonius(const point<T>& a, const point<T>& b, T n, T m) {
    assert(sgn(n - m) != 0);
    point<T> c = (a * (-m * m) + b * (n * n)) / (n * n - m * m);
    return circle<T>(c, std::sqrt(abs(a - c) * abs(b - c)));
}
}  // namespace geo
#line 14 "test/unit/geometry.test.cpp"
using P = geo::point<long long>;
using D = geo::point<double>;
using i128 = __int128;
// independent segment intersection via exact parametric solve
bool brute_intersect(P a, P b, P c, P d) {
    auto cr = [](P p, P q) { return (i128)p.x * q.y - (i128)p.y * q.x; };
    auto on = [&](P p, P q, P x) {
        return cr(q - p, x - p) == 0 && std::min(p.x, q.x) <= x.x && x.x <= std::max(p.x, q.x) &&
               std::min(p.y, q.y) <= x.y && x.y <= std::max(p.y, q.y);
    };
    if (a == b) return on(c, d, a);
    if (c == d) return on(a, b, c);
    P r = b - a, s = d - c;
    i128 den = cr(r, s);
    if (den != 0) {
        i128 tn = cr(c - a, s), un = cr(c - a, r);
        if (den < 0) den = -den, tn = -tn, un = -un;
        return 0 <= tn && tn <= den && 0 <= un && un <= den;
    }
    if (cr(c - a, r) != 0) return false;  // parallel, not collinear
    return on(a, b, c) || on(a, b, d) || on(c, d, a) || on(c, d, b);
}
bool in_triangle_closed(P a, P b, P c, P p) {
    int s1 = geo::ccw(a, b, p), s2 = geo::ccw(b, c, p), s3 = geo::ccw(c, a, p);
    if (geo::ccw(a, b, c) == 0) return false;
    return (s1 >= 0 && s2 >= 0 && s3 >= 0) || (s1 <= 0 && s2 <= 0 && s3 <= 0);
}
int main() {
    std::mt19937 rng(2024);
    auto rp = [&](int w) { return P(int(rng() % (2 * w + 1)) - w, int(rng() % (2 * w + 1)) - w); };
    for (int it = 0; it < 3000; it++) {
        // ccw sign == sign of cross, and matches double version
        P a = rp(5), b = rp(5), c = rp(5);
        D da(a.x, a.y), db(b.x, b.y), dc(c.x, c.y);
        assert(geo::ccw(a, b, c) == geo::sgn(geo::cross(b - a, c - a)));
        assert(geo::ccw(a, b, c) == geo::ccw(da, db, dc));
        // segment intersection vs independent implementation
        P d = rp(5);
        geo::segment<long long> s(a, b), t(c, d);
        assert(geo::intersect(s, t) == brute_intersect(a, b, c, d));
        assert(geo::intersect(s, t) == geo::intersect(geo::segment<double>(da, db), geo::segment<double>(dc, D(d.x, d.y))));
        // on_segment
        assert(geo::on_segment(s, c) == brute_intersect(a, b, c, c));
    }
    for (int it = 0; it < 2000; it++) {
        int n = 1 + rng() % 9;
        std::vector<P> ps(n);
        for (auto& p : ps) p = rp(4);
        // strict hull: vertices are exactly the points not inside/on any triangle or segment of others
        geo::polygon<long long> ch = geo::convex_hull(ps);
        std::set<std::pair<long long, long long>> hull_set;
        for (auto& p : ch) hull_set.insert({p.x, p.y});
        std::set<std::pair<long long, long long>> uniq;
        for (auto& p : ps) uniq.insert({p.x, p.y});
        std::vector<P> us;
        for (auto [x, y] : uniq) us.emplace_back(x, y);
        int m = int(us.size());
        for (int i = 0; i < m; i++) {
            bool vertex = true;
            for (int j = 0; j < m && vertex; j++) {
                if (j == i) continue;
                for (int k = j + 1; k < m && vertex; k++) {
                    if (k == i) continue;
                    if (geo::on_segment(geo::segment<long long>(us[j], us[k]), us[i])) vertex = false;
                    for (int l = k + 1; l < m && vertex; l++) {
                        if (l == i) continue;
                        if (in_triangle_closed(us[j], us[k], us[l], us[i])) vertex = false;
                    }
                }
            }
            assert(vertex == hull_set.count({us[i].x, us[i].y}));
        }
        // ccw order and convexity
        int h = int(ch.size());
        if (h >= 3) {
            assert(geo::is_convex(ch));
            for (int i = 0; i < h; i++) assert(geo::ccw(ch[i], ch[(i + 1) % h], ch[(i + 2) % h]) > 0);
            assert(geo::area2(ch) > 0);
            // contains on convex polygon vs brute
            for (int q = 0; q < 20; q++) {
                P p = rp(5);
                int expect = 2;
                for (int i = 0; i < h; i++) {
                    int s = geo::ccw(ch[i], ch[(i + 1) % h], p);
                    if (s < 0) expect = 0;
                    if (s == 0 && expect != 0 && geo::on_segment(geo::segment<long long>(ch[i], ch[(i + 1) % h]), p)) expect = 1;
                }
                assert(geo::contains(ch, p) == expect);
            }
            // area2 vs double shoelace
            geo::polygon<double> dch;
            for (auto& p : ch) dch.emplace_back(p.x, p.y);
            assert(std::abs(geo::area(dch) * 2 - (double)geo::area2(ch)) < 1e-9);
        }
        // non-strict hull: superset of strict, all on boundary, contains everything
        geo::polygon<long long> ch2 = geo::convex_hull(ps, false);
        assert(ch2.size() >= ch.size());
        if (h >= 3) {
            for (auto& p : ch2) assert(geo::contains(ch, p) == 1);
            for (auto& p : ps) assert(geo::contains(ch2, p) >= 1);
            int h2 = int(ch2.size());
            for (int i = 0; i < h2; i++) assert(geo::ccw(ch2[i], ch2[(i + 1) % h2], ch2[(i + 2) % h2]) >= 0);
        }
        // closest / farthest pair vs brute
        if (n >= 2) {
            auto [i, j] = geo::closest_pair(ps);
            auto [k, l] = geo::farthest_pair(ps);
            assert(i != j && k != l);
            long long best = geo::norm(ps[i] - ps[j]), worst = geo::norm(ps[k] - ps[l]);
            for (int x = 0; x < n; x++) {
                for (int y = x + 1; y < n; y++) {
                    assert(geo::norm(ps[x] - ps[y]) >= best);
                    assert(geo::norm(ps[x] - ps[y]) <= worst);
                }
            }
        }
        // arg sort vs atan2 (-pi, pi], origin treated as 0
        std::vector<P> qs = ps;
        std::sort(qs.begin(), qs.end(), geo::arg_less{});
        auto ang = [](P p) { return p.x == 0 && p.y == 0 ? 0.0L : atan2l((long double)p.y, (long double)p.x); };
        for (int i = 0; i + 1 < n; i++) assert(ang(qs[i]) <= ang(qs[i + 1]) + 1e-15L);
    }
    // non-convex polygon containment (handmade)
    geo::polygon<long long> u = {{0, 0}, {4, 0}, {4, 4}, {3, 4}, {3, 1}, {1, 1}, {1, 4}, {0, 4}};
    assert(geo::contains(u, {2, 2}) == 0);
    assert(geo::contains(u, {2, 0}) == 1);
    assert(geo::contains(u, {3, 2}) == 1);
    assert(geo::contains(u, {2, 1}) == 1);
    assert(geo::contains(u, {3, 1}) == 1);
    assert(geo::contains(u, {1, 2}) == 1);
    assert(geo::contains(u, {0, 2}) == 1);
    assert(geo::contains(u, {2, 3}) == 0);
    assert(geo::contains(u, {3, 5}) == 0);
    for (long long x = -1; x <= 5; x++) {
        for (long long y = -1; y <= 5; y++) {
            P p(x, y);
            bool inside = (0 < x && x < 4 && 0 < y && y < 4) && !(1 <= x && x <= 3 && 1 <= y);
            bool boundary = !inside && ((0 <= x && x <= 4 && (y == 0 || y == 4) && !(y == 4 && 1 < x && x < 3)) ||
                                        ((x == 0 || x == 4) && 0 <= y && y <= 4) ||
                                        ((x == 1 || x == 3) && 1 <= y && y <= 4) || (y == 1 && 1 <= x && x <= 3));
            assert(geo::contains(u, p) == (inside ? 2 : boundary ? 1 : 0));
        }
    }
    // double-only: circle sanity
    geo::circle<double> o(D(0, 0), 1);
    assert(geo::contains(o, D(0.5, 0)) == 2 && geo::contains(o, D(1, 0)) == 1 && geo::contains(o, D(2, 0)) == 0);
    auto cp = geo::cross_points(o, geo::line<double>(D(-2, 0), D(2, 0)));
    assert(cp.size() == 2 && cp[0] == D(-1, 0) && cp[1] == D(1, 0));
    assert(std::abs(geo::area_of_intersection(o, geo::circle<double>(D(0, 0), 1)) - M_PI) < 1e-9);
    assert(std::abs(geo::area_of_intersection(o, geo::circle<double>(D(2, 0), 1))) < 1e-9);
    std::cout << "Hello World" << std::endl;
}
Back to top page