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 <string.h>
#include <algorithm>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <vector>
#include <queue>
#include <string>
#include <map>
#include <set>

using namespace std;
typedef long long ll;
typedef unsigned int uint;
typedef unsigned long long ull;
static const long double EPS = 1e-9;
static const long double PI = acos(-1.0);

#define REP(i, n) for (int i = 0; i < (int)(n); i++)
#define FOR(i, s, n) for (int i = (s); i < (int)(n); i++)
#define FOREQ(i, s, n) for (int i = (s); i <= (int)(n); i++)
#define FORIT(it, c) for (__typeof((c).begin())it = (c).begin(); it != (c).end(); it++)
#define MEMSET(v, h) memset((v), h, sizeof(v))

#include <complex>

typedef complex<long double> Point;
typedef vector<Point> Polygon;

static const long double INF = 1e+5;

#define CURR(P, i) (P[i])
#define NEXT(P, i) (P[(i + 1) % P.size()])

struct Line : public vector<Point> {
  Line() {;}
  Line(Point a, Point b) { push_back(a); push_back(b); }
};

struct Circle {
  Point p;
  long double r;
  Circle() {;}
  Circle(Point p, long double r) : p(p), r(r) {;}
};

namespace std {
  bool operator<(const Point &lhs, const Point &rhs) {
    return lhs.real() == rhs.real() ? lhs.imag() < rhs.imag() : lhs.real() < rhs.real();
  }
}

inline long double cross(const Point &a, const Point &b) {
  return imag(conj(a) * b);
}

inline long double dot(const Point &a, const Point &b) {
  return real(conj(a) * b);
}

inline int ccw(Point a, Point b, Point c) {
  b -= a;
  c -= a;
  if (cross(b, c) > 0) { return 1; }
  if (cross(b, c) < 0) { return -1; }
  if (dot(b, c) < 0) { return 2; }
  if (norm(b) < norm(c)) { return -2; }
  return 0;
}


bool intersectSS(const Line &s, const Line &t) {
  return ccw(s[0], s[1], t[0]) * ccw(s[0], s[1], t[1]) <= 0 &&
    ccw(t[0], t[1], s[0]) * ccw(t[0], t[1], s[1]) <= 0;
}

bool intersectSP(const Line &s, const Point &p) {
  return abs(s[0] - p) + abs(s[1] - p) - abs(s[1] - s[0]) < EPS;
}

Point crosspointSS(const Line &l, const Line &m) {
  long double A = cross(l[1] - l[0], m[1] - m[0]);
  long double B = cross(l[1] - l[0], l[1] - m[0]);
  if (abs(A) < EPS && abs(B) < EPS) {
    assert(false);
    if (intersectSP(l, m[0])) { return m[0]; }
    if (intersectSP(l, m[1])) { return m[1]; }
    if (intersectSP(m, l[0])) { return l[0]; }
    if (intersectSP(m, l[1])) { return l[1]; }
    return m[0];
  }
  if (abs(A) < EPS) { assert(false); }
  return m[0] + B / A * (m[1] - m[0]);
}

Polygon ConvexCut(const Polygon &P, const Line &l) {
  Polygon Q;
  for (int i = 0; i < (int)P.size(); i++) {
    Point A = CURR(P, i), B = NEXT(P, i);
    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(crosspointSS(Line(A, B), l));
    }
  }
  return Q;
}

Line Bisect(const Point &p, const Point &q) {
  Point center = (p + q) / 2.0L;
  Point vect = (q - p) * Point(0, 1);
  vect /= abs(vect);
  return Line(center - vect * INF, center + vect * INF);
}

Polygon VoronoiCut(Polygon poly, const vector<Point> &points, int index) {
  const int n = points.size();
  REP(i, n) {
    if (i == index) { continue; }
    Line l = Bisect(points[index], points[i]);
    poly = ConvexCut(poly, l);
  }
  return poly;
}

long double Area(const Polygon &poly) {
  const int n = poly.size();
  long double ret = 0.0;
  REP(i, n) {
    Point prev = poly[i];
    Point next = poly[(i + 1) % n];
    ret += cross(prev, next);
  }
  return ret / 2.0;
}

long double V(long double a, long double b, long double x) {
  long double x2 = x * x;
  long double x3 = x * x * x;
  long double x4 = x * x * x * x;
  return 1.0 / 12.0 * a * (3 + a * a) * x4 + 1.0 / 3.0 * b * (1 + a * a) * x3 + 1.0 / 2.0 * a * b * b * x2 + 1.0 / 3.0 * b * b * b * x;
}

long double Int(long double a, long double b, long double l, long double r) {
  assert(l <= r);
  return V(a, b, r) - V(a, b, l);
}

long double calc(Point p, Point q) {
  if (fabs(p.real() - q.real()) < EPS) { return 0; }
  assert(fabs(p.imag()) < EPS || fabs(q.imag()) < EPS);
  if (p.real() < q.real()) { swap(p, q); }
  long double a = (p.imag() - q.imag()) / (p.real() - q.real());
  long double b = p.imag() - a * p.real();
  return Int(a, b, q.real(), p.real());
}

Polygon ini;
vector<Point> shelters;

int m, n;
int main() {
  while (scanf("%d %d", &m, &n) > 0) {
    long double ans = 0.0;
    ini.clear();
    shelters.clear();
    REP(i, m) {
      long double x, y;
      scanf("%Lf %Lf", &x, &y);
      ini.push_back(Point(x, y));
    }
    REP(i, n) {
      long double x, y;
      scanf("%Lf %Lf", &x, &y);
      shelters.push_back(Point(x, y));
    }
    REP(i, n) {
      Polygon voronoi = VoronoiCut(ini, shelters, i);
      int k = voronoi.size();
      long double lans  = 0.0;
      //REP(j, k) {
      //  cout << voronoi[j] << " ";
      //}
      //puts("");
      REP(j, k) { voronoi[j] -= shelters[i]; }
      REP(j, k) {
        int prev = j;
        int next = (j + 1) % k;
        long double theta = -arg(voronoi[prev]);
        Point p = voronoi[prev] * Point(cos(theta), sin(theta));
        Point q = voronoi[next] * Point(cos(theta), sin(theta));
        long double sum = 0.0;
        if (q.real() < -EPS) {
          sum -= calc(Point(0, 0), q);
        } else if (q.real() > EPS) {
          sum += calc(Point(0, 0), q);
        }
        if (q.real() > p.real()) {
          sum -= calc(p, q);
        } else {
          sum += calc(p, q);
        }
        lans += sum;
      }
      ans += lans;
    }
    printf("%.8Lf\n", ans / Area(ini));
  }
}