[SCOI2018]游泳池(計算幾何+分數規劃+最大權閉合子圖)
題目連結
https://www.luogu.org/problemnew/show/U56187
注:題面參考了網上的其他部落格,並非原題題面,因此資料範圍可能有誤。資料為原創資料。
題解
其實就是許多板子碼到一起。
首先對於邊緣上的任意一點 \(u\),假設離它最遠的頂點為 \(A\),那麼我們稱點 \(u\) 位於頂點 \(A\) 的控制範圍之中。我們考慮在沒有石雕的情況下怎麼求出每個頂點的控制範圍。對於除頂點 \(A\) 之外的任意一個頂點 \(B\),連線 \(AB\) 並作 \(AB\) 的中垂線 \(l\),顯然若僅考慮頂點 \(A\) 與頂點 \(B\),那麼直線 \(l\) 兩側中距頂點 \(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\)
我們將式子 \(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;
}