fork(1) download
copy 
  1. #include<bits/stdc++.h>
  2. using namespace std;
  3.  
  4. //Works for both directed, undirected and with negative cost too
  5. //doesn't work for negative cycles
  6. //for undirected edges just make the directed flag false
  7. //Complexity: O(min(E^2 *V log V, E logV * flow))
  8. // -> no, it's Ω (2^{V/2}) in the worst case
  9. using T = long long;
  10. const T inf = 1LL << 61;
  11. struct MCMF {
  12.   struct edge {
  13.     int u, v;
  14.     T cap, cost;
  15.     int id;
  16.     edge(int _u, int _v, T _cap, T _cost, int _id) {
  17.       u = _u;
  18.       v = _v;
  19.       cap = _cap;
  20.       cost = _cost;
  21.       id = _id;
  22.     }
  23.   };
  24.   int n, s, t, mxid;
  25.   T flow, cost;
  26.   vector<vector<int>> g;
  27.   vector<edge> e;
  28.   vector<T> d, potential, flow_through;
  29.   vector<int> par;
  30.   uint64_t steps = 0;
  31.   bool neg;
  32.   MCMF() {}
  33.   MCMF(int _n) { // 0-based indexing
  34.     n = _n + 10;
  35.     g.assign(n, vector<int> ());
  36.     neg = false;
  37.     mxid = 0;
  38.   }
  39.   void add_edge(int u, int v, T cap, T cost, int id = -1, bool directed = true) {
  40.     if(cost < 0) neg = true;
  41.     g[u].push_back(e.size());
  42.     e.push_back(edge(u, v, cap, cost, id));
  43.     g[v].push_back(e.size());
  44.     e.push_back(edge(v, u, 0, -cost, -1));
  45.     mxid = max(mxid, id);
  46.     if(!directed) add_edge(v, u, cap, cost, -1, true);
  47.   }
  48.   bool dijkstra() {
  49.       ++steps;
  50.     par.assign(n, -1);
  51.     d.assign(n, inf);
  52.     priority_queue<pair<T, T>, vector<pair<T, T>>, greater<pair<T, T>> > q;
  53.     d[s] = 0;
  54.     q.push(pair<T, T>(0, s));
  55.     while (!q.empty()) {
  56.       int u = q.top().second;
  57.       T nw = q.top().first;
  58.       q.pop();
  59.       if(nw != d[u]) continue;
  60.       for (int i = 0; i < (int)g[u].size(); i++) {
  61.         int id = g[u][i];
  62.         int v = e[id].v;
  63.         T cap = e[id].cap;
  64.         T w = e[id].cost + potential[u] - potential[v];
  65.         if (d[u] + w < d[v] && cap > 0) {
  66.           d[v] = d[u] + w;
  67.           par[v] = id;
  68.           q.push(pair<T, T>(d[v], v));
  69.         }
  70.       }
  71.     }
  72.     for (int i = 0; i < n; i++) { // update potential
  73.       if(d[i] < inf) potential[i] += d[i];
  74.     }
  75.     return d[t] != inf;
  76.   }
  77.   T send_flow(int v, T cur) {
  78.     if(par[v] == -1) return cur;
  79.     int id = par[v];
  80.     int u = e[id].u;
  81.     T w = e[id].cost;
  82.     T f = send_flow(u, min(cur, e[id].cap));
  83.     cost += f * w;
  84.     e[id].cap -= f;
  85.     e[id ^ 1].cap += f;
  86.     return f;
  87.   }
  88.   //returns {maxflow, mincost}
  89.   pair<T, T> solve(int _s, int _t, T goal = inf) {
  90.     s = _s;
  91.     t = _t;
  92.     flow = 0, cost = 0;
  93.     potential.assign(n, 0);
  94.     if (neg) {
  95.       // run Bellman-Ford to find starting potential
  96.       d.assign(n, inf);
  97.       for (int i = 0, relax = true; i < n && relax; i++) {
  98.         for (int u = 0; u < n; u++) {
  99.           for (int k = 0; k < (int)g[u].size(); k++) {
  100.             int id = g[u][k];
  101.             int v = e[id].v;
  102.             T cap = e[id].cap, w = e[id].cost;
  103.             if (d[v] > d[u] + w && cap > 0) {
  104.               d[v] = d[u] + w;
  105.               relax = true;
  106.             }
  107.           }
  108.         }
  109.       }
  110.       for(int i = 0; i < n; i++) if(d[i] < inf) potential[i] = d[i];
  111.     }
  112.     while (flow < goal && dijkstra()) flow += send_flow(t, goal - flow);
  113.     flow_through.assign(mxid + 10, 0);
  114.     for (int u = 0; u < n; u++) {
  115.       for (auto v : g[u]) {
  116.         if (e[v].id >= 0) flow_through[e[v].id] = e[v ^ 1].cap;
  117.       }
  118.     }
  119.     return make_pair(flow, cost);
  120.   }
  121. };
  122. void run(int k){
  123. 	// A bad network problem for the simplex method and other minimum cost flow algorithms, 1973, Norman Zadeh  
  124.     const int inf = 5 * (2<<(k-2));
  125.     const int n = 2*k+2, m = k*(k+1)+1;
  126.  
  127.     MCMF fl(n+1);
  128.     for(int i=0; i<k; ++i){
  129.         fl.add_edge(1, i+2, i==0 ? 1 : i==1 ? 3 : 5<<(i-2), 0);
  130.     }
  131.     for(int i=0; i<k; ++i){
  132.         for(int j=0; j<k; ++j) if(j != i){
  133.             fl.add_edge(i+2, k+j+2, inf, (1 << max(i, j))-1);
  134.         }
  135.     }
  136.     for(int i=0; i<k; ++i){
  137.         fl.add_edge(i+2+k, 2*k+2, i<2 ? 2 : 5<<(i-2), 0);
  138.     }
  139.     fl.add_edge(2, k+2, inf, 0);
  140.     auto t1 = chrono::high_resolution_clock::now();
  141.     auto ans = fl.solve(1, n);
  142.     assert(ans.first != ans.second); // avoid optimization
  143.     auto t2 = chrono::high_resolution_clock::now();
  144.     cout << n << " " << m << " -> " << chrono::duration_cast<chrono::nanoseconds>(t2-t1).count()*1e-9 << "\n";
  145.     cerr << "steps: " << fl.steps << "\n";
  146. }
  147. int main() {
  148.   for(int k=13; k<19; ++k){
  149.       run(k);
  150.   }
  151. }
  152.  
Success #stdin #stdout #stderr 4.27s 5600KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
28 183 -> 0.0425663
30 211 -> 0.095778
32 241 -> 0.21483
34 273 -> 0.489516
36 307 -> 1.07289
38 343 -> 2.36075
stderr
copy 
steps: 10239
steps: 20479
steps: 40959
steps: 81919
steps: 163839
steps: 327679
https://ideone.com/8juY9Y
language:
C++14 (gcc 8.3)
created:
3 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!