// https://yukicoder.me/submissions/592136 改変 山登り法 // #pragma GCC optimize("Ofast") // #pragma GCC optimize("unroll-loops") #include #include #include #include #include #include #include #include #include ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // MinCostFlow by https://github.com/yosupo06/library-checker-problems/blob/master/graph/min_cost_b_flow/sol/correct.cpp // Copyright 2020 yosupo06/library-checker-problems ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// enum Objective { MINIMIZE = 1, MAXIMIZE = -1, }; enum class Status { OPTIMAL, INFEASIBLE, }; template class MinCostFlow { using V_id = uint32_t; using E_id = uint32_t; class Edge { friend class MinCostFlow; V_id src, dst; Flow flow, cap; Cost cost; E_id rev; public: Edge() = default; Edge(const V_id src, const V_id dst, const Flow cap, const Cost cost, const E_id rev) : src(src), dst(dst), flow(0), cap(cap), cost(cost), rev(rev) {} [[nodiscard]] Flow residual_cap() const { return cap - flow; } }; public: class EdgePtr { friend class MinCostFlow; const MinCostFlow *instance; V_id v; E_id e; EdgePtr(const MinCostFlow * const instance, const V_id v, const E_id e) : instance(instance), v(v), e(e) {} [[nodiscard]] const Edge &edge() const { return instance->g[v][e]; } [[nodiscard]] const Edge &rev() const { const Edge &e = edge(); return instance->g[e.dst][e.rev]; } public: EdgePtr() = default; [[nodiscard]] V_id src() const { return v; } [[nodiscard]] V_id dst() const { return edge().dst; } [[nodiscard]] Flow flow() const { return edge().flow; } [[nodiscard]] Flow lower() const { return -rev().cap; } [[nodiscard]] Flow upper() const { return edge().cap; } [[nodiscard]] Cost cost() const { return edge().cost; } [[nodiscard]] Cost gain() const { return -edge().cost; } }; private: V_id n; std::vector> g; std::vector b; public: MinCostFlow() : n(0) {} V_id add_vertex() { ++n; g.resize(n); b.resize(n); return n-1; } std::vector add_vertices(const size_t size) { std::vector ret(size); std::iota(std::begin(ret), std::end(ret), n); n += size; g.resize(n); b.resize(n); return ret; } EdgePtr add_edge(const V_id src, const V_id dst, const Flow lower, const Flow upper, const Cost cost) { const E_id e = g[src].size(), re = src == dst ? e + 1 : g[dst].size(); assert(lower <= upper); g[src].emplace_back(Edge{src, dst, upper, cost * objective, re}); g[dst].emplace_back(Edge{dst, src, -lower, -cost * objective, e}); return EdgePtr{this, src, e}; } void add_supply(const V_id v, const Flow amount) { b[v] += amount; } void add_demand(const V_id v, const Flow amount) { b[v] -= amount; } private: // Variables used in calculation const Cost unreachable = std::numeric_limits::max(); Cost farthest; std::vector potential; std::vector dist; std::vector parent; // out-forrest. std::priority_queue, std::vector>, std::greater<>> pq; // should be empty outside of dual() std::vector excess_vs, deficit_vs; Edge &rev(const Edge &e) { return g[e.dst][e.rev]; } void push(Edge &e, const Flow amount) { e.flow += amount; g[e.dst][e.rev].flow -= amount; } Cost residual_cost(const V_id src, const V_id dst, const Edge &e) { return e.cost + potential[src] - potential[dst]; } bool dual(const Flow delta) { dist.assign(n, unreachable); parent.assign(n, nullptr); excess_vs.erase(std::remove_if(std::begin(excess_vs), std::end(excess_vs), [&](const V_id v) { return b[v] < delta; }), std::end(excess_vs)); deficit_vs.erase(std::remove_if(std::begin(deficit_vs), std::end(deficit_vs), [&](const V_id v) { return b[v] > -delta; }), std::end(deficit_vs)); for (const auto v : excess_vs) pq.emplace(dist[v] = 0, v); farthest = 0; std::size_t deficit_count = 0; while (!pq.empty()) { Cost d; std::size_t u; std::tie(d, u) = pq.top(); // const auto [d, u] = pq.top(); pq.pop(); if (dist[u] < d) continue; farthest = d; if (b[u] <= -delta) ++deficit_count; if (deficit_count >= deficit_vs.size()) break; for (auto &e : g[u]) { if (e.residual_cap() < delta) continue; const auto v = e.dst; const auto new_dist = d + residual_cost(u, v, e); if (new_dist >= dist[v]) continue; pq.emplace(dist[v] = new_dist, v); parent[v] = &e; } } pq = decltype(pq)(); // pq.clear() doesn't exist. for (V_id v = 0; v < n; ++v) { potential[v] += std::min(dist[v], farthest); } return deficit_count > 0; } void primal(const Flow delta) { for (const auto t : deficit_vs) { if (dist[t] > farthest) continue; Flow f = -b[t]; V_id v; for (v = t; parent[v] != nullptr; v = parent[v]->src) { f = std::min(f, parent[v]->residual_cap()); } f = std::min(f, b[v]); f -= f % delta; if (f <= 0) continue; for (v = t; parent[v] != nullptr;) { auto &e = *parent[v]; push(e, f); int u = parent[v]->src; if (e.residual_cap() <= 0) parent[v] = nullptr; v = u; } b[t] += f; b[v] -= f; } } void saturate_negative(const Flow delta) { excess_vs.clear(); deficit_vs.clear(); for (auto &es : g) for (auto &e : es) { Flow rcap = e.residual_cap(); rcap -= rcap % delta; const Cost rcost = residual_cost(e.src, e.dst, e); if (rcost < 0 || rcap < 0) { push(e, rcap); b[e.src] -= rcap; b[e.dst] += rcap; } } for (V_id v = 0; v < n; ++v) if (b[v] != 0) { (b[v] > 0 ? excess_vs : deficit_vs).emplace_back(v); } } public: std::pair solve() { potential.resize(n); Flow inf_flow = 1; for (const auto t : b) inf_flow = std::max({inf_flow, t, -t}); for (const auto &es : g) for (const auto &e : es) inf_flow = std::max({inf_flow, e.residual_cap(), -e.residual_cap()}); Flow delta = 1; while (delta < inf_flow) delta *= SCALING_FACTOR; for (; delta; delta /= SCALING_FACTOR) { saturate_negative(delta); while (dual(delta)) primal(delta); } Cost value = 0; for (const auto &es : g) for (const auto &e : es) { value += e.flow * e.cost; } value /= 2; if (excess_vs.empty() && deficit_vs.empty()) { return { Status::OPTIMAL, value / objective }; } else { return { Status::INFEASIBLE, value / objective }; } } std::vector get_potential() { // Not strictly necessary, but re-calculate potential to bound the potential values, // plus make them somewhat canonical so that it is robust for the algorithm chaneges. std::fill(std::begin(potential), std::end(potential), 0); for (size_t i = 0; i < n; ++i) for (const auto &es : g) for (const auto &e : es) if (e.residual_cap() > 0) potential[e.dst] = std::min(potential[e.dst], potential[e.src] + e.cost); return potential; } template T get_result_value() { T value = 0; for (const auto &es : g) for (const auto &e : es) { value += (T)(e.flow) * (T)(e.cost); } value /= (T)2; return value; } std::vector get_cut() { std::vector res; if (excess_vs.empty()) return res; for (size_t v = 0; v < n; ++v) { if (deficit_vs.empty() || (dist[v] < unreachable)) res.emplace_back(v); } return res; } }; ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // MinCostFlow END ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// using namespace std; constexpr long long INFll = (1ll << 60) - 1; int n, m, k; vector> E; vector depth, parent, edge; vector depth_, parent_, edge_; void build(const vector &v) { depth.assign(n, 0); parent.assign(n, 0); edge.assign(n, 0); // depth[0] = 0; vector>> g(n); for(int i : v){ int a = E[i][0], b = E[i][1]; g[a].emplace_back(b, i); g[b].emplace_back(a, i); } auto dfs = [&](int from, int at, auto dfs) -> void { const int d2 = depth[at] + 1; for(const auto& [to, id] : g[at]) if(to != from){ depth[to] = d2; parent[to] = at; edge[to] = id; dfs(at, to, dfs); } }; dfs(-1, 0, dfs); } vector exchangable_edges(int e){ vector ans; int a = E[e][0], b = E[e][1]; while(a != b){ if(depth[a] < depth[b]) swap(a, b); ans.push_back(edge[a]); a = parent[a]; } return ans; } int64_t flow_us; // 「G の最小全域木として vl が採用されうる」という条件の下で x を動かしたときの最大値を求める. long long subsolve(const vector& vl, const vector& vr){ long long res = 0; for (auto i : vl) res += E[i][2]; res *= k; build(vl); MinCostFlow mcf; const auto vs = mcf.add_vertices(m + 2); int s = m, t = m + 1; const long long BIG = 1e10; for (auto j : vr) { // vl に含まれない辺 j を追加したとき,代わりに取り除ける辺 i を列挙 for (auto i : exchangable_edges(j)) mcf.add_edge(i, j, 0, BIG, E[j][2] - E[i][2]); } for (auto i : vl) mcf.add_edge(s, i, max(k - E[i][3], 0), BIG, 0); for (auto j : vr) mcf.add_edge(j, t, 0, E[j][3], 0); mcf.add_edge(t, s, 0, BIG, 0); auto START = std::chrono::system_clock::now(); auto [status, f] = mcf.solve(); flow_us += std::chrono::duration_cast(std::chrono::system_clock::now() - START).count(); if (status == Status::INFEASIBLE) return INFll; res += f; return res; } struct UnionFind{ vector data; UnionFind(int n): data(n, -1){} bool unite(int a, int b){ a = root(a); b = root(b); if(a == b) return 0; if(data[a] > data[b]) swap(a, b); data[a] += data[b]; data[b] = a; return 1; } bool find(int a, int b){ return root(a) == root(b); } int root(int a){ return data[a] < 0 ? a : data[a] = root(data[a]); } int size(int a){ return -data[root(a)]; } int operator[](int a){ return root(a); } }; struct Xorshift64{ using result_type = uint32_t; static constexpr result_type min(){ return 0; } static constexpr result_type max(){ return -1; } uint64_t x = 1; result_type operator()(){ x ^= (x << 13); x ^= (x >> 7); x ^= (x << 17); return static_cast(x); } }rnd; int main(){ cin.tie(nullptr); ios::sync_with_stdio(false); auto start = chrono::system_clock::now(); cin >> n >> m >> k; E.resize(m); for (auto& [a, b, c, d] : E){ cin >> a >> b >> c >> d; a--; b--; } sort(E.begin(), E.end(), [](const auto& a, const auto& b){ return a[2] < b[2]; }); UnionFind uf(n); vector vl, vr; for (int i = 0; i < m; i++) { int a = E[i][0], b = E[i][1]; if(uf.unite(a, b)) vl.push_back(i); else vr.push_back(i); } long long ans = subsolve(vl, vr); int nbupd = 0, nb_loop = 0; if(vr.size()){ while (chrono::duration_cast(chrono::system_clock::now() - start).count() < 1900) { nb_loop++; auto p = vr.begin() + rnd() % vr.size(); auto v = exchangable_edges(*p); auto q = find(vl.begin(), vl.end(), v[rnd() % v.size()]); iter_swap(p, q); swap(depth, depth_); swap(parent, parent_); swap(edge, edge_); const long long ans2 = subsolve(vl, vr); if (ans < ans2) nbupd++; ans = max(ans, ans2); if(ans != ans2){ iter_swap(p, q); depth = move(depth_); parent = move(parent_); edge = move(edge_); } } } cerr << nb_loop << ' ' << nbupd << ' ' << flow_us / 1000 << '\n'; if (ans == INFll) ans = -1; cout << ans << '\n'; }