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 repl(i,a,b) for(int i=(int)(a);i<(int)(b);i++)
#define rep(i,n) repl(i,0,n)
#define all(x) (x).begin(),(x).end()
#define dbg(x) cout<<#x"="<<x<<endl
#define mmax(x,y) (x>y?x:y)
#define mmin(x,y) (x<y?x:y)
#define maxch(x,y) x=mmax(x,y)
#define minch(x,y) x=mmin(x,y)
#define uni(x) x.erase(unique(all(x)),x.end())
#define exist(x,y) (find(all(x),y)!=x.end())
#define bcnt __builtin_popcountll

typedef double D;
typedef complex<D> P;

#define X real()
#define Y imag()

const D eps=1e-8;
const D inf=1e12;
struct L : public vector<P> { // line and segment
  L(const P& a,const P &b){
    push_back(a);
    push_back(b);
  }
};
D dot(P a,P b){ return (conj(a)*b).X; }
D cross(P a,P b){ return (conj(a)*b).Y; }
P crosspoint(L l,L m) {
  D A=cross(l[1]-l[0],m[1]-m[0]);
  D B=cross(l[1]-l[0],l[1]-m[0]);
  if (abs(A)<eps&&abs(B)<eps) return m[0]; // same line
  if (abs(A)<eps) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
  return m[0]+B/A*(m[1]-m[0]);
}
int ccw(P a,P b,P c) {
  b=b-a;c=c-a;
  if(cross(b,c)>0) return +1;       // counter clockwise
  if(cross(b,c)<0) return -1;       // clockwise
  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;
}
vector<P> convex_cut(const vector<P>& poly,L l) {
  vector<P> Q;
  int sz=poly.size();
  for (int i = 0; i < sz; i++) {
    P A = poly[i], B = poly[(i+1)%sz];
    if(ccw(l[0], l[1], A) != -1)Q.push_back(A);
    if(ccw(l[0], l[1], A)*ccw(l[0], l[1], B) < 0){
      Q.push_back(crosspoint(L(A, B), l));
    }
  }
  return Q;
}
L half_line(P a,P b){
  P m=(a+b)*(D)0.5;
  P dir=b-a;
  return L(m,m+P(-dir.Y,dir.X));
}
D calc_area(vector<P>& poly) {
  D res = 0;
  int sz=poly.size();
  for (int i = 0; i < sz; i++){
    res += cross(poly[i], poly[(i+1)%sz]);
  }
  return res/2.0;
}

int M,N;
vector<P> poly,ps;

int main(){
  cin>>M>>N;
  rep(i,M){
    int x,y;
    cin>>x>>y;
    poly.push_back(P(x,y));
  }
  rep(i,N){
    int x,y;
    cin>>x>>y;
    ps.push_back(P(x,y));
  }
  D res=0.0;
  rep(i,N){
    vector<P> crtpoly=poly;
    rep(j,N){
      if(i==j)continue;
      L h=half_line(ps[i],ps[j]);
      if(ccw(h[0],h[1],ps[i])==-1)swap(h[0],h[1]);
      crtpoly=convex_cut(crtpoly,h);
    }
    assert(crtpoly.size()>2);
    D totarea=calc_area(poly);
    rep(j,crtpoly.size()){
      P p0=crtpoly[j]-ps[i],p1=crtpoly[(j+1)%crtpoly.size()]-ps[i];
      vector<P> tri;
      tri.push_back(P(0,0)); tri.push_back(p0); tri.push_back(p1);
      D area=calc_area(tri);
      D p=p1.X-p0.X,q=p1.Y-p0.Y;
      D a=p*p+q*q, b=2*(p*p0.X+q*p0.Y), c=p0.X*p0.X+p0.Y*p0.Y;
      res+=((2*a+3*b+6*c)/12.0)*(area/totarea);
    }
  }
  printf("%.10f\n", res);

  return 0;
}