Problem: Example

Input

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


Output

2.0000000000

Solution (in CPP):

Solution:
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
using namespace std;
const double EPS = 1e-10;
const double INF = 1e+10;
const double PI = acos(-1);
inline double ABS(double a){return max(a,-a);}
int sig(double r) { return (r < -EPS) ? -1 : (r > +EPS) ? +1 : 0; }
struct Pt {
	double x, y;
	Pt() {}
	Pt(double x, double y) : x(x), y(y) {}
	Pt operator+(const Pt &a) const { return Pt(x + a.x, y + a.y); }
	Pt operator-(const Pt &a) const { return Pt(x - a.x, y - a.y); }
	Pt operator*(const Pt &a) const { return Pt(x * a.x - y * a.y, x * a.y + y * a.x); }
	Pt operator-() const { return Pt(-x, -y); }
	Pt operator*(const double &k) const { return Pt(x * k, y * k); }
	Pt operator/(const double &k) const { return Pt(x / k, y / k); }
	double abs() const { return sqrt(x * x + y * y); }
	double abs2() const { return x * x + y * y; }
	double arg() const { return atan2(y, x); }
	double dot(const Pt &a) const { return x * a.x + y * a.y; }
	double det(const Pt &a) const { return x * a.y - y * a.x; }
};
double tri(const Pt &a, const Pt &b, const Pt &c) { return (b - a).det(c - a); }
int iSP(Pt a, Pt b, Pt c) {
	int s = sig((b - a).det(c - a));
	if (s) return s;
	if (sig((b - a).dot(c - a)) < 0) return -2; // c-a-b
	if (sig((a - b).dot(c - b)) < 0) return +2; // a-b-c
	return 0;
}
int iLL(Pt a, Pt b, Pt c, Pt d) {
	if (sig((b - a).det(d - c))) return 1; // intersect
	if (sig((b - a).det(c - a))) return 0; // parallel
	return -1; // correspond
}
bool iLS(Pt a, Pt b, Pt c, Pt d) {
	return (sig(tri(a, b, c)) * sig(tri(a, b, d)) <= 0);
}
bool iSS(Pt a, Pt b, Pt c, Pt d) {
	return (iSP(a, b, c) * iSP(a, b, d) <= 0 && iSP(c, d, a) * iSP(c, d, b) <= 0);
}
bool iSSstrict(Pt a, Pt b, Pt c, Pt d) {
	return (sig(tri(a, b, c)) * sig(tri(a, b, d)) < 0 && sig(tri(c, d, a)) * sig(tri(c, d, b)) < 0);
}
Pt pLL(Pt a, Pt b, Pt c, Pt d) {
	b = b - a; d = d - c; return a + b * (c - a).det(d) / b.det(d);
}
double dLP(Pt a, Pt b, Pt c) {
	return ABS(tri(a, b, c)) / (b - a).abs();
}
double dSP(Pt a, Pt b, Pt c) {
	if (sig((b - a).dot(c - a)) <= 0) return (c - a).abs();
	if (sig((a - b).dot(c - b)) <= 0) return (c - b).abs();
	return ABS(tri(a, b, c)) / (b - a).abs();
}
double dLL(Pt a, Pt b, Pt c, Pt d) {
	return iLL(a, b, c, d) ? 0 : dLP(a, b, c);
}
double dLS(Pt a, Pt b, Pt c, Pt d) {
	return iLS(a, b, c, d) ? 0 : min(dLP(a, b, c), dLP(a, b, d));
}
double dSS(Pt a, Pt b, Pt c, Pt d) {
	return iSS(a, b, c, d) ? 0 : min(min(dSP(a, b, c), dSP(a, b, d)), min(dSP(c, d, a), dSP(c, d, b)));
}
vector<Pt> convexCut(int n, vector<Pt>p, Pt a, Pt b) {
	int m = 0, i;
	p.push_back(p[0]);
	vector<Pt>q;
	for (i = 0; i < n; ++i) {
	//	printf("%d %d\n",i,q.size());
		if (sig(tri(a, b, p[i])) >= 0) q.push_back(p[i]);
		if (sig(tri(a, b, p[i])) * sig(tri(a, b, p[i + 1])) < 0)
			q.push_back(pLL(a, b, p[i], p[i + 1]));
	}
	return q;
}
Pt c[123];
Pt d[123];
vector<vector<Pt> > poly[112];
double calc(vector<Pt>v){
	double ret=0;
	vector<double> x;
	vector<double> y;
	for(int i=0;i<v.size();i++){
		x.push_back(v[i].x-EPS);
		x.push_back(v[i].x+EPS);
		y.push_back(v[i].y-EPS);
		y.push_back(v[i].y+EPS);
	}
	std::sort(x.begin(),x.end());
	std::sort(y.begin(),y.end());
	double L,R;
	int n=v.size();
	for(int i=0;i<v.size()*2-1;i++){
		if(x[i]+1e-9>x[i+1])continue;
		L=99999999;
		R=-99999999;
		for(int j=0;j<n;j++){
			if(min(v[j].x,v[(j+1)%n].x)<x[i]&&x[i]<max(v[j].x,v[(j+1)%n].x)){
				double py=v[j].y+(v[(j+1)%n].y-v[j].y)*(x[i]-v[j].x)/(v[(j+1)%n].x-v[j].x);
				L=min(L,py);R=max(R,py);
			}
		}
		double t1=R-L;
		L=99999999;
		R=-99999999;
		for(int j=0;j<n;j++){
			if(min(v[j].x,v[(j+1)%n].x)<x[i+1]&&x[i+1]<max(v[j].x,v[(j+1)%n].x)){
				double py=v[j].y+(v[(j+1)%n].y-v[j].y)*(x[i+1]-v[j].x)/(v[(j+1)%n].x-v[j].x);
				L=min(L,py);R=max(R,py);
			}
		}
		double t2=R-L;
		double val=(t1*x[i]*x[i]+t2*x[i+1]*x[i+1]+(t1+t2)*(x[i]+x[i+1])*(x[i]+x[i+1])/2)*(x[i+1]-x[i])/6;
		//printf("%f %f: %f %f %f\n",x[i],x[i+1],t1,t2,val);
		ret+=val;
	}
	for(int i=0;i<v.size()*2-1;i++){
		if(y[i]+1e-9>y[i+1])continue;
		L=99999999;
		R=-99999999;
		for(int j=0;j<n;j++){
			if(min(v[j].y,v[(j+1)%n].y)<y[i]&&y[i]<max(v[j].y,v[(j+1)%n].y)){
				double px=v[j].x+(v[(j+1)%n].x-v[j].x)*(y[i]-v[j].y)/(v[(j+1)%n].y-v[j].y);
				L=min(L,px);R=max(R,px);
			}
		}
		double t1=R-L;
		L=99999999;
		R=-99999999;
		for(int j=0;j<n;j++){
			if(min(v[j].y,v[(j+1)%n].y)<y[i+1]&&y[i+1]<max(v[j].y,v[(j+1)%n].y)){
				double px=v[j].x+(v[(j+1)%n].x-v[j].x)*(y[i+1]-v[j].y)/(v[(j+1)%n].y-v[j].y);
				L=min(L,px);R=max(R,px);
			}
		}
		double t2=R-L;
		double val=(t1*y[i]*y[i]+t2*y[i+1]*y[i+1]+(t1+t2)*(y[i]+y[i+1])*(y[i]+y[i+1])/2)*(y[i+1]-y[i])/6;
		ret+=val;
	}
	return ret;
}
int main(){
	int a,b;
	scanf("%d%d",&a,&b);
	vector<Pt>s;
	for(int i=0;i<a;i++){
		double x,y;scanf("%lf%lf",&x,&y);
		c[i]=Pt(x,y);
		s.push_back(c[i]);
	}
	double S=0;
	for(int i=0;i<s.size();i++){
		S+=s[i].det(s[(i+1)%s.size()]);
	}
	S=ABS(S)/2;
	for(int i=0;i<b;i++){
		double x,y;scanf("%lf%lf",&x,&y);
		d[i]=Pt(x,y);
	}
	//	s.push_back(c[0]);
	double val=0;
	for(int i=0;i<b;i++){
		vector<Pt> now=s;
		for(int j=0;j<b;j++){
			if(i==j)continue;
			pair<Pt,Pt> L=make_pair((d[i]+d[j])/2.0,(d[i]+d[j])/2.0+(d[i]-d[j])*Pt(0,1));
			now=convexCut(now.size(),now,L.second,L.first);
		}
		for(int j=0;j<now.size();j++)now[j]=now[j]-d[i];
		val+=calc(now);
	}
	printf("%.12f\n",val/S);
}