网络流模板-网络单纯形

最小费用最大流，但能过 HLPP 板子题，还能处理负环

namespace flow{ // }{{{
typedef long long ll;
constexpr int V = 5e3, E = 5e4;
constexpr int EDGE_NIL = 0;
struct Edge{
	int to;
	ll lf, cost;
	int nxt;
} es[E*2+4];
ll sumcost = 0, sumflow = 0;
int is, it, iv;
ll minc, maxf;
int head[V], cnt = (EDGE_NIL|1)+1;
ll pi[V]; 
int fe[V], mark[V], time = 2; // fe: father edge
int fa[V]; // 200ms improve
void init(int v, int s, int t){
	is = s, it = t, iv = v;
	sd fill(head, head+iv, EDGE_NIL);
	sd fill(mark, mark+iv, 0);
	sd fill(fe, fe+iv, EDGE_NIL);
	time = 2;
	cnt = (EDGE_NIL|1)+1;
	sumcost = sumflow = 0;
	minc = maxf = 0;
}
void addflow(int s, int t, ll f, ll c){
	es[cnt] = (Edge){t, f, c, head[s]}, head[s] = cnt++;
	es[cnt] = (Edge){s, 0, -c, head[t]}, head[t] = cnt++;
	sumflow += f, sumcost += sd abs(c);
}
void mktree(int x, int from_e){
	fe[x] = from_e;
	fa[x] = es[from_e^1].to;
	mark[x] = 1;
	for(int i=head[x]; i!=EDGE_NIL; i=es[i].nxt){
		if(mark[es[i].to] == 1 || es[i].lf == 0) continue;
		mktree(es[i].to, i);
	}
}
ll getpi(int x){
	if(mark[x] == time) return pi[x];
	mark[x] = time;
	pi[x] = getpi(fa[x]) - es[fe[x]].cost;
	return pi[x];
}
ll pushflow(int e){ // return delta-cost
	int rt = es[e].to, lca = es[e^1].to;
	time++;
	while(rt){
		mark[rt] = time;
		rt = fa[rt];
	}
	while(mark[lca] != time){
		mark[lca] = time;
		lca = fa[lca];
	}
	ll df = es[e].lf;
   	int todel = e, dir = -1; // dir: direction, 0->es[e].to
	for(int i=es[e^1].to; i!=lca; i=fa[i]){
		if(es[fe[i]].lf < df){
			df = es[fe[i]].lf;
			todel = fe[i];
			dir = 1;
		}
	}
	for(int i=es[e].to; i!=lca; i=fa[i]){
		if(es[fe[i]^1].lf < df){
			df = es[fe[i]^1].lf;
			todel = fe[i];
			dir = 0;
		}
	}
	ll dcst = 0; // delta cost
	if(df) {
		for(int i=es[e].to; i!=lca; i=fa[i]){
			es[fe[i]].lf += df;
			es[fe[i]^1].lf -= df;
			dcst += es[fe[i]^1].cost * df;
		}
		for(int i=es[e^1].to; i!=lca; i=fa[i]){
			es[fe[i]].lf -= df;
			es[fe[i]^1].lf += df;
			dcst += es[fe[i]].cost * df;
		}
		es[e].lf -= df;
		es[e^1].lf += df;
		dcst += es[e].cost * df;
	}
	if(todel == e) return dcst;
	int last = e^dir, lastu = es[e^dir^1].to;
	for(int i=es[e^dir].to; i!=es[todel^1].to; ){
		mark[i]=time-1;
		int i_ = fa[i];
		fa[i] = lastu;
		lastu = i;
		sd swap(fe[i], last);
		last ^= 1;
		i=i_;
	}
	return dcst;
}
void mcmf(){
	ll sfl_ = sumflow, scs_ = sumcost;
	addflow(it, is, sumflow+1, -(sumcost+1));
	sumflow = sfl_, sumcost = scs_;
	mktree(it, EDGE_NIL);
	mark[it] = ++time;
	fa[it] = 0;
	bool run = true;
	while(run){
		run = false;
		UP(i, (EDGE_NIL|1)+1, cnt){
			int s = es[i^1].to, t = es[i].to;
			if(es[i].lf && es[i].cost + getpi(t) - getpi(s) < 0){
				run = true;
				minc += pushflow(i);
			}
		}
	}
	maxf = es[cnt-1].lf;
	minc += maxf * (scs_+1);
}
} // {}}}