最大流 & 最小割

  • Dinic 网络流 - 图1，求解二分图最大匹配时间复杂度为 网络流 - 图2 ```cpp template struct simple_queue { std::vector payload; int pos = 0; void reserve(int n) { payload.reserve(n); } int size() const { return int(payload.size()) - pos; } bool empty() const { return pos == int(payload.size()); } void push(const T& t) { payload.push_back(t); } T& front() { return payload[pos]; } void clear() {
    1.   payload.clear();
    2.   pos = 0;
    } void pop() { pos++; } };

template struct mf_graph { public: mf_graph() : _n(0) {} explicit mf_graph(int n) : _n(n), g(n) {}

  1. // 添加 from 到 to，容量为 cap 的边
  2. int add_edge(int from, int to, Cap cap) {
  3.     assert(0 <= from && from < _n);
  4.     assert(0 <= to && to < _n);
  5.     assert(0 <= cap);
  6.     int m = int(pos.size());
  7.     pos.push_back({from, int(g[from].size())});
  8.     int from_id = int(g[from].size());
  9.     int to_id = int(g[to].size());
  10.     if (from == to) to_id++;
  11.     g[from].push_back(_edge{to, to_id, cap});
  12.     g[to].push_back(_edge{from, from_id, 0});
  13.     return m;
  14. }
  16. struct edge {
  17.     int from, to;
  18.     Cap cap, flow;
  19. };
  21. // 获取第 i 次 add 的边，返回该边的 {from, to, cap, cur_cap}
  22. edge get_edge(int i) {
  23.     int m = int(pos.size());
  24.     assert(0 <= i && i < m);
  25.     auto _e = g[pos[i].first][pos[i].second];
  26.     auto _re = g[_e.to][_e.rev];
  27.     return edge{pos[i].first, _e.to, _e.cap + _re.cap, _re.cap};
  28. }
  29. // 返回所有边
  30. std::vector<edge> edges() {
  31.     int m = int(pos.size());
  32.     std::vector<edge> result;
  33.     for (int i = 0; i < m; i++) {
  34.         result.push_back(get_edge(i));
  35.     }
  36.     return result;
  37. }
  39. // 修改某条边的当前容量 new_cap - new_flow，对应反向边容量为 new_flow
  40. void change_edge(int i, Cap new_cap, Cap new_flow) {
  41.     int m = int(pos.size());
  42.     assert(0 <= i && i < m);
  43.     assert(0 <= new_flow && new_flow <= new_cap);
  44.     auto& _e = g[pos[i].first][pos[i].second];
  45.     auto& _re = g[_e.to][_e.rev];
  46.     _e.cap = new_cap - new_flow;
  47.     _re.cap = new_flow;
  48. }
  50. // 求解 s 点到 t 号点的最大流
  51. Cap flow(int s, int t) {
  52.     return flow(s, t, std::numeric_limits<Cap>::max());
  53. }
  55. // dinic，求解 s 到 t 的最大流
  56. Cap flow(int s, int t, Cap flow_limit) {
  57.     assert(0 <= s && s < _n);
  58.     assert(0 <= t && t < _n);
  59.     assert(s != t);
  61.     std::vector<int> level(_n), iter(_n);
  62.     simple_queue<int> que;
  64.     auto bfs = [&]() {
  65.         std::fill(level.begin(), level.end(), -1);
  66.         level[s] = 0;
  67.         que.clear();
  68.         que.push(s);
  69.         while (!que.empty()) {
  70.             int v = que.front();
  71.             que.pop();
  72.             for (auto e : g[v]) {
  73.                 if (e.cap == 0 || level[e.to] >= 0) continue;
  74.                 level[e.to] = level[v] + 1;
  75.                 if (e.to == t) return;
  76.                 que.push(e.to);
  77.             }
  78.         }
  79.     };
  80.     auto dfs = [&](auto self, int v, Cap up) {
  81.         if (v == s) return up;
  82.         Cap res = 0;
  83.         int level_v = level[v];
  84.         for (int& i = iter[v]; i < int(g[v].size()); i++) {
  85.             _edge& e = g[v][i];
  86.             if (level_v <= level[e.to] || g[e.to][e.rev].cap == 0) continue;
  87.             Cap d =
  88.                 self(self, e.to, std::min(up - res, g[e.to][e.rev].cap));
  89.             if (d <= 0) continue;
  90.             g[v][i].cap += d;
  91.             g[e.to][e.rev].cap -= d;
  92.             res += d;
  93.             if (res == up) return res;
  94.         }
  95.         level[v] = _n;
  96.         return res;
  97.     };
  99.     Cap flow = 0;
  100.     while (flow < flow_limit) {
  101.         bfs();
  102.         if (level[t] == -1) break;
  103.         std::fill(iter.begin(), iter.end(), 0);
  104.         Cap f = dfs(dfs, t, flow_limit - flow);
  105.         if (!f) break;
  106.         flow += f;
  107.     }
  108.     return flow;
  109. }
  111. //返回最小割，需要提前调用 flow
  112. std::vector<bool> min_cut(int s) {
  113.     std::vector<bool> visited(_n);
  114.     simple_queue<int> que;
  115.     que.push(s);
  116.     while (!que.empty()) {
  117.         int p = que.front();
  118.         que.pop();
  119.         visited[p] = true;
  120.         for (auto e : g[p]) {
  121.             if (e.cap && !visited[e.to]) {
  122.                 visited[e.to] = true;
  123.                 que.push(e.to);
  124.             }
  125.         }
  126.     }
  127.     return visited;
  128. }

private: int _n; struct _edge { int to, rev; Cap cap; }; std::vector> pos; std::vector> g; };

  1. <a name="jgyE2"></a>
  2. ## 最小费用流
  3. ![](https://cdn.nlark.com/yuque/__latex/98921346cebd0385c2c3194e13a7c3b3.svg#card=math&code=Primal-Dual&id=AEYYn)最小费用最大流，复杂度![](https://cdn.nlark.com/yuque/__latex/1185efd542e780997fe7c2848acd595a.svg#card=math&code=O%28F%20%28n%20%2B%20m%29%20%5Clog%20%28n%20%2B%20m%29%29&id=uJPKE)。<br />以下模版中，边费用必须为非零整数。所以只能用来求解最小费用。
  4. ```cpp
  5. template <class E> struct csr {
  6.     std::vector<int> start;
  7.     std::vector<E> elist;
  8.     explicit csr(int n, const std::vector<std::pair<int, E>>& edges)
  9.         : start(n + 1), elist(edges.size()) {
  10.         for (auto e : edges) {
  11.             start[e.first + 1]++;
  12.         }
  13.         for (int i = 1; i <= n; i++) {
  14.             start[i] += start[i - 1];
  15.         }
  16.         auto counter = start;
  17.         for (auto e : edges) {
  18.             elist[counter[e.first]++] = e.second;
  19.         }
  20.     }
  21. };
  23. template <class Cap, class Cost> 
  24. struct mcf_graph {
  25. public:
  26.     mcf_graph() {}
  27.     explicit mcf_graph(int n) : _n(n) {}
  28.     // 添加一条从 from 到 to，容量为 cap，单位费用为 cost 的边
  29.     int add_edge(int from, int to, Cap cap, Cost cost) {
  30.         assert(0 <= from && from < _n);
  31.         assert(0 <= to && to < _n);
  32.         assert(0 <= cap);
  33.         assert(0 <= cost);
  34.         int m = int(_edges.size());
  35.         _edges.push_back({from, to, cap, 0, cost});
  36.         return m;
  37.     }
  39.     struct edge {
  40.         int from, to;
  41.         Cap cap, flow;
  42.         Cost cost;
  43.     };
  44.     // 获取第 i 条边（按照 add_edge 的顺序）
  45.     edge get_edge(int i) {
  46.         int m = int(_edges.size());
  47.         assert(0 <= i && i < m);
  48.         return _edges[i];
  49.     }
  50.     // 获取所有边
  51.     std::vector<edge> edges() { return _edges; }
  53.     std::pair<Cap, Cost> flow(int s, int t) {
  54.         return flow(s, t, std::numeric_limits<Cap>::max());
  55.     }
  56.     std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
  57.         return slope(s, t, flow_limit).back();
  58.     }
  59.     std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
  60.         return slope(s, t, std::numeric_limits<Cap>::max());
  61.     }
  63.     std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
  64.         assert(0 <= s && s < _n);
  65.         assert(0 <= t && t < _n);
  66.         assert(s != t);
  68.         int m = int(_edges.size());
  69.         std::vector<int> edge_idx(m);
  71.         auto g = [&]() {
  72.             std::vector<int> degree(_n), redge_idx(m);
  73.             std::vector<std::pair<int, _edge>> elist;
  74.             elist.reserve(2 * m);
  75.             for (int i = 0; i < m; i++) {
  76.                 auto e = _edges[i];
  77.                 edge_idx[i] = degree[e.from]++;
  78.                 redge_idx[i] = degree[e.to]++;
  79.                 elist.push_back({e.from, {e.to, -1, e.cap - e.flow, e.cost}});
  80.                 elist.push_back({e.to, {e.from, -1, e.flow, -e.cost}});
  81.             }
  82.             auto _g = csr<_edge>(_n, elist);
  83.             for (int i = 0; i < m; i++) {
  84.                 auto e = _edges[i];
  85.                 edge_idx[i] += _g.start[e.from];
  86.                 redge_idx[i] += _g.start[e.to];
  87.                 _g.elist[edge_idx[i]].rev = redge_idx[i];
  88.                 _g.elist[redge_idx[i]].rev = edge_idx[i];
  89.             }
  90.             return _g;
  91.         }();
  93.         auto result = slope(g, s, t, flow_limit);
  95.         for (int i = 0; i < m; i++) {
  96.             auto e = g.elist[edge_idx[i]];
  97.             _edges[i].flow = _edges[i].cap - e.cap;
  98.         }
  100.         return result;
  101.     }
  103. private:
  104.     int _n;
  105.     std::vector<edge> _edges;
  107.     // inside edge
  108.     struct _edge {
  109.         int to, rev;
  110.         Cap cap;
  111.         Cost cost;
  112.     };
  114.     std::vector<std::pair<Cap, Cost>> slope(csr<_edge>& g, int s, int t, Cap flow_limit) { 
  115.         // variants (C = maxcost):
  116.         // -(n-1)C <= dual[s] <= dual[i] <= dual[t] = 0
  117.         // reduced cost (= e.cost + dual[e.from] - dual[e.to]) >= 0 for all edge
  119.         // dual_dist[i] = (dual[i], dist[i])
  120.         std::vector<std::pair<Cost, Cost>> dual_dist(_n);
  121.         std::vector<int> prev_e(_n);
  122.         std::vector<bool> vis(_n);
  123.         struct Q {
  124.             Cost key;
  125.             int to;
  126.             bool operator<(Q r) const { return key > r.key; }
  127.         };
  128.         std::vector<int> que_min;
  129.         std::vector<Q> que;
  130.         auto dual_ref = [&]() {
  131.             for (int i = 0; i < _n; i++) {
  132.                 dual_dist[i].second = std::numeric_limits<Cost>::max();
  133.             }
  134.             std::fill(vis.begin(), vis.end(), false);
  135.             que_min.clear();
  136.             que.clear();
  138.             // que[0..heap_r) was heapified
  139.             size_t heap_r = 0;
  141.             dual_dist[s].second = 0;
  142.             que_min.push_back(s);
  143.             while (!que_min.empty() || !que.empty()) {
  144.                 int v;
  145.                 if (!que_min.empty()) {
  146.                     v = que_min.back();
  147.                     que_min.pop_back();
  148.                 } else {
  149.                     while (heap_r < que.size()) {
  150.                         heap_r++;
  151.                         std::push_heap(que.begin(), que.begin() + heap_r);
  152.                     }
  153.                     v = que.front().to;
  154.                     std::pop_heap(que.begin(), que.end());
  155.                     que.pop_back();
  156.                     heap_r--;
  157.                 }
  158.                 if (vis[v]) continue;
  159.                 vis[v] = true;
  160.                 if (v == t) break;
  161.                 // dist[v] = shortest(s, v) + dual[s] - dual[v]
  162.                 // dist[v] >= 0 (all reduced cost are positive)
  163.                 // dist[v] <= (n-1)C
  164.                 Cost dual_v = dual_dist[v].first, dist_v = dual_dist[v].second;
  165.                 for (int i = g.start[v]; i < g.start[v + 1]; i++) {
  166.                     auto e = g.elist[i];
  167.                     if (!e.cap) continue;
  168.                     // |-dual[e.to] + dual[v]| <= (n-1)C
  169.                     // cost <= C - -(n-1)C + 0 = nC
  170.                     Cost cost = e.cost - dual_dist[e.to].first + dual_v;
  171.                     if (dual_dist[e.to].second - dist_v > cost) {
  172.                         Cost dist_to = dist_v + cost;
  173.                         dual_dist[e.to].second = dist_to;
  174.                         prev_e[e.to] = e.rev;
  175.                         if (dist_to == dist_v) {
  176.                             que_min.push_back(e.to);
  177.                         } else {
  178.                             que.push_back(Q{dist_to, e.to});
  179.                         }
  180.                     }
  181.                 }
  182.             }
  183.             if (!vis[t]) {
  184.                 return false;
  185.             }
  187.             for (int v = 0; v < _n; v++) {
  188.                 if (!vis[v]) continue;
  189.                 // dual[v] = dual[v] - dist[t] + dist[v]
  190.                 //         = dual[v] - (shortest(s, t) + dual[s] - dual[t]) +
  191.                 //         (shortest(s, v) + dual[s] - dual[v]) = - shortest(s,
  192.                 //         t) + dual[t] + shortest(s, v) = shortest(s, v) -
  193.                 //         shortest(s, t) >= 0 - (n-1)C
  194.                 dual_dist[v].first -= dual_dist[t].second - dual_dist[v].second;
  195.             }
  196.             return true;
  197.         };
  198.         Cap flow = 0;
  199.         Cost cost = 0, prev_cost_per_flow = -1;
  200.         std::vector<std::pair<Cap, Cost>> result = {{Cap(0), Cost(0)}};
  201.         while (flow < flow_limit) {
  202.             if (!dual_ref()) break;
  203.             Cap c = flow_limit - flow;
  204.             for (int v = t; v != s; v = g.elist[prev_e[v]].to) {
  205.                 c = std::min(c, g.elist[g.elist[prev_e[v]].rev].cap);
  206.             }
  207.             for (int v = t; v != s; v = g.elist[prev_e[v]].to) {
  208.                 auto& e = g.elist[prev_e[v]];
  209.                 e.cap += c;
  210.                 g.elist[e.rev].cap -= c;
  211.             }
  212.             Cost d = -dual_dist[s].first;
  213.             flow += c;
  214.             cost += c * d;
  215.             if (prev_cost_per_flow == d) {
  216.                 result.pop_back();
  217.             }
  218.             result.push_back({flow, cost});
  219.             prev_cost_per_flow = d;
  220.         }
  221.         return result;
  222.     }
  223. };

最大费用流

  1. template <class T> struct simple_queue {
  2.     std::vector<T> payload;
  3.     int pos = 0;
  4.     void reserve(int n) { payload.reserve(n); }
  5.     int size() const { return int(payload.size()) - pos; }
  6.     bool empty() const { return pos == int(payload.size()); }
  7.     void push(const T& t) { payload.push_back(t); }
  8.     T& front() { return payload[pos]; }
  9.     void clear() {
  10.         payload.clear();
  11.         pos = 0;
  12.     }
  13.     void pop() { pos++; }
  14. };
  17. template <class E> struct csr {
  18.     std::vector<int> start;
  19.     std::vector<E> elist;
  20.     explicit csr(int n, const std::vector<std::pair<int, E>>& edges)
  21.         : start(n + 1), elist(edges.size()) {
  22.         for (auto e : edges) {
  23.             start[e.first + 1]++;
  24.         }
  25.         for (int i = 1; i <= n; i++) {
  26.             start[i] += start[i - 1];
  27.         }
  28.         auto counter = start;
  29.         for (auto e : edges) {
  30.             elist[counter[e.first]++] = e.second;
  31.         }
  32.     }
  33. };
  35. template <class Cap, class Cost> 
  36. struct mcf_graph {
  37. public:
  38.     mcf_graph() {}
  39.     explicit mcf_graph(int n) : _n(n) {}
  40.     // 添加一条从 from 到 to，容量为 cap，单位费用为 cost 的边
  41.     int add_edge(int from, int to, Cap cap, Cost cost) {
  42.         assert(0 <= from && from < _n);
  43.         assert(0 <= to && to < _n);
  44.         assert(0 <= cap);
  45.         // assert(0 <= cost);
  46.         int m = int(_edges.size());
  47.         _edges.push_back({from, to, cap, 0, cost});
  48.         return m;
  49.     }
  51.     struct edge {
  52.         int from, to;
  53.         Cap cap, flow;
  54.         Cost cost;
  55.     };
  56.     // 获取第 i 条边（按照 add_edge 的顺序）
  57.     edge get_edge(int i) {
  58.         int m = int(_edges.size());
  59.         assert(0 <= i && i < m);
  60.         return _edges[i];
  61.     }
  62.     // 获取所有边
  63.     std::vector<edge> edges() { return _edges; }
  65.     std::pair<Cap, Cost> flow(int s, int t) {
  66.         return flow(s, t, std::numeric_limits<Cap>::max());
  67.     }
  68.     std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
  69.         return slope(s, t, flow_limit).back();
  70.     }
  71.     std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
  72.         return slope(s, t, std::numeric_limits<Cap>::max());
  73.     }
  75.     std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
  76.         assert(0 <= s && s < _n);
  77.         assert(0 <= t && t < _n);
  78.         assert(s != t);
  80.         int m = int(_edges.size());
  81.         std::vector<int> edge_idx(m);
  83.         auto g = [&]() {
  84.             std::vector<int> degree(_n), redge_idx(m);
  85.             std::vector<std::pair<int, _edge>> elist;
  86.             elist.reserve(2 * m);
  87.             for (int i = 0; i < m; i++) {
  88.                 auto e = _edges[i];
  89.                 edge_idx[i] = degree[e.from]++;
  90.                 redge_idx[i] = degree[e.to]++;
  91.                 elist.push_back({e.from, {e.to, -1, e.cap - e.flow, e.cost}});
  92.                 elist.push_back({e.to, {e.from, -1, e.flow, -e.cost}});
  93.             }
  94.             auto _g = csr<_edge>(_n, elist);
  95.             for (int i = 0; i < m; i++) {
  96.                 auto e = _edges[i];
  97.                 edge_idx[i] += _g.start[e.from];
  98.                 redge_idx[i] += _g.start[e.to];
  99.                 _g.elist[edge_idx[i]].rev = redge_idx[i];
  100.                 _g.elist[redge_idx[i]].rev = edge_idx[i];
  101.             }
  102.             return _g;
  103.         }();
  105.         auto result = slope(g, s, t, flow_limit);
  107.         for (int i = 0; i < m; i++) {
  108.             auto e = g.elist[edge_idx[i]];
  109.             _edges[i].flow = _edges[i].cap - e.cap;
  110.         }
  112.         return result;
  113.     }
  115. private:
  116.     int _n;
  117.     std::vector<edge> _edges;
  119.     // inside edge
  120.     struct _edge {
  121.         int to, rev;
  122.         Cap cap;
  123.         Cost cost;
  124.     };
  126.     std::vector<Cost> spfa(csr<_edge>& g, int s, int t) {
  127.         std::vector<Cost> dist(_n, std::numeric_limits<Cost>::max() / 2);
  128.         std::vector<bool> vis(_n);
  129.         simple_queue<int> q;
  130.         q.push(s); dist[s] = 0;
  131.         while(q.size()) {
  132.             int x = q.front(); q.pop();
  133.             vis[x] = 0;
  134.             for(int i = g.start[x]; i < g.start[x + 1]; i ++) {
  135.                 auto e = g.elist[i];
  136.                 int y = e.to;
  137.                 if(e.cap > 0 && dist[y] - dist[x] > e.cost) {
  138.                     dist[y] = dist[x] + e.cost;
  139.                     if(vis[y]) continue;
  140.                     vis[y] = 1; q.push(y);
  141.                 }
  142.             }
  143.         }
  144.         return dist;
  145.     }
  147.     std::vector<std::pair<Cap, Cost>> slope(csr<_edge>& g, int s, int t, Cap flow_limit) { 
  148.         std::vector<std::pair<Cost, Cost>> dual_dist(_n);
  149.         std::vector<int> prev_e(_n);
  150.         std::vector<bool> vis(_n);
  151.         struct Q {
  152.             Cost key;
  153.             int to;
  154.             bool operator<(Q r) const { return key > r.key; }
  155.         };
  156.         std::vector<int> que_min;
  157.         std::vector<Q> que;
  158.         auto h = spfa(g, s, t);
  159.         for(int i = 0; i < _n; i++) dual_dist[i].first = h[i];
  161.         auto dual_ref = [&]() {
  162.             for (int i = 0; i < _n; i++) {
  163.                 dual_dist[i].second = std::numeric_limits<Cost>::max();
  164.             }
  165.             std::fill(vis.begin(), vis.end(), false);
  166.             que_min.clear();
  167.             que.clear();
  169.             // que[0..heap_r) was heapified
  170.             size_t heap_r = 0;
  172.             dual_dist[s].second = 0;
  173.             que_min.push_back(s);
  174.             while (!que_min.empty() || !que.empty()) {
  175.                 int v;
  176.                 if (!que_min.empty()) {
  177.                     v = que_min.back();
  178.                     que_min.pop_back();
  179.                 } else {
  180.                     while (heap_r < que.size()) {
  181.                         heap_r++;
  182.                         std::push_heap(que.begin(), que.begin() + heap_r);
  183.                     }
  184.                     v = que.front().to;
  185.                     std::pop_heap(que.begin(), que.end());
  186.                     que.pop_back();
  187.                     heap_r--;
  188.                 }
  189.                 if (vis[v]) continue;
  190.                 vis[v] = true;
  191.                 Cost dual_v = dual_dist[v].first, dist_v = dual_dist[v].second;
  192.                 for (int i = g.start[v]; i < g.start[v + 1]; i++) {
  193.                     auto e = g.elist[i];
  194.                     if (!e.cap) continue;
  195.                     // |-dual[e.to] + dual[v]| <= (n-1)C
  196.                     // cost <= C - -(n-1)C + 0 = nC
  197.                     Cost cost = e.cost - dual_dist[e.to].first + dual_v;
  198.                     if (dual_dist[e.to].second - dist_v > cost) {
  199.                         Cost dist_to = dist_v + cost;
  200.                         assert(dist_to >= 0);
  201.                         dual_dist[e.to].second = dist_to;
  202.                         prev_e[e.to] = e.rev;
  203.                         if (dist_to == dist_v) {
  204.                             que_min.push_back(e.to);
  205.                         } else {
  206.                             que.push_back(Q{dist_to, e.to});
  207.                         }
  208.                     }
  209.                 }
  210.             }
  211.             if (!vis[t]) {
  212.                 return false;
  213.             }
  215.             for (int v = 0; v < _n; v++) {
  216.                 dual_dist[v].first += dual_dist[v].second;
  217.             }
  218.             return true;
  219.         };
  220.         Cap flow = 0;
  221.         Cost cost = 0, prev_cost_per_flow = -1;
  222.         std::vector<std::pair<Cap, Cost>> result = {{Cap(0), Cost(0)}};
  223.         while (flow < flow_limit) {
  224.             if (!dual_ref()) break;
  225.             Cap c = flow_limit - flow;
  226.             for (int v = t; v != s; v = g.elist[prev_e[v]].to) {
  227.                 c = std::min(c, g.elist[g.elist[prev_e[v]].rev].cap);
  228.             }
  229.             for (int v = t; v != s; v = g.elist[prev_e[v]].to) {
  230.                 auto& e = g.elist[prev_e[v]];
  231.                 e.cap += c;
  232.                 g.elist[e.rev].cap -= c;
  233.             }
  234.             Cost d = dual_dist[t].first;
  235.             flow += c;
  236.             cost += c * d;
  237.             if (prev_cost_per_flow == d) {
  238.                 result.pop_back();
  239.             }
  240.             result.push_back({flow, cost});
  241.             prev_cost_per_flow = d;
  242.         }
  243.         return result;
  244.     }
  245. };