Geometry: Circle
(geometry/circle.hpp)

View this file on GitHub
Last update: 2026-10-22 23:24:07+09:00
Include: #include "geometry/circle.hpp"

概要

円。浮動小数点型のみ。

使い方

geo::circle<T>(c, r): 中心 c、半径 r の円。
area(o): 面積を返す。
contains(o, p): 点 p が円の内部なら $2$、円周上なら $1$、外部なら $0$ を返す。
intersect(o, l), intersect(o, s), intersect(a, b): 直線、線分、円と共有点があるかを返す。
count_common_tangents(a, b): 円 a, b の共通接線の本数を返す。
cross_points(o, l), cross_points(o, s), cross_points(a, b): 直線、線分、円との交点を x, y の辞書順で返す。接するときは同じ点を 2 つ返す。
tangent_points(o, p): 点 p から円 o に引いた接線の接点を返す。
common_tangents(a, b): 円 a, b の共通接線の、円 a 上の接点を返す。
incircle(a, b, c), circumcircle(a, b, c): 三角形 abc の内接円、外接円を返す。
smallest_enclosing_circle(ps): 点集合 ps を含む最小の円を返す。
area_of_intersection(o, poly), area_of_intersection(a, b): 円と多角形、円と円の共通部分の面積を返す。
circle_of_apollonius(a, b, n, m): 点 a, b からの距離の比が $n : m$ である点の軌跡を返す。

制約

tangent_points(o, p) では p は円の外部または円周上
circle_of_apollonius では $n \ne m$

計算量

smallest_enclosing_circle: 期待 $\mathcal{O}(N)$
area_of_intersection(o, poly): $\mathcal{O}(N)$
その他: $\mathcal{O}(1)$

Depends on

Verified with

Code

#pragma once
#include <algorithm>
#include <cassert>
#include <cmath>
#include <concepts>
#include <numbers>
#include <random>
#include <vector>

#include "geometry/line.hpp"
#include "geometry/point.hpp"
#include "geometry/polygon.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 2 "geometry/circle.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <concepts>
#include <numbers>
#include <random>
#include <vector>

#line 6 "geometry/line.hpp"
#include <iostream>

#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

Geometry: Circle (geometry/circle.hpp)

概要

使い方

制約

計算量

Depends on

Verified with

Code

Geometry: Circle
(geometry/circle.hpp)