Problem: Example

Input

4 1
0 0
3 0
3 3
0 3
1 1


Output

2.0000000000

Solution (in CPP):

Solution:
#include <bits/stdc++.h>
using namespace std;
#define FOR(i,k,n) for(int i = (int)(k); i < (int)(n); i++)
#define REP(i,n) FOR(i,0,n)
#define ALL(a) a.begin(), a.end()
#define MS(m,v) memset(m,v,sizeof(m))
typedef long long ll;
typedef long double ld;
typedef vector<int> vi;
typedef vector<string> vs;
typedef pair<int, int> pii;
const int MOD = 1000000007;
const int INF = MOD + 1;
const ld EPS = 1e-12;
template<class T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template<class T> T &chmax(T &a, const T &b) { return a = max(a, b); }

/*--------------------template--------------------*/

ld calc(ld h, ld a, ld b)
{
	ld res = 0.25 * pow(h, 4) * ((1.0 / 3.0) * (pow(a, 3) - pow(b, 3)) + (a - b));
	cout << res << endl;
	return res;
}

typedef long double ld;
const ld PI = acos(-1.0);
bool eq(ld a, ld b) { return abs(a - b) < EPS; }
typedef complex<ld> Point;
typedef vector<Point> Polygon;

namespace std
{
	bool operator < (const Point& a, const Point& b)
	{
		return real(a) != real(b) ? real(a) < real(b) : imag(a) < imag(b);
	}
}

struct Line
{
	Point a, b;
	Line() {};
	Line(Point p, Point q) :a(p), b(q) {};
	Line(ld x1, ld y1, ld x2, ld y2) :a(Point(x1, y1)), b(Point(x2, y2)) {};
};

struct Circle
{
	Point p; ld r;
	Circle() {};
	Circle(Point a, ld b) :p(a), r(b) {};
};

ld dot(Point a, Point b)
{
	return real(conj(a)*b);
}

ld cross(Point a, Point b)
{
	return imag(conj(a)*b);
}

int ccw(Point a, Point b, Point c)
{
	b -= a; c -= a;
	if (cross(b, c) > EPS) return 1;    //counter cloclwise
	if (cross(b, c) < -EPS) return -1;  //cloclwise
	if (dot(b, c) < 0) return 2;        //c--a--b on line 
	if (norm(b) < norm(c)) return -2;   //a--b--c on line
	return 0;                           //a--c--b on line
}

bool isis_ll(Line l, Line m)
{
	return abs(cross(l.b - l.a, m.b - m.a)) > EPS;
}

bool isis_ls(Line l, Line s)
{
	return cross(l.b - l.a, s.a - l.a)*cross(l.b - l.a, s.b - l.a) < EPS;
}

bool isis_ss(Line s, Line t)
{
	return (ccw(s.a, s.b, t.a)*ccw(s.a, s.b, t.b) <= 0 &&
		ccw(t.a, t.b, s.a)*ccw(t.a, t.b, s.b) <= 0);
}

bool isis_lp(Line l, Point p)
{
	return (abs(cross(l.b - p, l.a - p)) < EPS);
}

bool isis_sp(Line s, Point p)
{
	return (abs(s.a - p) + abs(s.b - p) - abs(s.b - s.a)) < EPS;
}

Point projection(Line l, Point p)
{
	Point base = l.b - l.a;
	ld r = dot(p - l.a, base) / norm(base);
	return l.a + base*r;
}

Point mirror(Line l, Point p)
{
	return Point(2, 0)*projection(l, p) - p;
}

ld dist_lp(Line l, Point p)
{
	return abs(p - projection(l, p));
}

ld dist_ll(Line l, Line m)
{
	return isis_ll(l, m) ? 0 : dist_lp(l, m.a);
}

ld dist_ls(Line l, Line s)
{
	if (isis_ls(l, s)) return 0;
	return min(dist_lp(l, s.a), dist_lp(l, s.b));
}

ld dist_sp(Line s, Point p)
{
	Point r = projection(s, p);
	if (isis_sp(s, r)) return abs(r - p);
	return min(abs(s.a - p), abs(s.b - p));
}

ld dist_ss(Line s, Line t)
{
	if (isis_ss(s, t)) return 0;
	return min(min(dist_sp(s, t.a), dist_sp(s, t.b)), min(dist_sp(t, s.a), dist_sp(t, s.b)));
}

Point is_ll(Line s, Line t)
{
	ld a = cross(s.b - s.a, t.b - t.a);
	ld b = cross(s.b - s.a, s.b - t.a);
	return t.a + b / a*(t.b - t.a);
}

ld area(const Polygon &p) {
	ld res = 0;
	int n = p.size();
	REP(j, n) res += cross(p[j], p[(j + 1) % n]);
	return res / 2;
}


Polygon ConvexCut(const Polygon& g, const Line l)
{
	Polygon res;
	REP(i, g.size())
	{
		Point a = g[i], b = g[(i + 1) % g.size()];
		if (ccw(l.a, l.b, a) != -1) res.push_back(a);
		if (ccw(l.a, l.b, a)*ccw(l.a, l.b, b) < 0) res.push_back(is_ll(Line(a, b), l));
	}
	return res;
}

Line bisecter(Point p, Point q)
{
	Point m = (p + q) / Point(2, 0);
	Point v = p - q;
	v = Point(v.imag(), -v.real());
	return Line(m, m + v);
}


Polygon voronoi(const Polygon& g, const vector<Point>& ps, const int p)
{
	Polygon res = g;
	REP(i, ps.size())
	{
		if (i == p) continue;
		Line l = bisecter(ps[p], ps[i]);
		res = ConvexCut(res, l);
	}
	return res;
}

ld calc(Point a, Point b, Point c)
{
	if (abs(a - b) < EPS || abs(b - c) < EPS || abs(c - a) < EPS) return 0;
	ld p = abs(a - b), q = abs(a - c), r = abs(b - c);
	ld s = (p + q + r) / 2;
	ld area = sqrt(s * (s - p) * (s - q) * (s - r));
	ld h = area * 2 / r;
	if (p < q) swap(p, q);
	ld pp = (p - h > EPS ? sqrt(p*p - h*h) : 0);
	ld rp = pp / h;
	ld qq = (q - h > EPS ? sqrt(q*q - h*h) : 0);
	ld rq = qq / h;
	if (r*r > p * p - h * h) rq *= -1;
	ld res = pow(h, 4) * ((1.0 / 3.0) * (pow(rp, 3) - pow(rq, 3)) + rp - rq) / 4;
	return res;
}

ld solve(const Polygon& g, const Point p)
{
	ld res = 0;
	REP(i, g.size())
	{
		Point a = g[i], b = g[(i + 1) % g.size()];
		res += calc(p, a, b);
	}
	return res;
}

int main()
{
	cin.sync_with_stdio(false); cout << fixed << setprecision(10);
	int m, n;
	cin >> m >> n;
	Polygon g;
	REP(i, m)
	{
		ld x, y; cin >> x >> y;
		g.emplace_back(x, y);
	}
	vector<Point> ps;
	REP(i, n)
	{
		ld x, y; cin >> x >> y;
		ps.emplace_back(x, y);
	}
	ld ans = 0;
	REP(i, n)
	{
		Polygon tmp = voronoi(g, ps, i);
		ans += solve(tmp, ps[i]);
	}
	cout << ans / area(g) << endl;
	return 0;
}