rcpl

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

View the Project on GitHub ruthen71/rcpl

:heavy_check_mark: verify/geometry/common_area_cp.test.cpp

Depends on

Code

#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_H"
#define ERROR 0.00001

#include <iostream>
#include <iomanip>

#include "geometry/common_area.hpp"

int main() {
    int n;
    double r;
    std::cin >> n >> r;
    Polygon<double> p(n);
    std::cin >> p;
    // 適当に平行移動しておく
    Circle<double> c(Point<double>(100, 100), r);
    for (auto&& pt : p) pt.x += 100, pt.y += 100;
    std::cout << std::fixed << std::setprecision(15) << common_area(c, p) << '\n';
    return 0;
}
#line 1 "verify/geometry/common_area_cp.test.cpp"
#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_H"
#define ERROR 0.00001

#include <iostream>
#include <iomanip>

#line 2 "geometry/common_area.hpp"

#line 2 "geometry/circle.hpp"

#line 2 "geometry/point.hpp"

#line 2 "geometry/geometry_template.hpp"

#include <type_traits>

// Constants (EPS, PI)
// EPS の変更は Constants<T>::set_eps(new_eps) で
template <class T> struct Constants {
    static T EPS;
    static void set_eps(const T e) { EPS = e; }
    static constexpr T PI = 3.14159'26535'89793L;
};

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

// 汎用関数
template <class T> inline int sign(const T x) { return x < -Constants<T>::EPS ? -1 : (x > Constants<T>::EPS ? 1 : 0); }
template <class T> inline bool equal(const T a, const T b) { return sign(a - b) == 0; }
template <class T> inline T radian_to_degree(const T r) { return r * 180.0 / Constants<T>::PI; }
template <class T> inline T degree_to_radian(const T d) { return d * Constants<T>::PI / 180.0; }

// type traits
template <class T> using is_geometry_floating_point = typename std::conditional<std::is_same<T, double>::value || std::is_same<T, long double>::value, std::true_type, std::false_type>::type;
template <class T> using is_geometry_integer = typename std::conditional<std::is_same<T, long long>::value, std::true_type, std::false_type>::type;
template <class T> using is_geometry = typename std::conditional<is_geometry_floating_point<T>::value || is_geometry_integer<T>::value, std::true_type, std::false_type>::type;
#line 4 "geometry/point.hpp"

#include <cmath>
#include <cassert>

// 点
template <class T> struct Point {
    T x, y;

    Point() = default;
    Point(const T x, const T y) : x(x), y(y) {}
    template <class U> Point(const Point<U> p) : x(p.x), y(p.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) {
        static_assert(is_geometry_floating_point<T>::value == true);
        return *this = Point(x * p.x - y * p.y, x * p.y + y * p.x);
    }
    Point& operator/=(const Point& p) {
        static_assert(is_geometry_floating_point<T>::value == true);
        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 T k) {
        static_assert(is_geometry_floating_point<T>::value == true);
        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& 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& p, const T k) { return Point(p) /= k; }

    // 辞書式順序
    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; }
    friend bool operator==(const Point& a, const Point& b) { return a.x == b.x and a.y == b.y; }

    // 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 <class T> inline bool equal(const Point<T>& a, const Point<T>& b) { return equal(a.x, b.x) and equal(a.y, b.y); }
// 内積
template <class T> inline T dot(const Point<T>& a, const Point<T>& b) { return a.x * b.x + a.y * b.y; }
// 外積
template <class T> inline T cross(const Point<T>& a, const Point<T>& b) { return a.x * b.y - a.y * b.x; }
// rad ラジアンだけ反時計回りに回転
template <class T> inline Point<T> rotate(const Point<T>& p, const T theta) {
    static_assert(is_geometry_floating_point<T>::value == true);
    return p * Point<T>(std::cos(theta), std::sin(theta));
}
// (x, y) の辞書式順序 (誤差許容)
template <class 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; }
// (y, x) の辞書式順序 (誤差許容)
template <class 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; }
// 整数のまま行う偏角ソート
// 無限の精度をもつ arg(p) = atan2(y, x) で比較し, 同じ場合は norm(p) で比較 (atan2(0, 0) = 0 とする)
// 基本的に (-PI, PI] でソートされ, 点 (0, 0) は (-PI, 0) と [0, PI] の間に入る
// https://ngtkana.hatenablog.com/entry/2021/11/13/202103
// https://judge.yosupo.jp/problem/sort_points_by_argument
template <class T> inline bool compare_atan2(const Point<T>& a, const Point<T>& b) {
    static_assert(is_geometry_integer<T>::value == true);
    if ((Point<T>(a.y, -a.x) > Point<T>(0, 0)) == (Point<T>(b.y, -b.x) > Point<T>(0, 0))) {  // a, b in (-PI, 0] or a, b in (0, PI]
        if (a.x * b.y != a.y * b.x) return a.x * b.y > a.y * b.x;                            // cross(a, b) != 0
        return a == Point<T>(0, 0) ? b.x > 0 and b.y == 0 : (b == Point<T>(0, 0) ? a.y < 0 : norm(a) < norm(b));
    }
    return Point<T>(a.y, -a.x) < Point<T>(b.y, -b.x);
}
// 絶対値の 2 乗
template <class T> inline T norm(const Point<T>& p) { return p.x * p.x + p.y * p.y; }
// 絶対値
template <class T> inline T abs(const Point<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    return std::sqrt(norm(p));
}
// 偏角
template <class T> inline T arg(const Point<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    return std::atan2(p.y, p.x);  // (-PI, PI]
}
// 共役複素数 (x 軸について対象な点)
template <class T> inline Point<T> conj(const Point<T>& p) { return Point(p.x, -p.y); }
// 極座標系 -> 直交座標系 (絶対値が rho で偏角が theta ラジアン)
template <class T> inline Point<T> polar(const T rho, const T theta = T(0)) {
    static_assert(is_geometry_floating_point<T>::value == true);
    assert(sign(rho) >= 0);
    return Point<T>(std::cos(theta), std::sin(theta)) * rho;
}
// ccw の戻り値
enum class Ccw {
    COUNTER_CLOCKWISE = 1,  // a->b->c 反時計回り
    CLOCKWISE = -1,         // a->b->c 時計回り
    ONLINE_BACK = 2,        // c->a->b 直線
    ONLINE_FRONT = -2,      // a->b->c 直線
    ON_SEGMENT = 0          // a->c->b 直線
};
// 点 a, b, c の位置関係
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C
template <class T> Ccw ccw(const Point<T>& a, const Point<T>& b, const Point<T>& c) {
    Point<T> ab = b - a;
    Point<T> ac = c - a;
    if (sign(cross(ab, ac)) == 1) return Ccw::COUNTER_CLOCKWISE;
    if (sign(cross(ab, ac)) == -1) return Ccw::CLOCKWISE;
    if (sign(dot(ab, ac)) == -1) return Ccw::ONLINE_BACK;
    if (sign(norm(ab) - norm(ac)) == -1) return Ccw::ONLINE_FRONT;
    return Ccw::ON_SEGMENT;
}
// 線分 a -> b から 線分 a -> c までの角度 (ラジアンで -PI より大きく PI 以下)
template <class T> T get_angle(const Point<T>& a, const Point<T>& b, const Point<T>& c) {
    Point<T> ab = b - a;
    Point<T> ac = c - a;
    // a-bが x 軸になるように回転
    ac *= conj(ab) / norm(ab);
    return arg(ac);  // (-PI, PI]
}
#line 4 "geometry/circle.hpp"

// 円
template <class 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; }
    friend std::ostream& operator<<(std::ostream& os, const Circle& c) { return os << c.o << ", " << c.r; }
};

// 共通接線の本数
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_A
template <class T> int tangent_number(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 ^ 2 > (c1.r + c2.r) ^ 2
    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 ^ 2 = (c1.r + c2.r) ^ 2
    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 ^ 2 > (c1.r - c2.r) ^ 2
    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 ^ 2 = (c1.r - c2.r) ^ 2
    return 0;
}
#line 2 "geometry/line.hpp"

#line 4 "geometry/line.hpp"

// 直線
template <class 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) {
        static_assert(is_geometry_floating_point<T>::value == true);
        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), -A / B);
        } 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; }
};

// 線分
template <class T> struct Segment : Line<T> {
    Segment() = default;

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

// 点 p から直線 l に下ろした垂線と直線 l の交点
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_A
template <class T> Point<T> projection(const Line<T>& l, const Point<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    T t = dot(p - l.a, l.b - l.a) / norm(l.b - l.a);
    return l.a + (l.b - l.a) * t;
}

// 直線 l に関して点 p と対象な位置にある点
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_B
template <class T> Point<T> reflection(const Line<T>& l, const Point<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    return p + (projection(l, p) - p) * T(2);
}

// 直線 l1, l2 の垂直判定
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <class 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; }

// 直線 l1, l2 の平行判定
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_A
template <class 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 2 "geometry/polygon.hpp"

#line 4 "geometry/polygon.hpp"

#include <vector>
#line 7 "geometry/polygon.hpp"

// 多角形
template <class T> using Polygon = std::vector<Point<T>>;
template <class T> std::istream& operator>>(std::istream& is, Polygon<T>& p) {
    for (auto&& pi : p) is >> pi;
    return is;
}
template <class T> std::ostream& operator<<(std::ostream& os, const Polygon<T>& p) {
    for (auto&& pi : p) os << pi << " -> ";
    return os;
}

// 多角形の面積 (符号付き)
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_A
// return area * 2
template <class T> T area2(const Polygon<T>& p) {
    const int n = (int)(p.size());
    assert(n >= 2);
    T res = T(0);
    for (int i = 0; i < n; i++) res += cross(p[i], p[i + 1 == n ? 0 : i + 1]);
    // counter clockwise: res > 0
    // clockwise: res < 0
    return res;
}
template <class T> T area(const Polygon<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    return area2(p) / T(2);
}

// 多角形の凸判定 (角度が 0 でも PI でも許容)
// 許容したくないときには ON_SEGMENT, ONLINE_FRONT, ONLINE_BACK が出てきたら false を返せば OK
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_B
template <class T> bool is_convex(const Polygon<T>& p) {
    const int n = (int)(p.size());
    assert(n >= 3);
    bool okccw = true, okcw = true;
    for (int i = 0; i < n; i++) {
        auto res = ccw(p[i], p[i + 1 >= n ? i + 1 - n : i + 1], p[i + 2 >= n ? i + 2 - n : i + 2]);
        if (res == Ccw::CLOCKWISE) okccw = false;
        if (res == Ccw::COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    return true;
}
#line 2 "geometry/cross_point.hpp"

#line 2 "geometry/is_intersect.hpp"

#line 6 "geometry/is_intersect.hpp"

// 交差判定 (直線, 線分, 円, 点)
// 直線 l1, l2 の交差判定
template <class T> bool is_intersect(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);
    // sign(d12) != 0 -> 平行でないので交差する
    // sign(d12) == 0 and sign(d1) == 0 -> 一致するので交差する
    return sign(d12) != 0 or sign(d1) == 0;
}
// 直線 l, 点 p の交差判定
template <class T> inline bool is_intersect(const Line<T>& l, const Point<T>& p) {
    auto res = ccw(l.a, l.b, p);
    return res == Ccw::ONLINE_BACK or res == Ccw::ONLINE_FRONT or res == Ccw::ON_SEGMENT;
}
template <class T> bool is_intersect(const Point<T>& p, const Line<T>& l) { return is_intersect(l, p); }

// 線分 s, 点 p の交差判定
template <class T> inline bool is_intersect(const Segment<T>& s, const Point<T>& p) { return ccw(s.a, s.b, p) == Ccw::ON_SEGMENT; }
template <class T> inline bool is_intersect(const Point<T>& p, const Segment<T>& s) { return ccw(s.a, s.b, p) == Ccw::ON_SEGMENT; }

// 直線 l, 線分 s の交差判定
template <class T> bool is_intersect(const Line<T>& l, const Segment<T>& s) {
    // s.a と s.b が直線 l に関して同じ側 (直線上を除き) にある場合に限り交差しない
    auto c1 = ccw(l.a, l.b, s.a);
    auto c2 = ccw(l.a, l.b, s.b);
    return !((c1 == c2) and (c1 == Ccw::CLOCKWISE or c1 == Ccw::COUNTER_CLOCKWISE));
}
template <class T> bool is_intersect(const Segment<T>& s, const Line<T>& l) { return is_intersect(l, s); }

// 線分 s1, s2 の交差判定
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_B
template <class T> inline bool is_intersect(const Segment<T>& s1, const Segment<T>& s2) {
    auto c1 = ccw(s1.a, s1.b, s2.a);
    auto c2 = ccw(s1.a, s1.b, s2.b);
    auto c3 = ccw(s2.a, s2.b, s1.a);
    auto c4 = ccw(s2.a, s2.b, s1.b);
    // 平行な場合: 重なるなら 1 次元で考えると必ず端点のどれかがもう一方の線分上にある
    if (c1 == Ccw::ON_SEGMENT or c2 == Ccw::ON_SEGMENT or c3 == Ccw::ON_SEGMENT or c4 == Ccw::ON_SEGMENT) return true;
    // 平行でない場合: 一方の線分の両側にもう一方の線分の端点がある, を線分を入れ替えても成立
    // 平行だが重ならない場合は以下の条件は成立しない
    bool ok1 = ((c1 == Ccw::CLOCKWISE and c2 == Ccw::COUNTER_CLOCKWISE) or (c1 == Ccw::COUNTER_CLOCKWISE and c2 == Ccw::CLOCKWISE));
    bool ok2 = ((c3 == Ccw::CLOCKWISE and c4 == Ccw::COUNTER_CLOCKWISE) or (c3 == Ccw::COUNTER_CLOCKWISE and c4 == Ccw::CLOCKWISE));
    return ok1 and ok2;
}

// 点 p1, p2 の交差判定
template <class T> inline bool is_intersect(const Point<T>& p1, const Point<T>& p2) { return equal(p1, p2); }

// 円 c1, c2 の交差判定
template <class T> inline bool is_intersect(const Circle<T>& c1, const Circle<T>& c2) {
    const int num = tangent_number(c1, c2);
    return 1 <= num and num <= 3;
}

// 円 c, 点 p の交差判定
template <class T> inline bool is_intersect(const Circle<T>& c, const Point<T>& p) { return equal(norm(p - c.o), c.r * c.r); }
template <class T> inline bool is_intersect(const Point<T>& p, const Circle<T>& c) { return equal(norm(p - c.o), c.r * c.r); }

// 円 c, 直線 l の交差判定
template <class T> inline bool is_intersect(const Circle<T>& c, const Line<T>& l) {
    static_assert(is_geometry_floating_point<T>::value == true);
    // norm(c.o - projection(l, c.o)) が直線 l と 円 c の中心 c.o の距離の 2 乗
    return sign(c.r * c.r - norm(c.o - projection(l, c.o))) >= 0;
}
template <class T> inline bool is_intersect(const Line<T>& l, const Circle<T>& c) { return is_intersect(c, l); }

// 円 c, 線分 s の交差判定
template <class T> inline bool is_intersect(const Circle<T>& c, const Segment<T>& s) {
    static_assert(is_geometry_floating_point<T>::value == true);
    if (!is_intersect(c, Line(s.a, s.b))) return false;  // 直線としても交差しない
    T d1 = abs(c.o - s.a), d2 = abs(c.o - s.b);
    if (sign(d1 - c.r) == -1 and sign(d2 - c.r) == -1) return false;  // 端点がどちらも円の内側
    if (sign(d1 - c.r) * sign(d2 - c.r) <= 0) return true;            // 円周をまたいでいる
    const Point<T> h = projection(s, c.o);
    return ccw(s.a, s.b, h) == Ccw::ON_SEGMENT;  // s.a -> h -> s.b の順で並んでいる
}
template <class T> inline bool is_intersect(const Segment<T>& s, const Circle<T>& c) { return is_intersect(c, s); }
#line 7 "geometry/cross_point.hpp"

#line 10 "geometry/cross_point.hpp"

// 交点 (直線, 線分, 円)
// 交点を持たない場合の挙動は未定義
// 直線 l1, l2 の交点 1 つ
template <class T> Point<T> cross_point(const Line<T>& l1, const Line<T>& l2) {
    static_assert(is_geometry_floating_point<T>::value == true);
    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) {
        assert(sign(d1) == 0);  // 平行かつ一致しない
        return l2.a;
    }
    return l2.a + (l2.b - l2.a) * (d1 / d12);
}

// 線分 s1, s2 の交点 1 つ
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_2_C
template <class T> Point<T> cross_point(const Segment<T>& s1, const Segment<T>& s2) {
    static_assert(is_geometry_floating_point<T>::value == true);
    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) {
        assert(sign(d1) == 0);  // 平行かつ一致しない
        // 端点のどれかを返す
        if (is_intersect(s1, s2.a)) return s2.a;
        if (is_intersect(s1, s2.b)) return s2.b;
        if (is_intersect(s2, s1.a)) return s1.a;
        assert(is_intersect(s2, s1.b));
        return s1.b;
    }
    return s2.a + (s2.b - s2.a) * (d1 / d12);
}

// 直線 l, 線分 s の交点 1 つ
template <class T> Point<T> cross_point(const Line<T>& l, const Segment<T>& s) {
    // cross_point(l1, l2) は重なるとき l2.a を返すので OK
    return cross_point(l, Line<T>(s.a, s.b));
}
template <class T> Point<T> cross_point(const Segment<T>& s, const Line<T>& l) { return cross_point(l, s); }

// 円 c, 直線 l の交点 1 or 2 つ
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_D
template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c, const Line<T>& l) {
    static_assert(is_geometry_floating_point<T>::value == true);
    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) / abs(l.b - l.a);
    auto k = sqrtl(c.r * c.r - norm(pr - c.o));
    return {pr - e * k, pr + e * k};
}
template <class T> std::vector<Point<T>> cross_point(const Line<T>& l, const Circle<T>& c) { return cross_point(c, l); }

// 円 c, 線分 s の交点 1 or 2 つ
template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c, const Segment<T>& s) {
    assert(is_intersect(c, s));
    auto ps = cross_point(c, Line<T>(s.a, s.b));
    std::vector<Point<T>> res;
    for (auto&& p : ps)
        if (ccw(s.a, s.b, p) == Ccw::ON_SEGMENT) res.emplace_back(p);
    return res;
}
template <class T> std::vector<Point<T>> cross_point(const Segment<T>& s, const Circle<T>& c) { return cross_point(c, s); }

// 円 c1, c2 の交点 1 or 2 つ
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_E
template <class T> std::vector<Point<T>> cross_point(const Circle<T>& c1, const Circle<T>& c2) {
    static_assert(is_geometry_floating_point<T>::value == true);
    assert(is_intersect(c1, c2));
    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, q)) return {p};
    return {p, q};
}
#line 7 "geometry/common_area.hpp"

// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_7_H
// https://drken1215.hatenablog.com/entry/2020/02/02/091000
// 円 c の中心 o, 線分 s のなす三角形と円 c の共通部分の面積
template <class T> T common_area(Circle<T> c, Segment<T> s) {
    // c.o を原点と見て平行移動
    s.a -= c.o;
    s.b -= c.o;
    c.o -= c.o;
    // c.o, s.a, s.b が一直線に並ぶ場合
    if (sign(norm(s.a - s.b)) == 0) return 0;
    bool a_in_c = sign(norm(s.a) - c.r * c.r) == -1;
    bool b_in_c = sign(norm(s.b) - c.r * c.r) == -1;
    // s.a, s.b がともに円 c の内部にある場合 -> 三角形
    if (a_in_c and b_in_c) return cross(s.a, s.b) / 2;
    // s.a, s.b がともに円 c の外部にある場合 -> 扇形
    if (!is_intersect(c, s)) return c.r * c.r * get_angle(c.o, s.a, s.b) / 2;
    auto cp = cross_point(c, s);
    // 交点を cp1, cp2 とする (1 つの場合は cp1 == cp2)
    auto cp1 = cp.front(), cp2 = cp.back();
    // 線分 s.a -> cp1, cp1 -> cp2, cp2 -> s.b について, 扇形か三角形か場合分けしながら求めていく
    T res = 0;
    // s.a -> cp1
    if (a_in_c) {
        res += cross(s.a, cp1) / 2;
    } else {
        res += c.r * c.r * get_angle(c.o, s.a, cp1) / 2;
    }
    // cp1 -> cp2
    res += cross(cp1, cp2) / 2;
    // cp2 -> s.b
    if (b_in_c) {
        res += cross(cp2, s.b) / 2;
    } else {
        res += c.r * c.r * get_angle(c.o, cp2, s.b) / 2;
    }
    return res;
}

// 円 c, 多角形 p の共通部分の面積 (符号付き)
template <class T> T common_area(const Circle<T>& c, const Polygon<T>& p) {
    static_assert(is_geometry_floating_point<T>::value == true);
    const int n = (int)(p.size());
    assert(n >= 2);
    T res = T(0);
    for (int i = 0; i < n; i++) res += common_area(c, Segment(p[i], p[(i + 1) % n]));
    // counter clockwise: res > 0
    // clockwise: res < 0
    return res;
}
template <class T> T common_area(const Polygon<T>& p, const Circle<T>& c) { return common_area(c, p); }

// 円 c1, c2 の共通部分の面積
// 扇形 2 つの面積の和からひし形の面積を引く
template <class T> T common_area(const Circle<T>& c1, const Circle<T>& c2) {
    static_assert(is_geometry_floating_point<T>::value == true);
    const int num = tangent_number(c1, c2);
    if (num >= 3) return 0;
    if (num <= 1) {
        // 一方に他方が完全に含まれる
        T r = std::min(c1.r, c2.r);
        return r * r * Constants<T>::PI;
    }
    auto cp = cross_point(c1, c2);
    T res = T(0);
    // get_angle(c1.o, cp[0], cp[1]) にすると扇形の中心角が直角より大きいときにダメ
    // 円 c1 を含む扇形
    res += c1.r * c1.r * std::abs(get_angle(c1.o, c2.o, cp[0]));
    // 円 c2 を含む扇形
    res += c2.r * c2.r * std::abs(get_angle(c2.o, c1.o, cp[0]));
    // ひし形
    res -= abs(cp[1] - cp[0]) * abs(c1.o - c2.o) / 2;
    return res;
}
#line 8 "verify/geometry/common_area_cp.test.cpp"

int main() {
    int n;
    double r;
    std::cin >> n >> r;
    Polygon<double> p(n);
    std::cin >> p;
    // 適当に平行移動しておく
    Circle<double> c(Point<double>(100, 100), r);
    for (auto&& pt : p) pt.x += 100, pt.y += 100;
    std::cout << std::fixed << std::setprecision(15) << common_area(c, p) << '\n';
    return 0;
}
Back to top page