ferinの競プロ帳

競プロについてのメモ

Educational Codeforces Round 2 D. Area of Two Circles' Intersection

問題ページ

まず2円が離れているか外接の関係にある場合、答えは0です。2円が内包か内接の関係にある場合、答えは \pi \min(r_1, r_2)^2 です。2円が2点で交差する場合について、2通りに場合分けを行います。

  • 交点からなる線分から見て、(x_1,y_1),(x_2,y_2) が反対側にある
    扇形ABCから三角形ABCを引くことで、各円について斜線部の領域を求めます。斜線部の和が答えの面積になります。 f:id:ferin_tech:20200315220846p:plainf:id:ferin_tech:20200315220850p:plain
  • 交点からなる線分から見て、(x_1,y_1),(x_2,y_2) が同じ側にある
    半径が大きい方の円については先程と同様に、扇形から三角形を引きます。半径が小さい方については扇形の面積に三角形を"足す"ことで該当する領域の面積を求められます。
    f:id:ferin_tech:20200315220857p:plainf:id:ferin_tech:20200315220906p:plain

どちらかの円の中心が交点からなる線分上に乗っている場合もありますが、その場合足し引きする三角形の面積が0なのでどちらでも大丈夫です。どちらの場合にあたるか?を調べるにはccw関数などを用いればよいです。

扇形と三角形の面積を求めます。扇形については r^2 \angle BAC / 2、三角形については r^2 \sin \angle BAC / 2 で求められます。角度については内積を求める式 \vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}|\cos \angle BAC より \angle BAC = \arccos (\vec{a} \cdot \vec{b}) / |\vec{a}||\vec{b}| と求められます。

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG_   
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
const ll INF = 1LL<<60;  
  
using R = long double; // Rにmint渡せるようにする  
using P = complex<R>;  
struct C {  
    P c; R r;  
    C() {}  
    C(const P &a, const R &b) : c(a), r(b) {}  
};  
  
const R EPS = 1e-8;  
const R PI = atanl(1)*4;  
  
inline int sgn(const R& r) { return (r>EPS) - (r<-EPS); }  
inline R dot(const P& a, const P& b) {  
    return real(a)*real(b) + imag(a)*imag(b);  
}  
inline R det(const P& a, const P& b) {  
    return real(a)*imag(b) - imag(a)*real(b);  
}  
namespace std {  
    bool operator < (const P& a, const P& b) {  
        return sgn(real(a-b)) ? real(a-b) < 0 : sgn(imag(a-b)) < 0;  
    }  
    bool operator == (const P& a, const P& b) {  
        return sgn(real(a-b)) == 0 && sgn(imag(a-b)) == 0;  
    }  
    bool cmp_y (const P& a, const P& b) {  
        return sgn(imag(a-b)) ? imag(a-b) < 0 : sgn(real(a-b)) < 0;  
    }  
}  
  
// 線分abから見たcの位置  
enum CCW{LEFT=1, RIGHT=2, BACK=4, FRONT=8, ON_SEG=16};  
int ccw(P a, P b, P c) {  
    P p = (c-a)/(b-a);  
    if(sgn(imag(p)) > 0) return LEFT;  
    if(sgn(imag(p)) < 0) return RIGHT;  
    if(sgn(real(p)) < 0) return BACK;  
    if(sgn(real(p)-1) > 0) return FRONT;  
    return ON_SEG;  
}  
  
// 交差  
int intersect(const C& a, const C& b) {  
    R dist = norm(a.c-b.c), r1 = (a.r+b.r)*(a.r+b.r), r2 = (a.r-b.r)*(a.r-b.r);  
    if(sgn(r1-dist) < 0)  return 4;    // 円が離れている  
    if(sgn(r1-dist) == 0) return 3;   // 外接  
    if(sgn(r2-dist) < 0 && sgn(dist-r1) < 0) return 2; // 交差  
    if(sgn(dist-r2) == 0) return 1; // 内接  
    return 0;    // 内部に含む  
}  
  
// 交点 交差の確認を先にすること!!!  
// intersect=2 を確認すること  
vector<P> crosspoint(C a, C b) {  
    R d = abs(a.c-b.c);  
    R t = (a.r*a.r-b.r*b.r+d*d)/2/d, h = sqrt(a.r*a.r-t*t);  
    P m = t/abs(b.c-a.c)*(b.c-a.c)+a.c;  
    P n(-(a.c-b.c).imag()/abs(a.c-b.c), (a.c-b.c).real()/abs(a.c-b.c));  
    vector<P> ret(2, m);  
    ret[0] -= h*n; ret[1] += h*n;  
    return ret;  
}  
  
// 2円の共通面積を求める  
R common_area(C c1, C c2) {  
    ll type = intersect(c1, c2);  
    if(type >= 3) return 0;  
    if(type <= 1) {  
        ll r = min(c1.r, c2.r);  
        return PI*r*r;  
    }  
    if(c1.r > c2.r) swap(c1, c2);  
    vector<P> ps = crosspoint(c1, c2);  
    R ang1 = acosl(dot(ps[0]-c1.c, ps[1]-c1.c) / (c1.r * c1.r));  
    R ang2 = acosl(dot(ps[0]-c2.c, ps[1]-c2.c) / (c2.r * c2.r));  
    if(ccw(ps[0], ps[1], c1.c) == ccw(ps[0], ps[1], c2.c)) ang1 = max(ang1, 2*PI - ang1);  
    else ang1 = min(ang1, 2*PI - ang1);  
    ang2 = min(ang2, 2*PI-ang2);  
    return c1.r*c1.r*0.5*(ang1-sinl(ang1)) + c2.r*c2.r*0.5*(ang2-sinl(ang2));  
}  
  
int main(void) {  
    ll x1, y1, r1, x2, y2, r2;  
    cin >> x1 >> y1 >> r1 >> x2 >> y2 >> r2;  
    C c1({{x1, y1}, r1}), c2({{x2, y2}, r2});  
  
    cout << fixed << setprecision(9) << common_area(c1, c2) << endl;  
  
    return 0;  
}