[SCOI2018]游泳池（计算几何＋分数规划＋最大权闭合子图）

题目链接

https://www.luogu.org/problemnew/show/U56187

注：题面参考了网上的其他博客，并非原题题面，因此数据范围可能有误。数据为原创数据。

题解

其实就是许多板子码到一起。

首先对于边缘上的任意一点 (u)，假设离它最远的顶点为 (A)，那么我们称点 (u) 位于顶点 (A) 的控制范围之中。我们考虑在没有石雕的情况下怎么求出每个顶点的控制范围。对于除顶点 (A) 之外的任意一个顶点 (B)，连接 (AB) 并作 (AB) 的中垂线 (l)，显然若仅考虑顶点 (A) 与顶点 (B)，那么直线 (l) 两侧中距顶点 (A) 较远的那一侧内的点会位于顶点 (A) 的控制范围之中。因此，若需要在考虑其余所有顶点的情况下求出顶点 (A) 的控制范围，只需要对其余每一个顶点都求出顶点 (A) 与其连线的中垂线，然后对所有直线与凸多边形本身求半平面交即可。当有石雕时，我们还需要求出每个点对应的控制范围中被石雕遮挡的范围。对于每个顶点 (A)，我们可以直接暴力枚举石雕，找到顶点 (A) 与石雕对应的圆的两条切线，那么切线之间的范围就是顶点 (A) 不能直接到达的范围。我们对所有石雕遮挡的范围求交后对顶点 (A) 本身的控制范围求并即可，这一步可以通过对每个顶点做一次扫描线来实现。

接下来考虑会移除石雕的情况。

首先，求答案式 (frac{E}{t + kx}) 的最值显然可以使用分数规划。设原本在没有移走任何石雕的情况下，所有权值为 (1) 的顶点的未被石雕遮挡的控制范围占凸多边形总边长的比值为 (p_0)，那么显然有 (E = p_0t)。如果移走了 (x) 个石雕后，新增加的控制范围占泳池总边长的比值为 (p')，那么显然新的期望收益 (E' = (p_0 + p')t)，故我们二分答案 ({ m ans}) 后只需判断是否存在一种移除方式，满足移除后有 (frac{(p_0 + p')t}{t + kx} geq { m ans})，即 (p_0t + p't - { m ans}cdot t - kxcdot { m ans} geq 0)。当 ({ m ans}) 确定时，(p_0t - { m ans} cdot t) 为定值，因此我们只需考虑求 (p't - kxcdot { m ans}) 的最大值。

我们将式子 (p't - kxcdot { m ans}) 中的两部分分开考虑，整个问题其实可视为一个最大权闭合子图模型。具体地，若我们将一个石雕看做权值为 (-k cdot { m ans}) 的结点，将一段被若干石雕覆盖的区域看做是权值为 (p_xt)（(p_x) 为该段区域的长度占凸多边形总边长的比值）的结点，由于凸多边形上的一段区域可能同时被多个石雕所遮挡，那么选择一段被覆盖的区域就意味着要选择覆盖该区域的所有石雕。我们的目的就是保证合法的情况下选择出权值和最大的若干个点。因此，我们可以建图跑最小割来判断 ({ m ans}) 的合法性。对于如何求出每一段区域被哪些石雕覆盖，我们只需在最开始做扫描线时一并记录即可。

时间复杂度为 (O(n^2m + (n(n + m) + { m maxflow})log w))，但常数较大。

代码

我的实现略显复杂。

#include<bits/stdc++.h>

using namespace std;

const int N = 510, max_node = 3e4 + 10;
const double eps = 1e-10, inf = 1e9;

struct point {
  int id;
  double x, y;

  point () {}
  point (double x, double y): x(x), y(y) {}

  point operator + (point a) {
    return point (x + a.x, y + a.y);
  }

  point operator - (point a) {
    return point (x - a.x, y - a.y);
  }

  point operator * (double a) {
    return point (x * a, y * a);
  }

  point operator / (double a) {
    return point (x / a, y / a);
  }

  bool operator == (const point& a) const {
    return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
  }

  bool operator != (const point& a) const {
    return !(*this == a);
  }
} points[N], point_q[N << 1];

double dot(point a, point b) {
  return a.x * b.x + a.y * b.y;
}

double cross(point a, point b) {
  return a.x * b.y - a.y * b.x;
}

double dist(point a, point b) {
  return sqrt(dot(a - b, a - b));
}

point rotate(point a, double rad) {
  return point (a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
}

bool on_segment(point p, point a, point b) {
  double d = dist(a, b);
  return dist(a, p) <= d && dist(b, p) <= d;
}

struct circle {
  point center;
  double r;

  circle () {}
  circle (point center, double r): center(center), r(r) {}
} stone[N];

struct line {
  int id;
  point p, v;
  double angle;

  line () {}
  line (point p, point v, int id): p(p), v(v), id(id) {
    angle = atan2(v.y, v.x);
  }

  bool operator < (const line& a) const {
    return angle < a.angle;
  }
} line_q[N << 1];

int n, m, days, cost, type[N];
double p0, all_dist;
vector<pair<point, point>> control[N][N];

point line_intersection(line a, line b) {
  if (fabs(cross(a.v, b.v)) > eps) {
    double t = cross(b.v, a.p - b.p) / cross(a.v, b.v);
    return a.p + a.v * t;
  } else {
    return point (inf, inf);
  }
}

bool on_left(point p, line x) {
  return cross(x.v, p - x.p) > 0;
}

vector<line> half_plane_intersection(vector<line>& s) {
  sort(s.begin(), s.end());
  int first = 0, last = 0;
  line_q[0] = s[0];
  for (int i = 1; i < s.size(); ++i) {
    for (; first < last && !on_left(point_q[last - 1], s[i]); --last);
    for (; first < last && !on_left(point_q[first], s[i]); ++first);
    line_q[++last] = s[i];
    if (first < last && fabs(cross(line_q[last].v, line_q[last - 1].v)) < eps) {
      --last;
      if (on_left(s[i].p, line_q[last])) {
        line_q[last] = s[i];
      }
    }
    if (first < last) {
      point_q[last - 1] = line_intersection(line_q[last - 1], line_q[last]);
    }
  }
  for (; first < last && !on_left(point_q[last - 1], line_q[first]); --last);
  vector<line> result;
  for (int i = first; i <= last; ++i) {
    result.push_back(line_q[i]);
  }
  return result;
}

struct edge {
  int to;
  double refer, cap, flow;

  edge () {}
  edge (int to, double refer): to(to), refer(refer) {}
};

vector<edge> edges;
vector<int> graph[max_node];
int s, t, node_cnt, stone_id[N], cur[max_node], dis[max_node];
double dinic_flow;

void add_edge(int u, int v, double cap) {
  edges.emplace_back(v, cap);
  edges.emplace_back(u, 0);
  graph[u].push_back(edges.size() - 2);
  graph[v].push_back(edges.size() - 1);
}

bool bfs() {
  queue<int> que;
  memset(dis, 0, sizeof dis);
  memset(cur, 0, sizeof cur);
  que.push(s);
  dis[s] = 1;
  while (!que.empty()) {
    int x = que.front();
    que.pop();
    for (auto v : graph[x]) {
      edge& e = edges[v];
      if (e.cap - e.flow > eps && !dis[e.to]) {
        dis[e.to] = dis[x] + 1;
        que.push(e.to);
      }
    }
  }
  return dis[t] > 0;
}

double dfs(int u, double a) {
  if (u == t || a < eps) {
    return a;
  } else {
    double flow = 0, f;
    for (int& i = cur[u]; i < graph[u].size(); ++i) {
      edge& e = edges[graph[u][i]];
      if (dis[e.to] == dis[u] + 1 && (f = dfs(e.to, min(a, e.cap - e.flow))) > eps) {
        e.flow += f;
        edges[graph[u][i] ^ 1].flow -= f;
        flow += f;
        a -= f;
        if (a < eps) {
          break;
        }
      }
    }
    return flow;
  }
}

double dinic() {
  for (; bfs(); dinic_flow -= dfs(s, inf));
  return dinic_flow;
}

struct state {
  int p_id, stone_id;
  double d;
  point p;

  state () {}
  state (int p_id, int stone_id, double d, point p): p_id(p_id), stone_id(stone_id), d(d), p(p) {}

  bool operator < (const state& a) const {
    return p_id == a.p_id ? d < a.d : p_id < a.p_id;
  }
};

vector<state> cover[N];
bool appeared[N];

void add_edge(int i, int j, point a, point b) {
  if (control[i][j].size()) {
    point p = control[i][j][0].first;
    point q = control[i][j][0].second;
    double dist0 = dist(points[j], p);
    double dist1 = dist(points[j], a);
    double dist2 = dist(points[j], q);
    double dist3 = dist(points[j], b);
    if (max(dist0, dist1) <= min(dist2, dist3)) {
      point l = dist0 > dist1 ? p : a;
      point r = dist2 < dist3 ? q : b;
      ++node_cnt;
      double d = dist(l, r) / all_dist;
      p0 -= d;
      add_edge(s, node_cnt, d * days);
      dinic_flow += d * days;
      for (int k = 1; k <= m; ++k) {
        if (appeared[k]) {
          add_edge(node_cnt, stone_id[k], inf);
        }
      }
    }
  }
}

void add_edge() {
  s = ++node_cnt;
  t = ++node_cnt;
  for (int i = 1; i <= m; ++i) {
    stone_id[i] = ++node_cnt;
    add_edge(stone_id[i], t, -cost);
  }
  for (int i = 1; i <= n; ++i) {
    if (type[i]) {
      memset(appeared, false, sizeof appeared);
      for (int j = 0, k, last = -1; j < cover[i].size(); j = k) {
        if (~last) {
          for (int l = last == n ? 1 : last + 1; l != cover[i][j].p_id; l = l == n ? 1 : l + 1) {
            add_edge(i, l, points[l], points[l == n ? 1 : l + 1]);
          }
        }
        for (k = j; k < cover[i].size() && cover[i][j].p_id == cover[i][k].p_id; ++k);
        last = cover[i][j].p_id;
        add_edge(i, last, points[last], cover[i][j].p);
        for (int l = j; l < k; ++l) {
          int u = cover[i][l].stone_id;
          if (!appeared[u]) {
            appeared[u] = true;
          } else {
            appeared[u] = false;
          }
          if (l + 1 < k) {
            add_edge(i, last, cover[i][l].p, cover[i][l + 1].p);
          }
        }
        add_edge(i, last, cover[i][k - 1].p, points[last == n ? 1 : last + 1]);
      }
    }
  }
}

bool check(double answer) {
  dinic_flow = 0;
  for (auto& v : edges) {
    v.flow = 0;
    if (v.refer >= 0) {
      v.cap = v.refer;
      if (v.refer < inf) {
        dinic_flow += v.refer;
      }
    } else {
      v.cap = -answer * v.refer;
    }
  }
  return (p0 - answer) * days + dinic() >= 0;
}

int main() {
  scanf("%d%d%d%d", &n, &m, &days, &cost);
  for (int i = 1; i <= n; ++i) {
    scanf("%lf%lf%d", &points[i].x, &points[i].y, &type[i]);
  }
  for (int i = 1; i <= n; ++i) {
    if (type[i]) {
      vector<line> half_plane;
      for (int j = 1; j <= n; ++j) {
        int k = j == n ? 1 : j + 1;
        half_plane.emplace_back(points[j], points[k] - points[j], j);
        if (i != j) {
          point mid = point ((points[i].x + points[j].x) / 2, (points[i].y + points[j].y) / 2);
          point v = point (points[j].y - points[i].y, points[i].x - points[j].x);
          half_plane.emplace_back(mid, v, -1);
        }
      }
      vector<line> convex_hull = half_plane_intersection(half_plane);
      if (convex_hull.size() > 2) {
        for (int j = 0; j < convex_hull.size(); ++j) {
          if (~convex_hull[j].id) {
            int l = j == 0 ? convex_hull.size() - 1 : j - 1;
            int r = j + 1 == convex_hull.size() ? 0 : j + 1;
            point a = line_intersection(convex_hull[l], convex_hull[j]);
            point b = line_intersection(convex_hull[j], convex_hull[r]);
            p0 += dist(a, b);
            control[i][convex_hull[j].id].emplace_back(a, b);
          }
        }
      }
    }
  }
  for (int i = 1; i <= n; ++i) {
    int j = i == n ? 1 : i + 1;
    all_dist += dist(points[i], points[j]);
  }
  p0 /= all_dist;
  for (int i = 1; i <= m; ++i) {
    scanf("%lf%lf%lf", &stone[i].center.x, &stone[i].center.y, &stone[i].r);
  }
  for (int i = 1; i <= n; ++i) {
    if (type[i]) {
      for (int j = 1; j <= m; ++j) {
        point v = stone[j].center - points[i];
        double angle = asin(stone[j].r / dist(points[i], stone[j].center));
        line a = line (points[i], rotate(v, angle), -1);
        line b = line (points[i], rotate(v, -angle), -1);
        for (int k = 1; k <= n; ++k) {
          int p = k == n ? 1 : k + 1;
          point intersection = line_intersection(line(points[k], points[p] - points[k], -1), a);
          if (intersection != points[i] && on_segment(intersection, points[k], points[p])) {
            cover[i].emplace_back(k < i ? k + n : k, j, dist(points[k], intersection), intersection);
          }
          intersection = line_intersection(line(points[k], points[p] - points[k], -1), b);
          if (intersection != points[i] && on_segment(intersection, points[k], points[p])) {
            cover[i].emplace_back(k < i ? k + n : k, j, dist(points[k], intersection), intersection);
          }
        }
      }
      sort(cover[i].begin(), cover[i].end());
      for (auto& v : cover[i]) {
        if (v.p_id > n) {
          v.p_id -= n;
        }
      }
    }
  }
  add_edge();
  double l = 0, r = 1;
  while (r - l > 1e-5) {
    double mid = (l + r) / 2;
    if (check(mid)) {
      l = mid;
    } else {
      r = mid;
    }
  }
  printf("%.4lf
", l);
  return 0;
}