Update max flow code, fix MCMF

Author: ekzhang

maxflow.cpp

@@ -9,6 +9,8 @@ using namespace std;
  *   - O(V^2 E) for general graphs, but in practice ~O(E^1.5)
  *   - O(V^(1/2) E) for bipartite matching
  *   - O(min(V^(2/3), E^(1/2)) E) for unit capacity graphs
+ * 
+ * WARNING: prefer maxflow2.cpp when possible (which is faster)
  */
 template<int V, class T=long long>
 class max_flow {

maxflow2.cpp

@@ -58,8 +58,8 @@ class max_flow {
 
 public:
 	void add(int u, int v, T cap=1, T rcap=0) {
-		adj[u].push_back({ v, adj[v].size(), cap, 0 });
-		adj[v].push_back({ u, adj[u].size() - 1, rcap, 0 });
+		adj[u].push_back({ v, (int) adj[v].size(), cap, 0 });
+		adj[v].push_back({ u, (int) adj[u].size() - 1, rcap, 0 });
 	}
 
 	T calc(int s, int t) {

mincost_maxflow.cpp

@@ -6,105 +6,113 @@ using namespace std;
 
 /* Minimum-Cost, Maximum-Flow solver using Successive Shortest Paths with Dijkstra and SPFA-SLF.
  * Requirements:
- *   - No duplicate or antiparallel edges with different costs.
+ *   - Duplicate or antiparallel edges with different costs are allowed.
  *   - No negative cycles.
- * Time Complexity: O(Ef lg V) average-case, O(VE + Ef lg V) worst-case.
+ * Time Complexity: O(Ef lg V) average-case, O(VE + Ef lg V) worst-case with negative costs.
  */
 template<int V, class T=long long>
 class mcmf {
-	unordered_map<int, T> cap[V], cost[V];
+	/* making this static breaks compilation on -O0, but not on -O2; possible bug in GCC */
+	const T INF = numeric_limits<T>::max();
+
+	struct edge {
+		int t, rev;
+		T cap, cost, f;
+	};
+
+	vector<edge> adj[V];
 	T dist[V];
 	int pre[V];
-	bool visited[V];
+	bool vis[V];
 
-	void spfa(int s) {
-		static list<int> q;
+	void spfa(int s) { /* only needed if there are negative costs */
+		list<int> q;
 
 		memset(pre, -1, sizeof pre);
-		memset(dist, 63, sizeof dist);
-		memset(visited, 0, sizeof visited);
+		memset(vis, 0, sizeof vis);
+		fill(dist, dist + V, INF);
 
 		dist[s] = 0;
 		q.push_back(s);
 		while (!q.empty()) {
 			int v = q.front();
 			q.pop_front();
-			visited[v] = false;
-			for (auto p : cap[v]) if (p.second) {
-				int u = p.first;
-				T d = dist[v] + cost[v][u];
+			vis[v] = false;
+			for (auto e : adj[v]) if (e.cap != e.f) {
+				int u = e.t;
+				T d = dist[v] + e.cost;
 				if (d < dist[u]) {
-					dist[u] = d, pre[u] = v;
-					if (!visited[u]) {
+					dist[u] = d, pre[u] = e.rev;
+					if (!vis[u]) {
 						if (q.size() && d < dist[q.front()]) q.push_front(u);
 						else q.push_back(u);
-						visited[u] = true;
+						vis[u] = true;
 					}
 				}
 			}
 		}
 	}
 
-	void dijkstra(int s) {
-		static priority_queue<pair<T, int>, vector<pair<T, int> >,
-				greater<pair<T, int> > > pq;
+	priority_queue<pair<T, int>, vector<pair<T, int> >,
+		greater<pair<T, int> > > pq; /* for dijkstra */
 
+	void dijkstra(int s) {
 		memset(pre, -1, sizeof pre);
-		memset(dist, 63, sizeof dist);
-		memset(visited, 0, sizeof visited);
+		memset(vis, 0, sizeof vis);
+		fill(dist, dist + V, INF);
 
 		dist[s] = 0;
-		pq.push({0, s});
+		pq.emplace(0, s);
 		while (!pq.empty()) {
 			int v = pq.top().second;
 			pq.pop();
-			if (visited[v]) continue;
-			visited[v] = true;
-			for (auto p : cap[v]) if (p.second) {
-				int u = p.first;
-				T d = dist[v] + cost[v][u];
+			if (vis[v]) continue;
+			vis[v] = true;
+			for (auto e : adj[v]) if (e.cap != e.f) {
+				int u = e.t;
+				T d = dist[v] + e.cost;
 				if (d < dist[u]) {
-					dist[u] = d, pre[u] = v;
-					pq.push({d, u});
+					dist[u] = d, pre[u] = e.rev;
+					pq.emplace(d, u);
 				}
 			}
 		}
 	}
 
 	void reweight() {
-		for (int v = 0; v < V; v++) {
-			for (auto& p : cost[v]) {
-				p.second += dist[v] - dist[p.first];
-			}
-		}
+		for (int v = 0; v < V; v++)
+			for (auto& e : adj[v])
+				e.cost += dist[v] - dist[e.t];
 	}
 
 public:
 	unordered_map<int, T> flows[V];
 
-	void add(int u, int v, T f=1, T c=0) {
-		cap[u][v] += f;
-		cost[u][v] = c;
-		cost[v][u] = -c;
+	void add(int u, int v, T cap=1, T cost=0) {
+		adj[u].push_back({ v, (int) adj[v].size(), cap, cost, 0 });
+		adj[v].push_back({ u, (int) adj[u].size() - 1, 0, -cost, 0 });
 	}
 
 	pair<T, T> calc(int s, int t) {
-		spfa(s);
+		spfa(s); /* comment out if all costs are non-negative */
 		T totalflow = 0, totalcost = 0;
 		T fcost = dist[t];
 		while (true) {
 			reweight();
 			dijkstra(s);
 			if (~pre[t]) {
 				fcost += dist[t];
-				T flow = cap[pre[t]][t];
-				for (int v = t; ~pre[v]; v = pre[v])
-					flow = min(flow, cap[pre[v]][v]);
-				for (int v = t; ~pre[v]; v = pre[v]) {
-					cap[pre[v]][v] -= flow;
-					cap[v][pre[v]] += flow;
-					flows[pre[v]][v] += flow;
-					flows[v][pre[v]] -= flow;
+				T flow = INF;
+				for (int v = t; ~pre[v]; v = adj[v][pre[v]].t) {
+					edge& r = adj[v][pre[v]];
+					edge& e = adj[r.t][r.rev];
+					flow = min(flow, e.cap - e.f);
+				}
+				for (int v = t; ~pre[v]; v = adj[v][pre[v]].t) {
+					edge& r = adj[v][pre[v]];
+					edge& e = adj[r.t][r.rev];
+					e.f += flow;
+					r.f -= flow;
 				}
 				totalflow += flow;
 				totalcost += flow * fcost;
@@ -116,10 +124,8 @@ class mcmf {
 
 	void clear() {
 		for (int i = 0; i < V; i++) {
-			cap[i].clear();
-			cost[i].clear();
-			flows[i].clear();
-			dist[i] = pre[i] = visited[i] = 0;
+			adj[i].clear();
+			dist[i] = pre[i] = vis[i] = 0;
 		}
 	}
 };

treap.cpp

@@ -15,7 +15,7 @@ struct tnode {
 	int sz;
 
 	tnode() { update(); }
-	
+
 	void update() {
 		sz = 1 + (l ? l->sz : 0) + (r ? r->sz : 0);
 	}