Problem: Example

Input

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


Output

2.0000000000

Solution (in CPP):

Solution:
#include <iostream>
#include <cstdio>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <deque>
#include <stack>
#include <algorithm>
#include <cstring>
#include <functional>
#include <cmath>
#include <complex>
using namespace std;
#define rep(i,n) for(int i=0;i<(int)n;++i)
#define rep1(i,n) for(int i=1;i<=(int)n;++i)
#define all(c) (c).begin(),(c).end()
#define fs first
#define sc second
#define pb push_back
#define show(x) cout << #x << " " << x << endl
#define X real()
#define Y imag()
typedef long double D;
typedef complex<D> P;
typedef pair<P,P> L;
typedef vector<P> Pol;
D eps=1e-9;
D cro(P a,P b){return imag(conj(a)*b);}
inline int ccw(P a,P b,P c){
	if(cro(b-a,c-a)>eps) return 1;
	if(cro(b-a,c-a)<-eps) return -1;
	return 0;
}

bool iLSex(L l,L s){
	return cro(l.sc-l.fs,s.fs-l.fs)*cro(l.sc-l.fs,s.sc-l.fs)<-eps;
}
P intLL(L a,L b){
	D t=cro(a.sc-a.fs,a.sc-b.fs)/cro(a.sc-a.fs,b.sc-b.fs);
	return b.fs+t*(b.sc-b.fs);
}
Pol convexcut(Pol p,L l){
	Pol ret;
	rep(i,p.size()){
		if(ccw(l.fs,l.sc,p[i])!=-1) ret.pb(p[i]);
		L s=L(p[i],p[(i+1)%p.size()]);
		if(iLSex(l,s)) ret.pb(intLL(l,s));
	}
	return ret;
}
vector<Pol> makevoronoi(Pol vp,Pol pol){
	vector<Pol> ret;
	Pol polc=pol;
	rep(i,vp.size()){
		pol=polc;
		rep(j,vp.size()){
			if(i==j) continue;
			P p=(vp[i]+vp[j])/(long double)2.0,q=p+(vp[j]-vp[i])*P(0,1);
			pol=convexcut(pol,L(p,q));
		}
		ret.pb(pol);
	}
	return ret;
}
D aPol(Pol p){
	int n=p.size();
	D ret=0;
	rep(i,n) ret+=cro(p[i],p[(i+1)%n])/2;
	return ret;
}
int M,N;
Pol pol,vp;
D calc(D e,D f,D g,D h,D x){
	return e*x*x*x*x+f*x*x*x+g*x*x+h*x;
}
int main(){
	cin>>M>>N;
	rep(i,M){
		int x,y;
		cin>>x>>y;
		pol.pb(P(x,y));
	}
	rep(i,N){
		int x,y;
		cin>>x>>y;
		vp.pb(P(x,y));
	}
	vector<Pol> vor=makevoronoi(vp,pol);
	D sum=0;
	rep(i,N){
		Pol v=vor[i];
		P o=vp[i];
		D tmp=0;
//		show(i);
//		rep(j,v.size()) cout<<j<<":"<<v[j]<<endl;
//		cout<<endl;
		rep(j,v.size()){
			P p=v[j],q=v[(j+1)%v.size()];
			if(abs(p.X-q.X)<eps) continue;
			D a=(p.Y-q.Y)/(p.X-q.X),b=p.Y-a*p.X;
			b+=a*o.X-o.Y;
			D c=-o.Y;
			D e=a+a*a*a/3,f=b-c+a*a*b,g=a*b*b,h=(b*b*b-c*c*c)/3;
			e/=4,f/=3,g/=2;
			D val=calc(e,f,g,h,p.X-o.X)-calc(e,f,g,h,q.X-o.X);
			tmp+=val;
		}
		sum+=tmp;
	}
	printf("%.12f\n",(double)sum/(double)aPol(pol));
}