1. 程式人生 > >[SCOI2018]游泳池(計算幾何+分數規劃+最大權閉合子圖)

[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\),故我們二分答案 \(ans\) 後只需判斷是否存在一種移除方式,滿足移除後有 \(\frac{(p_0 + p')t}{t + kx} \geq ans\),即 \(p_0t + p't - ans\cdot t - kx\cdot ans \geq 0\)

。當 \(ans\) 確定時,\(p_0t - ans \cdot t\) 為定值,因此我們只需考慮求 \(p't - kx\cdot ans\) 的最大值。

我們將式子 \(p't - kx\cdot ans\) 中的兩部分分開考慮,整個問題其實可視為一個最大權閉合子圖模型。具體地,若我們將一個石雕看做權值為 \(-k \cdot ans\) 的結點,將一段被若干石雕覆蓋的區域看做是權值為 \(p_xt\)\(p_x\) 為該段區域的長度佔凸多邊形總邊長的比值)的結點,由於凸多邊形上的一段區域可能同時被多個石雕所遮擋,那麼選擇一段被覆蓋的區域就意味著要選擇覆蓋該區域的所有石雕。我們的目的就是保證合法的情況下選擇出權值和最大的若干個點。因此,我們可以建圖跑最小割來判斷 \(ans\) 的合法性。對於如何求出每一段區域被哪些石雕覆蓋,我們只需在最開始做掃描線時一併記錄即可。

時間複雜度約為 \(O(n^2m\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 cap, flow;

  edge () {}
  edge (int to, double cap, double flow): to(to), cap(cap), flow(flow) {}
};

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 reset() {
  edges.clear();
  dinic_flow = 0;
  for (int i = 1; i < max_node; ++i) {
    graph[i].clear();
  }
}

void add_edge(int u, int v, double cap) {
  edges.emplace_back(v, cap, 0);
  edges.emplace_back(u, 0, 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, bool get_dist) {
  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;
      if (get_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);
        }
      }
    }
  }
}

bool check(double answer, bool get_dist) {
  reset();
  node_cnt = 0;
  s = ++node_cnt;
  t = ++node_cnt;
  for (int i = 1; i <= m; ++i) {
    stone_id[i] = ++node_cnt;
    add_edge(stone_id[i], t, answer * 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], get_dist);
          }
        }
        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, get_dist);
        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, get_dist);
          }
        }
        add_edge(i, last, cover[i][k - 1].p, points[last == n ? 1 : last + 1], get_dist);
      }
    }
  }
  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;
        }
      }
    }
  }
  check(0, true);
  double l = 0, r = 1;
  while (r - l > 1e-5) {
    double mid = (l + r) / 2;
    if (check(mid, false)) {
      l = mid;
    } else {
      r = mid;
    }
  }
  printf("%.4lf\n", l);
  return 0;
}