Problem: Find the area of intersection between a circle $c$ and a polygon $g$. The center coordinate of the circle is ($0, 0$).

The polygon $g$ is represented by a sequence of points $p_1$, $p_2$,..., $p_n$ where line segments connecting $p_i$ and $p_{i+1}$ ($1 \leq i \leq n−1$) are sides of the polygon. The line segment connecting $p_n$ and $p_1$ is also a side of the polygon.

Note that the polygon is not necessarily convex.

Constraints

* $3 \leq n \leq 100$
* $1 \leq r \leq 100$
* $-100 \leq x_i, y_i \leq 100$

Input

The input is given in the following format.

$n$ $r$
$x_1$ $y_1$
$x_2$ $y_2$
:
$x_n$ $y_n$


In the first line, an integer n representing the number of points in the polygon is given. The coordinate of a point $p_i$ is given by two integers $x_i$ and $y_i$. The coordinates of the points are given in the order of counter-clockwise visit of them. All input values are given in integers.

Output

Print the area of intersection in a line. The output values should be in a decimal fraction with an error less than 0.00001.

Examples

Input

3 5
1 1
4 1
5 5


Output

4.639858417607


Input

4 5
0 0
-3 -6
1 -3
5 -4


Output

11.787686807576

Solution (in CPP):

Solution:
#include <iostream>
#include <vector>
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const int N = 200000 + 10;
const double EPS = 1e-7;
const double INF = 1e12;
int sgn(double x) {
    if (x < -EPS) return -1;
    if (x > +EPS) return +1;
    return 0;
}
struct Vec {
    double x,y;
    Vec() {}
    Vec(double x, double y) { 
        this->x = x, this->y = y; 
    }
    void read() { 
        scanf("%lf %lf", &x, &y); 
    }
    void prt() {
        printf("%.8f %.8f\n", x, y); 
    }
    double len() { 
        return sqrt(x * x + y * y); 
    }
    double len2() { 
        return x * x + y * y; 
    }
    bool operator <(const Vec& o)const {
        if(x == o.x)
            return y < o.y;
        return x < o.x;
    }
    Vec operator +(const Vec& o)const { 
        return Vec(x + o.x, y + o.y); 
    }
    Vec operator -(const Vec& o)const { 
        return Vec(x - o.x, y - o.y); 
    }
    Vec operator *(const double& k)const { 
        return Vec(k * x, k * y); 
    }
    double operator *(const Vec& o)const {
        return x * o.x + y * o.y; 
    }
    double operator ^(const Vec& o)const {
        return x * o.y - y * o.x; 
    }
    Vec rotate(double ang){ 
        return Vec(x * cos(ang) - y * sin(ang), x * sin(ang) + y * cos(ang)); 
    }
    Vec change(double l) { 
        if(len() < EPS) return *this; 
        return (*this) * (l/len()); 
    }
};
bool cmpang(Vec A, Vec B) {
    if ((A ^ B) == 0) {
        return A < B;
    }
    return (A ^ B) > 0;
}
double rad(Vec A, Vec B) {
    if (abs(A^B) < EPS) return 0;
    int sgn = (A^B)>0 ? 1: -1;
    return acos( (A * B) / A.len() / B.len() ) * sgn;
}

struct Line {
    Vec A1, A2;
    Line() {}
    Line(Vec A1, Vec A2) {
        this->A1 = A1, this->A2 = A2;
    }
    double len() {
        return (A2-A1).len();
    }
    double Len2() {
        return (A2-A1).len2();
    }
    double getAng() const {
        return atan2(A2.y-A1.y, A2.x-A1.x);
    }
    void prt() {
        printf("(%.5f %.5f) - (%.5f %.5f) %.7f\n", A1.x,A1.y,A2.x,A2.y,(*this).getAng());
    }
    bool operator <(const Line& o)const {
        Line l = o;
        double A = (*this).getAng();
        double B = l.getAng();
        if (abs(A-B) < EPS)
            return ((A2-A1)^(l.A1-A1)) > EPS;
        return A < B;
    }   
};
// P to L
double disToLine(Vec P, Line L) {
    return abs((P-L.A1)^(L.A2-L.A1)) / L.len();
}
Vec projection(Vec P, Line L) {
    double shadowLen = (P - L.A1) * (L.A2 - L.A1) / (L.A2 - L.A1).len();
    return L.A1 + (L.A2-L.A1).change(shadowLen);
}
Vec rejection(Vec P, Line L) {
    return projection(P,L) * 2 - P;
}
bool onseg(Vec P, Line L) {
    if(disToLine(P, L) > EPS) return 0;
    return (P-L.A1) * (P-L.A2) <= 0;
}
Vec Lineintersect(Line L1, Line L2) {
    // 0 ~ (L2.A1 - L1.A1) ^ (L2.A2 - L1.A1);
    // 1 ~ (L2.A1 - L1.A2) ^ (L2.A2 - L1.A2);
    double F0 = (L2.A1 - L1.A1) ^ (L2.A2 - L1.A1);
    double F1 = (L2.A1 - L1.A2) ^ (L2.A2 - L1.A2);
    if (abs(F1 - F0) < EPS) return Vec(INF, INF);
    return L1.A1 + (L1.A2 - L1.A1) * (- F0 / (F1 - F0));
}
Vec Segintersect(Line L1, Line L2) {
    Vec P = Lineintersect(L1, L2);
    if (P.x == INF) return P;
    if (onseg(P, L1) && onseg(P, L2)) return P;
    return Vec(INF, INF);
}
bool disToSeg(Vec P, Line L) {
    double ans = min((P-L.A1).len(), (P-L.A2).len());
    ans = min(ans, (P - projection(P,L)).len());
    return ans;
}

// 
vector<Vec> shiftConvexHull(vector<Vec> A) {
    double my = INF, mx = INF; int st = -1;
    for (int i = 0; i < A.size(); i ++) {
        if (A[i].y < my || (A[i].y == my && A[i].x < mx) ) {
            my = A[i].y; mx = A[i].x; st = i;
        }
    }
    vector<Vec> ans;
    for (int i = 0; i < A.size(); i ++) {
        ans.push_back(A[(i+st)%A.size()]);
    }
    return ans;
}

vector<Vec> convexHull(vector<Vec> v) {
    sort(v.begin(), v.end());
    vector<Vec> ans; ans.resize(2*v.size());
    int top = 0, n = v.size();
    if(n <= 1) return v;

    for(int i = 0; i < n; ans[top++] = v[i++]) {
        while (top > 1 && ((ans[top-1] - ans[top-2]) ^ (v[i] - ans[top-1])) < EPS) -- top;
    }
    for(int i = n-2, t = top; i >= 0; ans[top++] = v[i--]) {
        while (top > t && ((ans[top-1] - ans[top-2]) ^ (v[i] - ans[top-1])) < EPS) -- top;
    } 
    while(ans.size() >= top) ans.pop_back();
    return ans;
}
double PolyArea(vector<Vec> v) {
    double ans = 0;
    int n = v.size();
    for(int i = 0; i < v.size(); i ++) {
        ans = ans + (v[i] ^ v[(i+1)%n]);
    }
    return abs(ans) / 2;
}
vector<Vec> convexCut(vector<Vec> v, Vec P1, Vec P2) {
    vector<Vec> ans;
    for(int i=0;i<v.size();i++) {
        int d1 = sgn((P2-P1)^(v[i]-P1));
        int d2 = sgn((P2-P1)^(v[(i+1)%v.size()]-P1));
        if (d1 >= 0) ans.push_back(v[i]);
        if (d1 * d2 < 0) ans.push_back(Lineintersect(Line(P1,P2), Line(v[i],v[(i+1)%v.size()])));
    }
    return ans;
}
double cloestPair(vector<Vec> v) {
    if (v.size() <= 1) return INF;

    sort(v.begin(), v.end());
    vector<Vec> vL, vR;
    int mid = v.size() / 2;
    for (int i = 0; i < v.size(); i ++) {
        if (i < mid)
            vL.push_back(v[i]);
        else
            vR.push_back(v[i]);
    }
    double d = min(cloestPair(vL), cloestPair(vR));
    
    vector<Vec> strip;
    for (int i = 0; i < v.size(); i ++) {
        if (abs(v[mid].x - v[i].x) <= d + EPS) 
            strip.push_back(v[i]);
    }

    for (int i = 0; i < strip.size(); i ++) swap(strip[i].x, strip[i].y);
    sort(strip.begin(), strip.end());
    for (int i = 0; i < strip.size(); i ++) swap(strip[i].x, strip[i].y);

    for (int i = 0; i < strip.size(); i ++) {
        for (int j = i + 1; j < strip.size(); j ++) {
            if (strip[j].y - strip[i].y > d + EPS) break;
            d = min(d, (strip[j]-strip[i]).len());
        }
    }
    return d;
}
bool onRight(Line L1, Line L2, Line L3) {
    Vec P = Lineintersect(L2, L3);
    return ((L1.A2 - L1.A1) ^ (P - L1.A1)) < 0;
}
vector<Vec> halfPlaneIntersection(vector<Line> L) {
    sort(L.begin(), L.end());

    int head = 0, tail = 0, cnt = 0; //[head, tail)
    Line que[L.size()+1];
    for (int i = 0; i < (int)L.size() - 1; i ++) {
        if (abs(L[i].getAng() - L[i+1].getAng()) < EPS) continue;
        L[cnt ++] = L[i];
    }

    L[cnt ++] = L[(int)L.size() - 1];

    for (int i = 0; i < cnt; i ++) {
        while (tail - head > 1 && onRight(L[i], que[tail-1], que[tail-2])) -- tail;
        while (tail - head > 1 && onRight(L[i], que[head], que[head+1])) ++ head;
        que[tail++] = L[i];
    }
    while (tail - head > 1 && onRight(que[tail-1], que[head], que[head+1])) ++ head;
    while (tail - head > 1 && onRight(que[head], que[tail-1], que[tail-2])) -- tail;
    
    vector<Vec> ans;
    for (int i = head; i < tail - 1; i ++) {
        ans.push_back(Lineintersect(que[i], que[i+1]));
    }
    ans.push_back(Lineintersect(que[tail - 1], que[head]));
    
    return ans;
}
vector<Vec> Minkowski(vector<Vec> a, vector<Vec> b) {
    vector<Line> L;
    for (int i = 0; i < a.size(); i ++) {
        L.push_back(Line(a[i], a[(i+1)%(int)a.size()]));
    }
    for (int i = 0; i < b.size(); i ++) {
        L.push_back(Line(b[i], b[(i+1)%(int)b.size()]));
    }
    sort(L.begin(), L.end());
    vector<Vec> ans;

    Vec now = Vec(0, 0); 
    ans.push_back(now);
    for (int i = 0; i < L.size(); i ++) {
        ans.push_back(now + L[i].A2 - L[i].A1);
        now = now + L[i].A2 - L[i].A1;
    }

    Vec LD = Vec(INF,INF);
    for (int i = 0; i < ans.size(); i ++) {
        LD = min(LD, ans[i]);
    }
    for (int i = 0; i < ans.size(); i ++) {
        ans[i] = ans[i] + (a[0] + b[0] - LD);        
    }
    return ans;
}
/*
Circle:
*/

struct Circle {
    Vec O; double r;
    void read() {
        O.read(); scanf("%lf", &r);
    }
};
vector<Vec> isCL(Circle C, Line L) {
    vector<Vec> ans;
    double dis = disToLine(C.O, L);
    if (dis > C.r + EPS) return ans;
    Vec P = projection(C.O, L);

    Vec v = (L.A1 - L.A2).change(sqrt(C.r * C.r - dis * dis));
    ans.push_back(P + v);
    ans.push_back(P - v);
    return ans;
}
vector<Vec> isCC(Circle C1, Circle C2) {
    vector<Vec> ans;
    double d = (C1.O - C2.O).len();

    if (C1.r + C2.r < d - EPS) return ans;
    if (abs(C1.r - C2.r) > d + EPS) return ans;

    double cos1 = (C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d);
    Vec v1 = (C2.O - C1.O).change(C1.r * cos1);
    Vec v2 = v1.rotate(acos(-1)/2).change(C1.r * sqrt(1 - cos1*cos1));
    ans.push_back(C1.O + v1 + v2);
    ans.push_back(C1.O + v1 - v2);
    return ans;
} 
vector<Vec> tanCP(Circle C, Vec P) {
    double x = (C.O - P).len();
    double d = x * x - C.r * C.r;
    vector<Vec> ans;
    if (d < -EPS) return ans;
    if (d < 0) d = 0;
    Vec q1 = (P - C.O) * ( (C.r * C.r) / (x * x) );
    Vec q2 = q1.rotate(acos(-1) / 2).change(C.r / x * (sqrt(x * x - C.r * C.r)));
    ans.push_back(C.O + q1 + q2);
    if (q2.len() > EPS) // unique
    ans.push_back(C.O + q1 - q2);
    return ans;
}
vector<Line> tanCC(Circle C1, Circle C2) {
    vector<Line> ans;
    if ( abs(C1.r - C2.r) < EPS ) {
        Vec dir = C2.O - C1.O;
        dir = (dir * (C1.r / dir.len())).rotate(acos(-1) / 2);
        ans.push_back(Line(C1.O + dir, C2.O + dir));
        ans.push_back(Line(C1.O - dir, C2.O - dir));
    } else { 
        Vec P = (C1.O * (- C2.r) + C2.O * C1.r) * (1.0/(C1.r - C2.r));
        vector<Vec> v1 = tanCP(C1, P);
        vector<Vec> v2 = tanCP(C2, P);
        for (int i = 0; i < v1.size() && i < v2.size(); i ++)
            ans.push_back(Line(v1[i], v2[i]));
    }
    Vec P = ((C1.O * C2.r) + (C2.O * C1.r)) * (1.0/(C1.r + C2.r));
    vector<Vec> v1 = tanCP(C1, P);
    vector<Vec> v2 = tanCP(C2, P);
    for (int i = 0; i < v1.size() && i < v2.size(); i ++)
        ans.push_back(Line(v1[i], v2[i]));
    return ans;
}
double AreaCC(Circle C1, Circle C2) {
    double d = (C1.O - C2.O).len();
    if (C1.r + C2.r > d + EPS) {
        return C1.r * C1.r * acos(-1) + C2.r * C2.r * acos(-1);
    }
    if (abs(C1.r - C2.r) < d - EPS) {
        return max(C1.r, C2.r) * max(C1.r, C2.r) * acos(-1);
    }
    double ans = C1.r * C1.r * acos(-1) + C2.r * C2.r * acos(-1);
    double cos1 = (C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d);
    double the1 = 2 * acos(cos1);

    double cos2 = (C2.r * C2.r + d * d - C1.r * C1.r) / (2 * C2.r * d);
    double the2 = 2 * acos(cos2);

    ans = ans - ( (2 * acos(-1) - the1) * C1.r * C1.r / 2 - sin(the1) * C1.r * C1.r / 2 );
    ans = ans - ( (2 * acos(-1) - the2) * C2.r * C2.r / 2 - sin(the2) * C2.r * C2.r / 2 );
    return ans;
}   

double AreaCT(Circle C1, Vec P1, Vec P2) {
    vector<Vec> is = isCL(C1, Line(P1,P2));
    
    if (is.size() == 0) return C1.r * C1.r * rad(P1, P2) / 2;
    bool b1 = P1.len() > C1.r + EPS, b2 = P2.len() > C1.r + EPS;

    if (b1 && b2) {
        if (onseg(is[0], Line(P1,P2)) && onseg(is[1], Line(P1,P2)))
            return (C1.r * C1.r * (rad(P1, P2) - rad(is[0], is[1])) + (is[0] ^ is[1])) / 2;
        return C1.r * C1.r * rad(P1, P2) / 2;
    } else if (b1) {
        return (C1.r * C1.r * rad(P1, is[0]) + (is[0] ^ P2)) / 2;
    } else if (b2) {
        return (C1.r * C1.r * rad(is[1], P2) + (P1 ^ is[1])) / 2;
    } 
    return (P1 ^ P2) / 2;
}

vector<Vec> v;
int n; Circle C;
int main(){
    scanf("%d %lf", &n, &C.r);
    C.O = Vec(0, 0);
    for (int i = 0; i < n; i ++) {
        Vec p; p.read(); v.push_back(p);
    }
    double ans = 0;
    for (int i = 0; i < n; i ++) {
        ans = ans + AreaCT(C, v[i], v[(i+1)%n]);
    }
    printf("%.8f\n", ans);
}