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<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define rep(i, a, b) for(int i = (a); i < (b); ++i)
#define per(i, a, b) for(int i = (b) - 1; i >= (a); --i)
#define sz(a) (int)a.size()
#define de(c) cout << #c << " = " << c << endl
#define dd(c) cout << #c << " = " << c << " "
#define all(a) a.begin(), a.end()
#define pw(x) (1ll<<(x))
#define endl "\n"
typedef long long ll;
typedef double db;
typedef pair<int, int> pii;
typedef vector<int> vi;

typedef db T;
const db eps = 1e-9 , pi = acosl(-1.);
int sgn(T x){return (x>eps)-(x<-eps);}
struct P{
  T x,y; P(){} P(T x,T y):x(x),y(y){}
  P operator - (const P&b) const {return P(x-b.x,y-b.y);}
  P operator + (const P&b) const {return P(x+b.x,y+b.y);}
  T operator * (const P&b) const {return x*b.x+y*b.y;}
  T operator / (const P&b) const {return x*b.y-y*b.x;}
  P operator * (const T&k) const {return P(x*k,y*k);}
  P operator / (const T&k) const {return P(x/k,y/k);}
	bool operator < (const P&b) const {return sgn(x-b.x)?x<b.x:y<b.y;}
	bool operator == (const P&b) const{return !sgn(x-b.x)&&!sgn(y-b.y);}
	bool operator != (const P&b) const{return !(*this == b);}
	P rot90(){return P(-y,x);}
};
T norm(P a){return a*a;}
T abs(P a) {return sqrtl(norm(a));}
P proj(P p,P a,P b){return (b-a)*((p-a)*(b-a)/norm(b-a))+a;}
P reflect(P p,P a,P b){return proj(p,a,b)*2-p;}
T cross(P o,P a,P b){return (a-o)/(b-o);}
int crossOp(P o,P a,P b){return sgn(cross(o,a,b));}
db rad(P p1,P p2){return atan2l(p1/p2,p1*p2);}
bool onPS(P p,P s,P t){return sgn((t-s)/(p-s))==0&&sgn((p-s)*(p-t))<=0;}
struct L{ P s,t;L(){} L(P s,P t):s(s),t(t){}};
P insLL(L a,L b){ // line x line
	P s = a.s - b.s , v = a.t - a.s , w = b.t - b.s;
	db k1 = s / w , k2 = w / v;
	if(sgn(k2) == 0) return abs(b.s - a.s) < abs(b.t - a.s) ? b.s : b.t;
	return a.s + v * (k1 / k2);
}
bool isSS(L a,L b){ // seg x seg , replace x->y to accelerate
	T c1=(a.t-a.s)/(b.s-a.s),c2=(a.t-a.s)/(b.t-a.s);
	T c3=(b.t-b.s)/(a.s-b.s),c4=(b.t-b.s)/(a.t-b.s);
	return sgn(c1) * sgn(c2) <= 0 && sgn(c3) * sgn(c4) <= 0 &&
		sgn(max(a.s.x,a.t.x) - min(b.s.x,b.t.x)) >= 0 &&
		sgn(max(b.s.x,b.t.x) - min(a.s.x,a.t.x)) >= 0 &&
		sgn(max(a.s.y,a.t.y) - min(b.s.y,b.t.y)) >= 0 &&
		sgn(max(b.s.y,b.t.y) - min(a.s.y,a.t.y)) >= 0;
}
db disPL(P p,L a){return fabs((a.t-a.s)/(p-a.s)) / abs(a.t-a.s);}
db disPS(P p,L a){ // p x seg dis
	if(sgn((a.t-a.s)*(p-a.s)) == -1) return abs(p-a.s);
	if(sgn((a.s-a.t)*(p-a.t)) == -1) return abs(p-a.t);
	return disPL(p,a);
}
db disSS(L a,L b){ // seg x seg dis
	if(isSS(a,b)) return 0;
	return min(min(disPS(a.s,b),disPS(a.t,b)),min(disPS(b.s,a),disPS(b.t,a)));
}
typedef vector<P> polygon;
polygon convex(polygon A){ // counter-clockwise , < : <=180 , <= : <180
  int n=sz(A),m=0;
  polygon B;B.resize(n<<1);
  sort(all(A));
  rep(i,0,n){
    while(m > 1 && sgn((B[m-1]-B[m-2])/(A[i]-B[m-2]))<0) --m;
    B[m++]=A[i];
  }
  int k = m;
  per(i,0,n-1){
    while(m > k && sgn((B[m-1]-B[m-2])/(A[i]-B[m-2]))<0) --m;
    B[m++]=A[i];
  }
  B.resize(m);
  if(sz(B) > 1) B.pop_back();
  return B;
}
T area(polygon A) { // multiple 2 with integer type
	T res=0;
	rep(i,0,sz(A)) res+=A[i]/(A[(i+1)%sz(A)]);
	return fabs(res) / 2;
}
bool isconvex(polygon A){ // counter-clockwise
	bool ok=1;int n=sz(A);
	rep(i,0,2) A.pb(A[i]);
	rep(i,0,n) ok&=((A[i+1]-A[i])/(A[i+2]-A[i]))>=0;
	return ok;
}
int inPpolygon(P p,polygon A){ // -1 : on , 0 : out , 1 : in
	int res=0;
	rep(i,0,sz(A)){
		P u=A[i],v=A[(i+1)%sz(A)];
		if(onPS(p,u,v)) return -1;
		T cross = sgn((v-u)/(p-u)) , d1 = sgn(u.y-p.y) , d2 = sgn(v.y-p.y);
		if(cross > 0 && d1 <= 0 && d2 > 0) ++res;
		if(cross < 0 && d2 <= 0 && d1 > 0) --res;
	}
	return res != 0;
}
T diameter(polygon A) { // longest distance
  int n=sz(A);if(n <= 1) return 0;
  int l=0,r=0;rep(i,1,n) (A[i]<A[l])&&(l=i),(A[r]<A[i])&&(r=i);
  db res=abs(A[l]-A[r]);int i=l,j=r;
  do (++((A[(i+1)%n]-A[i])/(A[(j+1)%n]-A[j])>=0?j:i))%=n,
    res=max(res,abs(A[i]-A[j]));
  while(i!=l||j!=r);
  return res;
}
polygon convexCut(polygon A,P s,P t){ // counter-clockwise , left hand of st
	int n=sz(A);
	polygon B;
	rep(i,0,n){
		P u=A[i],v=A[(i+1)%n];
		int d1 = sgn((t-s)/(u-s)) , d2 = sgn((t-s)/(v-s));
		if(d1 >= 0) B.pb(u);
		if(d1 * d2 < 0) B.pb(insLL(L(u,v),L(s,t)));
	}
	return B;
}
namespace NearestPoints{
  T solve(int l,int r,vector<P>&p){
    if(l == r) return 1e100;
    int m=(l+r)>>1;
    T Xm = p[m].x , lim = min(solve(l,m,p) , solve(m+1,r,p));
    inplace_merge(p.begin()+l,p.begin()+m+1,p.begin()+r+1,[&](P a,P b){return a.y<b.y;});
    vector<P> V;
    rep(i,l,r+1) if(fabs(p[i].x - Xm) <= lim) V.pb(p[i]);
    rep(i,0,sz(V)) rep(j,i+1,sz(V)){
      if(fabs(V[j].y - V[i].y) >= lim) break;
      T dis = abs(V[i]-V[j]);
      lim = min(lim,dis);
    }
    return lim;
  }
  T solve(vector<P> A){
    sort(all(A),[&](P a,P b){return a.x<b.x;});
    return solve(0,sz(A)-1,A);
  }
}
struct C{
	P o;T r;C(){} C(P o,T r):o(o),r(r){}
	bool operator == (const C&b) const {return o==b.o&&sgn(r-b.r)==0;}
};
int relCC(C A,C B){
	T dis = abs(A.o - B.o);
	if(sgn(dis - (A.r + B.r)) == 1) return 4;
	if(sgn(dis - (A.r + B.r)) == 0) return 3;
	if(sgn(dis - fabs(A.r - B.r)) == 1) return 2;
	if(sgn(dis - fabs(A.r - B.r)) == 0) return 1;
	return 0;
}
vector<P> insCL(C c,L a){
	db x = (a.s-c.o)*(a.t-a.s) , y = norm(a.t-a.s);
	db d = x * x - y * (norm(a.s-c.o) - c.r*c.r);
	vector<P> res;
	if(sgn(d) < 0) return res;
	d = max(d,0.);
	P mid = a.s - (a.t - a.s) * (x / y);
	P del = (a.t - a.s) * (sqrt(d) / y);
	return {mid - del,mid + del}; // dir : a.s -> a.t
}
vector<P> insCC(C a,C b){
	vector<P> res;
	T x = norm(a.o - b.o);
	if(sgn(x)==0) return res;
	T  y = ((a.r * a.r - b.r * b.r) / x + 1) / 2 ,
		 d = a.r * a.r / x - y * y;
	if(sgn(d) < 0) return res;
	d = max(d,0.);
	P mid = (b.o - a.o) * y + a.o ,
		del = ((b.o - a.o) * sqrt(d)).rot90();
	return {mid - del , mid + del};// counter-clockwise along a
}
vector<P> tanCP(C c,P p){
	db x = norm(p - c.o) , d = x - c.r * c.r;
	vector<P> res;
	if(sgn(d) < 0) return res;
	d = max(d,0.);
	P mid = c.o + (p - c.o) * (c.r * c.r / x) ,
		del = ((p-c.o)*(c.r*sqrt(d)/x)).rot90();
	return {mid - del ,mid + del}; // counter-clockwise
}
vector<pair<P,P> > tanCC(C c1,C c2){// need to unique
	vector<pair<P,P> > res;
	// extan
	if(!sgn(c1.r-c2.r)){
		P dir = c2.o-c1.o;
		dir = (dir*(c1.r/abs(dir))).rot90();
		res.pb({c1.o+dir,c2.o+dir});
		res.pb({c1.o-dir,c2.o-dir});
	} else {
		P p = (c2.o * c1.r - c1.o * c2.r) / (c1.r - c2.r);
		vector<P> ps = tanCP(c1 , p) , qs = tanCP(c2 , p);
		rep(i,0,min(sz(ps),sz(qs))) res.pb({ps[i],qs[i]});
	}
	// intan
	P p = (c1.o * c2.r + c2.o * c1.r) / (c1.r + c2.r);
	vector<P> ps = tanCP(c1 , p) , qs = tanCP(c2 , p);
	rep(i,0,min(sz(ps),sz(qs))) res.pb({ps[i],qs[i]});
	return res;
}
db areaCT(db r,P s,P t) { // need divide 2
	vector<P> p = insCL(C(P(0,0),r),L(s,t));
	if(!sz(p)) return r*r*rad(s,t);
	bool b1 = sgn(norm(s)-r*r) == 1 , b2 = sgn(norm(t)-r*r) == 1;
	if(b1 && b2) {
		if(sgn((s-p[0])*(t-p[0])) <= 0 && sgn((s-p[1])*(t-p[1])) <= 0)
			return r*r*(rad(s,p[0]) + rad(p[1],t)) + (p[0]/p[1]);
		else return r*r*rad(s,t);
	} else if(b1) return r*r*rad(s,p[0])+(p[0]/t);
	else if(b2) return r*r*rad(p[1],t)+(s/p[1]);
	return (s/t);
}
db areaCPoly(db r, polygon A) { // need divide 2, counter-clockwise
	db ans = 0;
	rep(i, 0, sz(A)) ans += areaCT(r, A[i], A[(i + 1) % sz(A)]);
	return ans;
}

int n, r;
P p;
polygon A;

int main() {
	std::ios::sync_with_stdio(0);
	std::cin.tie(0);
	//	freopen("a.in", "r", stdin);
	//	freopen("a.out", "w", stdout);
	cout << setiosflags(ios::fixed);
	cout << setprecision(10);
	cin >> n >> r;
	rep(i, 0, n) {
		cin >> p.x >> p.y;
		A.pb(p);
	}
	cout << areaCPoly(r, A) / 2 << endl;
	return 0;
}