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>
#include<vector>
#include<algorithm>
#include<iostream>
#include<deque>
#include<cmath>
#include<cstdlib>
using namespace std;
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#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(x) cout << #x << " = " << x << endl
#define dd(x) cout << #x << " = " << x << " "
#define endl "\n"
#define all(a) a.begin(), a.end()
#define pw(a) (1ll << (a))
typedef long long ll;
typedef pair<int, int> pii;
typedef double db;
typedef vector<int> vi;

const db eps = 1e-8;

int sign(db x) {
	return (x > eps) - (x < -eps);
}

struct P {
	db x, y;
	P() {}
	P(db x, db y) : x(x), y(y) {}
	void read() {
		cin >> x >> y;
	}
	P operator + (const P &c) const {
		return P(x + c.x, y + c.y);
	}
	P operator - (const P &c) const {
		return P(x - c.x, y - c.y);
	}
	P operator * (const db &c) const {
		return P(x * c, y * c);
	}
	P operator / (const db &c) const {
		return P(x / c, y / c);
	}
	db len2() {
		return x * x + y * y;
	}
	db len() {
		return sqrt(x * x + y * y);
	}
	P turn90() {
		return P(-y, x);
	}
	bool operator == (const P &c) const {
		return sign(x - c.x) == 0 && sign(y - c.y) == 0;
	}
	/*
	bool operator < (const P &c) const {
		int o = sign(x - c.x);
		return o ? o < 0 : sign(y - c.y) < 0;
	}
	*/
	int quad() const {
		return sign(y) == 1 || (sign(y) == 0 && sign(x) >= 0);
	}
	P unit() {
		return *this / len();
	}
};
db det(P a, P b) {
	return a.x * b.y - a.y * b.x;
}
struct L {
	P a, b;
	L() {}
	L(P a, P b) : a(a), b(b) {}
	void read() {
		a.read();
		b.read();
	}
	bool include(const P &p) const {
		return sign(det(b - a, p - a)) > 0;
	}
	L push() const {
		const db eps = 1e-6;
		P delta = (b - a).turn90().unit() * eps;
		return L(a - delta, b - delta);
	}
};
struct C {
	P o;
	db r;
	C() {}
	C(P o, db r) : o(o), r(r) {}
	void read() { 
		o.read();
		cin >> r;
	}
};
db det(P o, P a, P b) {
	return det(a - o, b - o);
}
db dot(P a, P b) {
	return a.x * b.x + a.y * b.y;
}
bool onSeg(const L &l, const P &p) {
	return sign(det(p - l.a, l.b - l.a)) == 0 && sign(dot(p - l.a, p - l.b)) <= 0;
}
P isLL(const L &l1, const L &l2) {
	db s1 =  det(l2.b - l2.a, l1.a - l2.a);
	db s2 = -det(l2.b - l2.a, l1.b - l2.a);
	return (l1.a * s2 + l1.b * s1) / (s1 + s2);
}
P proj(const L &l, const P &p) {
	return l.a + (l.b - l.a) * (dot(p - l.a, l.b - l.a) / (l.b - l.a).len2());
}
db disToLine(const L &l, const P &p) {
	return fabs(det(p - l.a, l.b - l.a) / (l.b - l.a).len());
}
db disToSeg(const L &l, const P &p) {
	return sign(dot(p - l.a, l.b - l.a)) * sign(dot(p - l.b, l.a - l.b)) == 1 ? disToLine(l, p) : min((p - l.a).len(), (p - l.b).len());
}
bool tanCP(const C &c, const P &p0, P &p1, P &p2) {
	db x = (p0 - c.o).len2(), d = x - c.r * c.r;
	if(d < eps) return 0;
	P p = (p0 - c.o) * (c.r * c.r / x);
	P delta = ((p0 - c.o) * (-c.r * sqrt(d) / x)).turn90();
	p1 = c.o + p + delta;
	p2 = c.o + p - delta;
	return 1;
}
vector<L> extanCC(const C &c1, const C &c2) {
	vector<L> ret;
	if(sign(c1.r - c2.r) == 0) {
		P dir = c2.o - c1.o;
		dir = (dir * (c1.r / dir.len())).turn90();
		ret.pb(L(c1.o + dir, c2.o + dir));
		ret.pb(L(c1.o - dir, c2.o - dir));
	} else {
		P p = (c1.o * -c2.r + c2.o * c1.r) / (c1.r - c2.r);
		P p1, p2, q1, q2;
		if(tanCP(c1, p, p1, p2) && tanCP(c2, p, q1, q2)) {
			if(c1.r < c2.r) swap(p1, p2), swap(q1, q2);
			ret.pb(L(p1, q1));
			ret.pb(L(p2, q2));
		}
	}
	return ret;
}
vector<L> intanCC(const C &c1, const C &c2) {
	vector<L> ret;
	P p = (c1.o * c2.r + c2.o * c1.r) / (c1.r + c2.r);
	P p1, p2, q1, q2;
	if(tanCP(c1, p, p1, p2) && tanCP(c2, p, q1, q2)) {
		ret.pb(L(p1, q1));
		ret.pb(L(p2, q2));
	}
	return ret;
}
vector<P> tanCC(const C &c1, const C &c2) { // 求圆与圆的切点
	vector<P> res;
	db dis = (c1.o - c2.o).len();
	if(sign(dis - (c1.r + c2.r)) == 0) {
		res.pb(c1.o + (c2.o - c1.o) * c1.r / (c1.r + c2.r));
	}
	if(sign(dis - fabs(c1.r - c2.r) == 0)) {
		res.pb(c1.o + (c2.o - c1.o) * c1.r / (c1.r - c2.r));
	}
	return res;
}
int contain(vector<P> poly, P p) {
	int ret = 0, n = sz(poly);
	rep(i, 0, n) {
		P u = poly[i], v = poly[(i + 1) % n];
		if(onSeg(L(u, v), p)) return 1;
		if(sign(u.y - v.y) <= 0) swap(u, v);
		if(sign(p.y - u.y) > 0 || sign(p.y - v.y) <= 0) continue;
		ret += sign(det(p, v, u)) > 0;
	}
	return (ret & 1) * 2;
}
vector<P> convexCut(const vector<P> &ps, L l) {
	vector<P> qs;
	int n = sz(ps);
	rep(i, 0, n) {
		P p1 = ps[i], p2 = ps[(i + 1) % n];
		int d1 = sign(det(l.a, l.b, p1)), d2 = sign(det(l.a, l.b, p2));
		if(d1 >= 0) qs.pb(p1);
		if(d1 * d2 < 0) qs.pb(isLL(L(p1, p2), l));
	}
	return qs;
}
db area(const vector<P> &poly) {
	int n = sz(poly);
	db res = 0;
	rep(i, 0, n) res += det(poly[i], poly[(i + 1) % n]);
	return fabs(res) / 2;
}
vector<P> convexHull(vector<P> ps) {
	int n = sz(ps);
	if(n <= 1) return ps;
	sort(all(ps));
	vector<P> qs;
	for(int i = 0; i < n; qs.pb(ps[i++])) {
		while(sz(qs) > 1 && sign(det(qs[sz(qs) - 2], qs.back(), ps[i])) < 0) qs.pop_back();
	}
	for(int i = n - 2, t = sz(qs); i >= 0; qs.pb(ps[i--])) {
		while(sz(qs) > t && sign(det(qs[sz(qs) - 2], qs.back(), ps[i])) < 0) qs.pop_back();
	}
	qs.pop_back();
	return qs;
}
bool sameDir(const L &l0, const L &l1) {
	return sign(det(l0.a - l0.b, l1.a - l1.b)) == 0 && sign(dot(l0.b - l0.a, l1.b - l1.a)) == 1;
}
bool operator < (const P &a, const P &b) {
	if(a.quad() != b.quad()) {
		return a.quad() < b.quad();
	} else {
		return sign(det(a, b)) > 0;
	}
}
bool operator < (const L &l0, const L &l1) {
	if(sameDir(l0, l1)) {
		return l1.include(l0.a);
	} else {
		return (l0.b - l0.a) < (l1.b - l1.a);
	}
}
bool check(const L &u, const L &v, const L &w) {
	return w.include(isLL(u, v));
}
vector<P> intersection(vector<L> &l) {
	sort(all(l));
	deque<L> q;
	rep(i, 0, sz(l)) {
		if(i && sameDir(l[i], l[i - 1])) continue;
		while(sz(q) > 1 && !check(q[sz(q) - 2], q.back(), l[i])) q.pop_back();
		while(sz(q) > 1 && !check(q[1], q[0], l[i])) q.pop_front();
		q.pb(l[i]);
	}
	while(sz(q) > 2 && !check(q[sz(q) - 2], q.back(), q[0])) q.pop_back();
	while(sz(q) > 2 && !check(q[1], q[0], q.back())) q.pop_front();
	vector<P> ret;
	rep(i, 0, sz(q)) ret.pb(isLL(q[i], q[(i + 1) % sz(q)]));
	return ret;
}
bool isCL(C a, L l, P &p1, P &p2) {
	db x = dot(l.a - a.o, l.b - l.a), y = (l.b - l.a).len2(), d = x * x - y * ((l.a - a.o).len2() - a.r * a.r);
	if(sign(d) < 0) return 0;
	d = max(d, 0.0);
	P p = l.a - ((l.b - l.a) * (x / y)), delta = (l.b - l.a) * (sqrt(d) / y);
	p1 = p - delta, p2 = p + delta;
	return 1;
}
db rad(P p1, P p2) {
	return atan2l(det(p1, p2), dot(p1, p2));
}
db areaCT(db r,P s,P t) { // need divide 2
	P p1, p2;
	bool f = isCL(C(P(0, 0), r), L(s, t), p1, p2);
	if(!f) return r * r * rad(s, t);
	bool b1 = sign(s.len2() - r * r) == 1 , b2 = sign(t.len2() - r * r) == 1;
	if(b1 && b2) {
		if(sign(dot(s - p1, t - p1)) <= 0 && sign(dot(s - p2, t - p2) <= 0))
			return r * r * (rad(s, p1) + rad(p2, t)) + det(p1, p2);
		else return r * r * rad(s, t);
	} else if(b1) return r * r * rad(s, p1) + det(p1, t);
	else if(b2) return r * r * rad(p2, t) + det(s, p2);
	return det(s, t);
}
db areaCPoly(C c, vector<P> p) {
	int n = sz(p);
	db ans = 0;
	rep(i, 0, n) {
		P u = p[i], v = p[(i + 1) % n];
		ans += areaCT(c.r, u - c.o, v - c.o);
	}
	return fabs(ans) / 2;
}

int main() {
	std::ios::sync_with_stdio(0);
	std::cin.tie(0);
	int n, r;
	cin >> n >> r;
	vector<P> p(n);
	rep(i, 0, n) cin >> p[i].x >> p[i].y;
	cout << setiosflags(ios::fixed);
	cout << setprecision(12);
	cout << areaCPoly(C(P(0, 0), r), p) << endl;
	return 0;
}