#pragma GCC optimize ("O3") #pragma GCC target ("avx") #include #include #include #include #include #include #include #include #include #include #include #include #define _rep(_1, _2, _3, _4, name, ...) name #define rep2(i, n) rep3(i, 0, n) #define rep3(i, a, b) rep4(i, a, b, 1) #define rep4(i, a, b, c) for (int i = int(a); i < int(b); i += int(c)) #define rep(...) _rep(__VA_ARGS__, rep4, rep3, rep2, _)(__VA_ARGS__) using namespace std; using i64 = long long; using u32 = unsigned; using u64 = unsigned long long; using f80 = long double; using i8 = signed char; template class MaximumWeightedMatching { /* Maximum Weighted Matching in General Graphs - O(n^3) time - O(n + m) space Note: each vertex is 1-indexed. */ public: using weight_t = WeightType; using weight_sum_t = i64; struct Edge { int from, to; weight_t weight; }; private: enum TreeLabelNumber { INNER = -1, UNUSED = 0, OUTER = 1 }; enum LabelNumber { SEPARATED = -2, DEFAULT = -1 }; enum EdgeNumber { UNDEFINED = 1 << 30 }; static constexpr weight_t INF = weight_t(1) << (sizeof(weight_t) * 8 - 2); struct Node { int next, from, to; }; struct Label { int from, to; }; struct LinkedList { int eid, next; }; class Queue { public: Queue() {} Queue(int N) : que(N), qh(0), qt(0) {} void clear() { qh = qt = 0; } int* data() { return que.data(); } bool empty() const { return qh == qt; } int dequeue() { return que[qh++]; } void enqueue(int u) { que[qt++] = u; } int operator [] (int i) const { return que[i]; } int size() const { return qt; } vector que; int qh, qt; }; public: MaximumWeightedMatching(int N, const vector& raw_edges) : N(N), B((N - 1) / 2), size(N + B + 1) { offsets.assign(N + 2, 0); for (auto& e : raw_edges) { offsets[e.from + 1]++; offsets[e.to + 1]++; } rep(i, 1, N + 1) offsets[i] += offsets[i - 1]; edges.resize(raw_edges.size() * 2); rep(i, raw_edges.size()) { auto& e = raw_edges[i]; edges[offsets[e.from]++] = {e.from, e.to, 2 * e.weight}; edges[offsets[e.to]++] = {e.to, e.from, 2 * e.weight}; } rep(i, N + 1) offsets[N + 1 - i] = offsets[N - i]; offsets[0] = 0; } weight_sum_t maximum_weighted_matching() { initialize(); set_potential(); rep(u, 1, N + 1) if (!mate[u]) { for (int s = 0; !augmented(u, s); s = adjust_dual_solutions()); fix_blossom_bases(); clear_label(); } weight_sum_t ret = 0; rep(u, 1, N + 1) if (mate[u] > u) { weight_t max_w = 0; rep(eid, offsets[u], offsets[u + 1]) if (edges[eid].to == mate[u]) { max_w = max(max_w, edges[eid].weight); } ret += max_w; } return ret >> 1; } private: inline int encode(int e) const { return e + size + 1; // should be >= 3 } inline weight_t reduced_cost(int u, int v, const Edge& e) const { return potential[u] + potential[v] - e.weight; } inline weight_t reduced_cost(int eid) const { return reduced_cost(edges[eid].from, edges[eid].to, edges[eid]); } void rematch(int v, int w) { auto t = mate[v]; mate[v] = w; if (mate[t] != v) return; if (label[v].to == 0) { mate[t] = label[v].from; rematch(mate[t], t); } else { int x = label[v].from, y = label[v].to; rematch(x, y); rematch(y, x); } } Label search_blossom_edge(int bid) const { int b = base[bid], bv = b; for (; node[bv].next != b; bv = node[node[bv].next].next); return {node[bv].from, node[bv].to}; } void label_blossom(int bid, int m, Label l) { label[bid] = {l.from, (l.to == surface[l.to]) ? 0 : l.to}; if (bid <= N) return; int b = base[bid]; label_blossom(b, mate[bid] = m, l); l = search_blossom_edge(bid); for (int bv = b, bw; node[bv].next != b; bv = node[bw].next) { label_blossom(bw = node[bv].next, 0, l); label_blossom(node[bw].next, node[bw].from, {node[bv].from, node[bv].to}); } } int find_mate(int bid) { return bid <= N ? mate[bid] : mate[bid] = find_mate(base[bid]); } void push_inner_blossom_rec(int bid, bool push=true) { tree_label[bid] = (bid <= N) ? INNER : UNUSED; if (bid > N) { int v = base[bid], u = v; do { push_inner_blossom_rec(v, push); } while ( (v = node[v].next) != u); } else if (push) inner_vertices[inner_vertices_size++] = bid; } void push_inner_blossom(int bid) { if (tree_label[bid] != UNUSED) return; bool push = label[bid].from != SEPARATED; if (bid > N) { if (push) inner_blossoms[inner_blossom_size++] = bid; push_inner_blossom_rec(bid, push); } else if (push) inner_vertices[inner_vertices_size++] = bid; tree_label[bid] = INNER; } void push_outer_blossom_rec(int bid) { tree_label[bid] = (bid <= N) ? OUTER : UNUSED; if (bid > N) { int v = base[bid], u = v; do { push_outer_blossom_rec(v); } while ( (v = node[v].next) != u ); } else outer_vertices.enqueue(bid); } void push_outer_blossom(int bid, bool push) { push_outer_blossom_rec(bid); if (bid <= N) return; if (push) outer_blossoms[outer_blossom_size++] = bid, tree_label[bid] = OUTER; else tree_label[bid] = UNUSED; } inline void merge_edge(int x, int bx, int eid) { auto& e = edges[eid]; int y = e.to, by = surface[y]; if (tree_label[by] != OUTER || bx == by) return; auto r_cost = reduced_cost(x, y, e); if (r_cost < best_cost[by].first) { if (best_cost[by].first == INF) merged_edges[merged_edge_size++] = by; best_cost[by] = {r_cost, eid}; } } inline void merge_vertex(int x, int bx) { rep(eid, offsets[x], offsets[x + 1]) merge_edge(x, bx, eid); best_edge[x] = UNDEFINED; } void clear_best_edges(int b) { if (b > N) { int v = b = base[b]; do { clear_best_edges(v); } while ( (v = node[v].next) != b ); } else best_edge[b] = UNDEFINED; } void merge_outer(int b, int bid) { if (b > N) { for (int head = blist_head[b]; head >= 0; head = bnode[head].next) { int eid = bnode[head].eid; merge_edge(edges[eid].from, bid, eid); next_bnode.push_back(head); } blist_head[b] = -1; clear_best_edges(b); } else merge_vertex(b, bid); } void merge_inner(int b, int bid) { if (b > N) { int v = b = base[b]; do { merge_inner(v, bid); } while ((v = node[v].next) != b); } else merge_vertex(b, bid); } void build_linked_list(int bid) { if (bid <= N) return; int last = -1; for (; merged_edge_size > 0; ) { int nid = next_bnode.back(); next_bnode.pop_back(); int by = merged_edges[--merged_edge_size], eid = best_cost[by].second; int x = edges[eid].from, y = edges[eid].to; bnode[nid] = {eid, last}; if (tree_label[y] == OUTER) update_best_edge(y, by, best_cost[by].first, eid); if (best_edge[x] == UNDEFINED || best_cost[by].first < reduced_cost(best_edge[x])) { best_edge[x] = eid; } best_cost[by] = {INF, UNDEFINED}; last = nid; } blist_head[bid] = last; } void merge_best_edges(int bid, int inner_count) { rep(i, inner_count) { int bv = outer_blossoms[outer_blossom_size + i]; if (bv >= 0) merge_outer(bv, bid), merge_inner(node[bv].next, bid); else merge_inner(~bv, bid), merge_outer(node[~bv].next, bid); } merge_outer(base[bid], bid); build_linked_list(bid); } void contract(int x, int y, int eid) { int s = surface[x], t = surface[y]; if (s == t) return; auto h = label[surface[mate[s]]].from = label[surface[mate[t]]].from = -encode(eid); int lca = -1; for (; ; label[surface[mate[s]]].from = h) { if (mate[t] != 0) swap(s, t); s = lca = surface[label[s].from]; if (label[surface[mate[s]]].from == h) break; } int inner_count = 0; for (int dir = 0; dir < 2; ++dir) { int v = (dir == 0) ? x : y; while (1) { int bv = surface[v], mv = mate[bv], bmv = surface[mv]; if (bv == lca) break; label[mv] = label[bmv] = {x, y}; auto n = node[bmv]; if (!dir) { node[bv] = {bmv, mate[mv], mv}; node[bmv].next = surface[n.to]; } else { node[surface[n.to]] = {bmv, n.to, n.from}; node[bmv] = {bv, mv, mate[mv]}; } push_outer_blossom(bmv, false); v = label[bv].from; // Caution: used as temporary array outer_blossoms[outer_blossom_size + (inner_count++)] = !dir ? bv : ~bmv; } } node[surface[y]] = {surface[x], y, x}; int bid = next_bid.back(); next_bid.pop_back(); base[bid] = lca, label[bid].from = label[lca].from, mate[bid] = mate[lca]; tree_label[bid] = OUTER; set_surface(bid, bid); merge_best_edges(bid, inner_count); outer_blossoms[outer_blossom_size++] = bid; } inline void update_best_edge(int y, int by, weight_t r_cost, int eid) { if (tree_label[by] != OUTER && best_edge[y] == UNDEFINED) { neighbors[neighbor_size++] = y; } if (best_edge[y] == UNDEFINED || r_cost < reduced_cost(best_edge[y])) { best_edge[y] = eid; } } void build_edge_list(int b) { if (b <= N) return; merge_inner(b, b); build_linked_list(b); } bool augmented(int root, int s) { if (s == 0) { int br = surface[root]; push_outer_blossom(br, true); label_blossom(br, 0, {0, 0}); build_edge_list(br); } for (; !outer_vertices.empty() || s > 0; s = 0) { auto x = (s > 0) ? s : outer_vertices.dequeue(); if (potential[x] == 0) { if (root != x) rematch(x, 0); return true; } rep(eid, offsets[x], offsets[x + 1]) { int bx = surface[x], y = edges[eid].to, by = surface[y]; if (bx == by) continue; auto r_cost = reduced_cost(x, y, edges[eid]); if (r_cost > 0 || tree_label[by] != OUTER) { update_best_edge(y, by, r_cost, eid); if (r_cost > 0) continue; } if (label[by].from >= 0) { contract(x, y, eid); continue; } if (tree_label[by] == UNUSED) { push_inner_blossom(by); if (by != y) label_blossom(by, find_mate(by), {DEFAULT, 0}); } int z = mate[by]; if (z == 0 && by != surface[root]) { rematch(x, y); rematch(y, x); return true; } int bz = surface[z]; if (label[bz].from < 0) { node[by] = {-1, y, x}; push_outer_blossom(bz, true); label_blossom(bz, mate[z], {x, y}); build_edge_list(bz); } } } return false; } void set_surface(int b, int bid) { for (int v = base[b]; surface[v] != bid; v = node[v].next) { if (v > N) tree_label[v] = UNUSED, set_surface(v, bid); surface[v] = bid; } } void reset_surface(int b, int bid) { surface[b] = bid; if (b <= N) return; for (b = base[b]; surface[b] != bid; b = node[b].next) reset_surface(b, bid); } void separate_blossom(int bid, bool push_blossom=true) { tree_label[bid] = UNUSED, label[bid].from = SEPARATED; if (bid <= N) return; if (push_blossom) inner_blossoms[inner_blossom_size++] = bid; for (int b = base[bid]; label[b].from != SEPARATED; b = node[b].next) { separate_blossom(b, false); } } void reverse_blossom(int b) { int v = b, fr = node[b].from, to = node[b].to; for (int nv = node[v].next; nv != b; ) { int nnext = node[nv].next, nfr = node[nv].from, nto = node[nv].to; node[nv].next = v, node[nv].from = to, node[nv].to = fr; fr = nfr, to = nto, v = nv, nv = nnext; } node[b].next = v, node[b].from = to, node[b].to = fr; } void expand_blossom(int bid) { next_bid.push_back(bid); tree_label[bid] = UNUSED; for (int b = base[bid]; surface[b] == bid; b = node[b].next) reset_surface(b, b); int old_base = base[bid], target = surface[node[bid].from]; if (mate[node[target].from] == node[target].to) reverse_blossom(old_base); for (int b = target; node[b].next != old_base; ) { separate_blossom(b = node[b].next); separate_blossom(b = node[b].next); } node[target] = node[bid]; for (int b = old_base; ; b = node[b].next) { label[b].from = DEFAULT, tree_label[b] = INNER; if (b > N) inner_blossoms[inner_blossom_size++] = b; int m = find_mate(b), bm = surface[m]; if (b != old_base) mate[bm] = mate[m]; label[m] = label[bm] = {node[b].to, node[b].from}; if (b == target) break; push_outer_blossom(b = node[b].next, true); build_edge_list(b); } base[bid] = bid, surface[bid] = bid; } void update_potential(int* vs, int s, weight_t delta, int label) { rep(i, s) { int x = vs[i]; if (tree_label[x] != label) continue; potential[x] += delta; } } int adjust_dual_solutions() { pair delta1(INF, 0), delta2(INF, 0), delta3(INF, 0), delta4(INF, 0); rep(i, outer_vertices.size()) { int y = outer_vertices[i], eid = best_edge[y]; delta1 = min(delta1, {potential[y], y}); if (eid != UNDEFINED) { delta3 = min(delta3, {reduced_cost(eid) >> 1, y}); } } rep(i, neighbor_size) { int y = neighbors[i]; if (tree_label[y] == UNUSED) { int eid = best_edge[y], x = edges[eid].from; delta2 = min(delta2, {reduced_cost(x, y, edges[eid]), x}); } } rep(i, inner_blossom_size) if (tree_label[inner_blossoms[i]] == INNER) { int b = inner_blossoms[i]; delta4 = min(delta4, {potential[b] >> 1, b}); } auto delta = min(min(delta1, delta2), min(delta3, delta4)); auto d = delta.first; update_potential(outer_vertices.data(), outer_vertices.size(), -1 * d, OUTER); update_potential(inner_vertices.data(), inner_vertices_size, 1 * d, INNER); update_potential(outer_blossoms.data(), outer_blossom_size, 2 * d, OUTER); update_potential(inner_blossoms.data(), inner_blossom_size, -2 * d, INNER); if (delta4.first == d) { expand_blossom(delta4.second); return -1; } else { return delta.second; } } void fix_blossom_bases() { int remain = size - next_bid.size() - (N + 1); for (int bid = N + 1; bid < size && remain > 0; ++bid) if (base[bid] != bid) { int b = base[bid]; for (int skipped = 0; skipped < 2;) { b = node[b].next; if (mate[node[b].from] == node[b].to) skipped = 0; else skipped++; } base[bid] = b; --remain; } } void free_edge_list(int x) { for (int head = blist_head[x]; head >= 0; head = bnode[head].next) { next_bnode.push_back(head); } blist_head[x] = -1; } void clear_vertices(int* vs, int size) { rep(i, size) { int v = vs[i]; label[v] = {DEFAULT, 0}; tree_label[v] = UNUSED; best_edge[v] = UNDEFINED; } } void clear_label() { label[0] = {DEFAULT, 0}; clear_vertices(outer_vertices.data(), outer_vertices.size()); outer_vertices.clear(); clear_vertices(inner_vertices.data(), inner_vertices_size); inner_vertices_size = 0; clear_vertices(outer_blossoms.data(), outer_blossom_size); rep(i, outer_blossom_size) if (blist_head[outer_blossoms[i]] >= 0) free_edge_list(outer_blossoms[i]); outer_blossom_size = 0; clear_vertices(inner_blossoms.data(), inner_blossom_size); inner_blossom_size = 0; rep(i, neighbor_size) best_edge[neighbors[i]] = UNDEFINED; neighbor_size = 0; } void set_potential() { potential.resize(size); rep(u, 1, N + 1) { weight_t max_w = 0; rep(eid, offsets[u], offsets[u + 1]) max_w = max(max_w, edges[eid].weight); potential[u] = max_w >> 1; } } void initialize() { mate.assign(size, 0); label.assign(size, {-1, 0}); surface.resize(size); rep(i, size) surface[i] = i; base.resize(size); rep(i, size) base[i] = i; node.resize(size); rep(i, size) node[i] = {i, i, i}; outer_vertices = Queue(N); inner_vertices.resize(N + 1); inner_vertices_size = 0; outer_blossoms.resize(B); outer_blossom_size = 0; inner_blossoms.resize(B); inner_blossom_size = 0; tree_label.assign(size, UNUSED); next_bid.resize(B); rep(i, B) next_bid[i] = size - 1 - i; merged_edges.resize(N + 1); merged_edge_size = 0; best_cost.assign(size, {INF, UNDEFINED}); neighbors.resize(N + 1); neighbor_size = 0; best_edge.assign(size, UNDEFINED); blist_head.assign(size, -1); next_bnode.resize(edges.size()); rep(i, edges.size()) next_bnode[i] = edges.size() - 1 - i; bnode.resize(edges.size()); } private: int N, B, size; vector edges; vector offsets; vector