diff --git a/docs/graph/general_weighted_matching.md b/docs/graph/general_weighted_matching.md new file mode 100644 index 0000000..a3308b2 --- /dev/null +++ b/docs/graph/general_weighted_matching.md @@ -0,0 +1,25 @@ +--- +title: General Weighted Matching +documentation_of: //graph/general_weighted_matching.hpp +--- + +一般グラフの最大重みマッチング(Blossom / 花アルゴリズム)を解く。 +重み付き辺を受け取り、重み和が最大となるマッチングを返す。 +辺重みは非負整数を想定する。 + +## 使い方 + +```cpp +mylib::GeneralWeightedMatching wm(N); +wm.add_edge(u, v, w); +auto [sum, mate] = wm.solve(); +``` + +- `sum`: 最大重み +- `mate[i]`: 頂点 `i` とマッチした頂点番号(未マッチは `-1`) + +同じ 2 頂点間に複数回 `add_edge` した場合は、最大の重みだけを使う。 + +## 計算量 + +- $O(N^3)$ diff --git a/graph/general_weighted_matching.hpp b/graph/general_weighted_matching.hpp new file mode 100644 index 0000000..eee20dd --- /dev/null +++ b/graph/general_weighted_matching.hpp @@ -0,0 +1,320 @@ +#ifndef HARUILIB_GRAPH_GENERAL_WEIGHTED_MATCHING_HPP +#define HARUILIB_GRAPH_GENERAL_WEIGHTED_MATCHING_HPP + +#include +#include +#include +#include +#include +#include + +namespace mylib { + +struct GeneralWeightedMatching { + int n; + std::vector> w; + + explicit GeneralWeightedMatching(int n_) : n(n_), w(n_, std::vector(n_, 0)) {} + + void add_edge(int u, int v, long long cost) { + if (u == v) return; + if (w[u][v] < cost) w[u][v] = w[v][u] = cost; + } + + std::pair> solve() const { + if (n == 0) return {0, {}}; + + struct Edge { + int u, v; + long long w; + Edge(int u_ = 0, int v_ = 0, long long w_ = 0) : u(u_), v(v_), w(w_) {} + }; + + const int N = n; + const int max_nodes = 2 * N; + int nodes = N; + + std::vector> g(max_nodes + 1, std::vector(max_nodes + 1)); + for (int i = 1; i <= max_nodes; i++) { + for (int j = 1; j <= max_nodes; j++) { + g[i][j] = Edge(i, j, 0); + } + } + for (int i = 1; i <= N; i++) { + for (int j = 1; j <= N; j++) { + g[i][j].w = w[i - 1][j - 1]; + } + } + + std::vector match(max_nodes + 1), slack(max_nodes + 1), st(max_nodes + 1), parent(max_nodes + 1); + std::vector state(max_nodes + 1), vis(max_nodes + 1); + std::vector label(max_nodes + 1); + std::vector> blossom_from(max_nodes + 1, std::vector(N + 1)); + std::vector> blossom(max_nodes + 1); + std::queue q; + int vis_timer = 0; + + auto delta = [&](const Edge &e) { + return label[e.u] + label[e.v] - e.w * 2; + }; + + auto update_slack = [&](int u, int x) { + if (!slack[x] || delta(g[u][x]) < delta(g[slack[x]][x])) { + slack[x] = u; + } + }; + + auto set_slack = [&](int x) { + slack[x] = 0; + for (int u = 1; u <= N; u++) { + if (g[u][x].w > 0 && st[u] != x && state[st[u]] == 0) update_slack(u, x); + } + }; + + std::function q_push = [&](int x) { + if (x <= N) { + q.push(x); + } else { + for (int v : blossom[x]) q_push(v); + } + }; + + std::function set_st = [&](int x, int b) { + st[x] = b; + if (x > N) { + for (int v : blossom[x]) set_st(v, b); + } + }; + + auto get_pr = [&](int b, int xr) { + int pr = std::find(blossom[b].begin(), blossom[b].end(), xr) - blossom[b].begin(); + if (pr % 2 == 1) { + std::reverse(blossom[b].begin() + 1, blossom[b].end()); + return (int)blossom[b].size() - pr; + } + return pr; + }; + + std::function set_match = [&](int u, int v) { + match[u] = g[u][v].v; + if (u <= N) return; + int xr = blossom_from[u][g[u][v].u]; + int pr = get_pr(u, xr); + for (int i = 0; i < pr; i++) { + set_match(blossom[u][i], blossom[u][i ^ 1]); + } + set_match(xr, v); + std::rotate(blossom[u].begin(), blossom[u].begin() + pr, blossom[u].end()); + }; + + auto augment = [&](int u, int v) { + while (true) { + int next = st[match[u]]; + set_match(u, v); + if (!next) return; + set_match(next, st[parent[next]]); + u = st[parent[next]]; + v = next; + } + }; + + auto lca = [&](int u, int v) { + ++vis_timer; + for (; u || v; std::swap(u, v)) { + if (!u) continue; + if (vis[u] == vis_timer) return u; + vis[u] = vis_timer; + u = st[match[u]]; + if (u) u = st[parent[u]]; + } + return 0; + }; + + auto add_blossom = [&](int u, int l, int v) { + int b = N + 1; + while (b <= nodes && st[b]) b++; + if (b > nodes) nodes++; + + label[b] = 0; + state[b] = 0; + match[b] = match[l]; + blossom[b].clear(); + blossom[b].push_back(l); + for (int x = u, y; x != l; x = st[parent[y]]) { + blossom[b].push_back(x); + y = st[match[x]]; + blossom[b].push_back(y); + q_push(y); + } + std::reverse(blossom[b].begin() + 1, blossom[b].end()); + for (int x = v, y; x != l; x = st[parent[y]]) { + blossom[b].push_back(x); + y = st[match[x]]; + blossom[b].push_back(y); + q_push(y); + } + set_st(b, b); + + for (int x = 1; x <= nodes; x++) { + g[b][x] = Edge(b, x, 0); + g[x][b] = Edge(x, b, 0); + } + for (int x = 1; x <= N; x++) blossom_from[b][x] = 0; + + for (int x : blossom[b]) { + for (int y = 1; y <= nodes; y++) { + if (g[b][y].w == 0 || delta(g[x][y]) < delta(g[b][y])) { + g[b][y] = g[x][y]; + g[y][b] = g[y][x]; + } + } + for (int y = 1; y <= N; y++) { + if (blossom_from[x][y]) blossom_from[b][y] = x; + } + } + set_slack(b); + }; + + auto expand_blossom = [&](int b) { + for (int x : blossom[b]) set_st(x, x); + int xr = blossom_from[b][g[b][parent[b]].u]; + int pr = get_pr(b, xr); + for (int i = 0; i < pr; i += 2) { + int xs = blossom[b][i], xns = blossom[b][i + 1]; + parent[xs] = g[xns][xs].u; + state[xs] = 1; + state[xns] = 0; + slack[xs] = 0; + set_slack(xns); + q_push(xns); + } + state[xr] = 1; + parent[xr] = parent[b]; + for (int i = pr + 1; i < (int)blossom[b].size(); i++) { + int x = blossom[b][i]; + state[x] = -1; + set_slack(x); + } + st[b] = 0; + }; + + auto on_found_edge = [&](const Edge &e) { + int u = st[e.u], v = st[e.v]; + if (state[v] == -1) { + parent[v] = e.u; + state[v] = 1; + int next = st[match[v]]; + slack[v] = slack[next] = 0; + state[next] = 0; + q_push(next); + } else if (state[v] == 0) { + int l = lca(u, v); + if (!l) { + augment(u, v); + augment(v, u); + return true; + } + add_blossom(u, l, v); + } + return false; + }; + + auto matching = [&]() { + std::fill(state.begin(), state.end(), -1); + std::fill(slack.begin(), slack.end(), 0); + q = std::queue(); + for (int x = 1; x <= nodes; x++) { + if (st[x] == x && !match[x]) { + parent[x] = 0; + state[x] = 0; + q_push(x); + } + } + if (q.empty()) return false; + + while (true) { + while (!q.empty()) { + int u = q.front(); + q.pop(); + if (state[st[u]] == 1) continue; + for (int v = 1; v <= N; v++) { + if (g[u][v].w > 0 && st[u] != st[v]) { + if (delta(g[u][v]) == 0) { + if (on_found_edge(g[u][v])) return true; + } else { + update_slack(u, st[v]); + } + } + } + } + + long long d = std::numeric_limits::max(); + for (int u = 1; u <= N; u++) { + if (state[st[u]] == 0) d = std::min(d, label[u]); + } + for (int b = N + 1; b <= nodes; b++) { + if (st[b] == b && state[b] == 1) d = std::min(d, label[b] / 2); + } + for (int x = 1; x <= nodes; x++) { + if (st[x] == x && slack[x]) { + if (state[x] == -1) d = std::min(d, delta(g[slack[x]][x])); + else if (state[x] == 0) d = std::min(d, delta(g[slack[x]][x]) / 2); + } + } + if (d == std::numeric_limits::max()) return false; + + for (int u = 1; u <= N; u++) { + if (state[st[u]] == 0) { + if (label[u] == d) return false; + label[u] -= d; + } else if (state[st[u]] == 1) { + label[u] += d; + } + } + for (int b = N + 1; b <= nodes; b++) { + if (st[b] == b) { + if (state[b] == 0) label[b] += d * 2; + else if (state[b] == 1) label[b] -= d * 2; + } + } + + q = std::queue(); + for (int x = 1; x <= nodes; x++) { + if (st[x] == x && slack[x] && st[slack[x]] != x && delta(g[slack[x]][x]) == 0) { + if (on_found_edge(g[slack[x]][x])) return true; + } + } + for (int b = N + 1; b <= nodes; b++) { + if (st[b] == b && state[b] == 1 && label[b] == 0) expand_blossom(b); + } + } + }; + + long long max_weight = 0; + for (int u = 1; u <= N; u++) { + st[u] = u; + blossom[u].clear(); + for (int v = 1; v <= N; v++) { + blossom_from[u][v] = (u == v ? u : 0); + max_weight = std::max(max_weight, g[u][v].w); + } + } + for (int u = 1; u <= N; u++) label[u] = max_weight; + + while (matching()) {} + + std::vector mate(n, -1); + long long total = 0; + for (int i = 1; i <= N; i++) { + if (match[i] && match[i] <= N) mate[i - 1] = match[i] - 1; + } + for (int i = 0; i < n; i++) { + if (mate[i] > i) total += w[i][mate[i]]; + } + return {total, mate}; + } +}; + +} // namespace mylib + +#endif diff --git a/test/verify/yosupo-general-weighted-matching.test.cpp b/test/verify/yosupo-general-weighted-matching.test.cpp new file mode 100644 index 0000000..922994e --- /dev/null +++ b/test/verify/yosupo-general-weighted-matching.test.cpp @@ -0,0 +1,29 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/general_weighted_matching" + +#include "template/template.hpp" +#include "graph/general_weighted_matching.hpp" + +int main() { + ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); + + int N, M; cin >> N >> M; + mylib::GeneralWeightedMatching wm(N); + for (int i = 0; i < M; i++) { + int u, v; + long long w; + cin >> u >> v >> w; + wm.add_edge(u, v, w); + } + + auto ans = wm.solve(); + vector> edges; + for (int i = 0; i < N; i++) { + if (ans.second[i] != -1 && i < ans.second[i]) edges.emplace_back(i, ans.second[i]); + } + + cout << edges.size() << " " << ans.first << "\n"; + for (auto [u, v] : edges) { + cout << u << " " << v << "\n"; + } + return 0; +}