cp-library

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

View the Project on GitHub rniya/cp-library

:heavy_check_mark: Contain
(src/geometry/contain.hpp)

Depends on

Verified with

Code

#pragma once
#include "Circle.hpp"
#include "Polygon.hpp"
#include "distance.hpp"

namespace geometry {

enum class Containment { OUT, ON, IN };

template <typename T> Containment contain(const Circle<T>& c, const Point<T>& p) {
    int x = sgn(distance(c.c, p) - c.r);
    return x == -1 ? Containment::IN : x == 1 ? Containment::OUT : Containment::ON;
}

template <typename T> Containment contain(const Polygon<T>& P, const Point<T>& p) {
    bool in = false;
    int n = P.size();
    for (int i = 0, j = 1; i < n; i++, j++) {
        if (j == n) j -= n;
        Point a = P[i] - p, b = P[j] - p;
        if (a.y > b.y) std::swap(a, b);
        if (sgn(a.y) <= 0 and sgn(b.y) > 0 and sgn(a.det(b)) < 0) in = not in;
        if (sgn(a.det(b)) == 0 and sgn(a.dot(b)) <= 0) return Containment::ON;
    }
    return in ? Containment::IN : Containment::OUT;
}

}  // namespace geometry
#line 2 "src/geometry/Point.hpp"
#include <cassert>
#include <cmath>
#include <iostream>
#include <type_traits>

namespace geometry {

template <typename T> struct Point {
    static T EPS;

    static void set_eps(T eps) { EPS = eps; }

    T x, y;

    Point() {}

    Point(T x, T y) : x(x), y(y) {}

    Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }

    Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }

    Point operator*(T t) const { return Point(x * t, y * t); }

    Point operator/(T t) const { return Point(x / t, y / t); }

    bool operator==(const Point& p) const { return x == p.x and y == p.y; }

    bool operator!=(const Point& p) const { return not((*this) == p); }

    bool operator<(const Point& p) const { return x != p.x ? x < p.x : y < p.y; }

    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; }

    T norm() { return std::sqrt(x * x + y * y); }

    T norm2() { return x * x + y * y; }

    T arg() { return std::atan2(y, x); }

    T dot(const Point& p) { return x * p.x + y * p.y; }

    T det(const Point& p) { return x * p.y - y * p.x; }

    Point perp() { return Point(-y, x); }

    Point unit() { return *this / norm(); }

    Point normal() { return perp().unit(); }

    Point rotate(T rad) { return Point(std::cos(rad) * x - std::sin(rad) * y, std::sin(rad) * x + std::cos(rad) * y); }
};

template <> double Point<double>::EPS = 1e-9;
template <> long double Point<long double>::EPS = 1e-12;
template <> int Point<int>::EPS = 0;
template <> long long Point<long long>::EPS = 0;

template <typename T> int sgn(T x) { return x < -Point<T>::EPS ? -1 : x > Point<T>::EPS ? 1 : 0; }

}  // namespace geometry
#line 3 "src/geometry/Circle.hpp"

namespace geometry {

template <typename T> struct Circle {
    Point<T> c;
    T r;

    Circle() {}

    Circle(Point<T> c, T r) : c(c), r(r) {}

    friend std::istream& operator>>(std::istream& is, Circle& c) { return is >> c.c >> c.r; }

    friend std::ostream& operator<<(std::ostream& os, const Circle& c) { return os << c.c << ' ' << c.r; }
};

}  // namespace geometry
#line 2 "src/geometry/Polygon.hpp"
#include <vector>
#line 3 "src/geometry/ccw.hpp"

namespace geometry {

enum Position { COUNTER_CLOCKWISE = +1, CLOCKWISE = -1, ONLINE_BACK = +2, ONLINE_FRONT = -2, ON_SEGMENT = 0 };

template <typename T> Position ccw(const Point<T>& a, Point<T> b, Point<T> c) {
    b = b - a;
    c = c - a;
    if (sgn(b.det(c)) == 1) return COUNTER_CLOCKWISE;
    if (sgn(b.det(c)) == -1) return CLOCKWISE;
    if (sgn(b.dot(c)) == -1) return ONLINE_BACK;
    if (b.norm2() < c.norm2()) return ONLINE_FRONT;
    return ON_SEGMENT;
}

}  // namespace geometry
#line 4 "src/geometry/Polygon.hpp"

namespace geometry {

template <typename T> struct Polygon : std::vector<Point<T>> {
    using std::vector<Point<T>>::vector;

    Polygon(int n) : std::vector<Point<T>>(n) {}

    T area2() {
        T sum = 0;
        int n = this->size();
        for (int i = 0; i < n; i++) sum += (*this)[i].det((*this)[i + 1 == n ? 0 : i + 1]);
        return sum < 0 ? -sum : sum;
    }

    T area() { return area2() / 2; }

    bool is_convex() {
        int n = this->size();
        for (int j = 0; j < n; j++) {
            int i = (j == 0 ? n - 1 : j - 1), k = (j == n - 1 ? 0 : j + 1);
            if (ccw((*this)[i], (*this)[j], (*this)[k]) == CLOCKWISE) return false;
        }
        return true;
    }
};

}  // namespace geometry
#line 2 "src/geometry/distance.hpp"
#include <algorithm>
#line 2 "src/geometry/crosspoint.hpp"
#include <utility>
#line 3 "src/geometry/Line.hpp"

namespace geometry {

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

    Line() {}

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

    // A * x + B * y + C = 0
    Line(T A, T B, T C) {}

    friend std::istream& operator>>(std::istream& is, Line& l) { return is >> l.a >> l.b; }

    friend std::ostream& operator<<(std::ostream& os, const Line& l) { return os << l.a << " to " << l.b; }
};

template <typename T> struct Segment : Line<T> {
    Segment() {}

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

}  // namespace geometry
#line 4 "src/geometry/projection.hpp"

namespace geometry {

template <typename T> Point<T> projection(const Line<T>& l, const Point<T>& p) {
    Point<T> a = p - l.a, b = l.b - l.a;
    return l.a + b * a.dot(b) / b.norm2();
}

}  // namespace geometry
#line 7 "src/geometry/crosspoint.hpp"

namespace geometry {

template <typename T> bool is_parallel(const Line<T>& l, const Line<T>& m) {
    return sgn((l.b - l.a).det(m.b - m.a)) == 0;
}

template <typename T> bool is_orthogonal(const Line<T>& l, const Line<T>& m) {
    return sgn((l.b - l.a).dot(m.b - m.a)) == 0;
}

template <typename T> bool has_crosspoint(const Segment<T>& s, const Segment<T>& t) {
    return ccw(s.a, s.b, t.a) * ccw(s.a, s.b, t.b) <= 0 and ccw(t.a, t.b, s.a) * ccw(t.a, t.b, s.b) <= 0;
}

template <typename T> int count_common_tangent(const Circle<T>& c, const Circle<T>& d) {
    T dist = (c.c - d.c).norm();
    int tmp = sgn(dist - (c.r + d.r));
    if (tmp > 0) return 4;
    if (tmp == 0) return 3;
    tmp = sgn(dist - (sgn(c.r - d.r) > 0 ? c.r - d.r : d.r - c.r));
    if (tmp > 0) return 2;
    if (tmp == 0) return 1;
    return 0;
}

template <typename T> Point<T> crosspoint(const Line<T>& l, const Line<T>& m) {
    assert(not is_parallel(l, m));
    Point<T> u = l.b - l.a, v = m.b - m.a;
    return l.a + u * v.det(m.a - l.a) / v.det(u);
}

template <typename T> Point<T> crosspoint(const Segment<T>& s, const Segment<T>& t) {
    assert(has_crosspoint(s, t));
    if (not is_parallel(s, t)) return crosspoint(Line(s.a, s.b), Line(t.a, t.b));
    std::vector<Point<T>> v = {s.a, s.b, t.a, t.b};
    for (int i = 0; i <= 2; i++) {
        for (int j = 2; j >= i; j--) {
            if (v[j + 1] < v[j]) {
                std::swap(v[j], v[j + 1]);
            }
        }
    }
    return v[1];
}

template <typename T> std::vector<Point<T>> crosspoint(const Circle<T>& c, const Line<T>& l) {
    Point<T> h = projection(l, c.c);
    T x = c.r * c.r - (c.c - h).norm2();
    if (sgn(x) < 0) return {};
    if (sgn(x) == 0) return {h};
    Point<T> v = (l.b - l.a).unit() * std::sqrt(x);
    return {h - v, h + v};
}

template <typename T> std::vector<Point<T>> crosspoint(const Circle<T>& c, const Segment<T>& s) {}

template <typename T> std::vector<Point<T>> crosspoint(const Circle<T>& c1, const Circle<T>& c2) {
    T r1 = c1.r, r2 = c2.r;
    if (r1 < r2) return crosspoint(c2, c1);
    T d = (c2.c - c1.c).norm();
    if (sgn(d - (r1 + r2)) > 0 or sgn(d - (r1 - r2)) < 0) return {};
    Point<T> v = c2.c - c1.c;
    if (sgn(d - (r1 + r2)) == 0 or sgn(d - (r1 - r2)) == 0) return {c1.c + v.unit() * r1};
    T p = ((r1 * r1 - r2 * r2) / d + d) / 2, q = std::sqrt(r1 * r1 - p * p);
    Point<T> h = c1.c + v.unit() * p;
    Point<T> i = v.normal();
    return {h + i * q, h - i * q};
}

}  // namespace geometry
#line 5 "src/geometry/distance.hpp"

namespace geometry {

template <typename T> T distance(const Point<T>& p, const Point<T>& q) { return (p - q).norm(); }

template <typename T> T distance(const Line<T>& l, const Point<T>& p) { return distance(p, projection(l, p)); }

template <typename T> T distance(const Segment<T>& s, const Point<T>& p) {
    Point<T> h = projection(s, p);
    return ccw(s.a, s.b, h) == ON_SEGMENT ? distance(p, h) : std::min(distance(p, s.a), distance(p, s.b));
}

template <typename T> T distance(const Segment<T>& s, const Segment<T>& t) {
    return has_crosspoint(s, t) ? 0
                                : std::min({distance(s, t.a), distance(s, t.b), distance(t, s.a), distance(t, s.b)});
}

}  // namespace geometry
#line 5 "src/geometry/contain.hpp"

namespace geometry {

enum class Containment { OUT, ON, IN };

template <typename T> Containment contain(const Circle<T>& c, const Point<T>& p) {
    int x = sgn(distance(c.c, p) - c.r);
    return x == -1 ? Containment::IN : x == 1 ? Containment::OUT : Containment::ON;
}

template <typename T> Containment contain(const Polygon<T>& P, const Point<T>& p) {
    bool in = false;
    int n = P.size();
    for (int i = 0, j = 1; i < n; i++, j++) {
        if (j == n) j -= n;
        Point a = P[i] - p, b = P[j] - p;
        if (a.y > b.y) std::swap(a, b);
        if (sgn(a.y) <= 0 and sgn(b.y) > 0 and sgn(a.det(b)) < 0) in = not in;
        if (sgn(a.det(b)) == 0 and sgn(a.dot(b)) <= 0) return Containment::ON;
    }
    return in ? Containment::IN : Containment::OUT;
}

}  // namespace geometry
Back to top page