Primal-Dual 原始对偶算法

想把 spfa 换成 dij，用 Johnson 里面的技巧，给予每个点一个势能 \(h_u\)，边 \((u,v,w)\) 的新边权为 \(w+h_u-h_v\)，为了保证其 \(\geq 0\) 以源点为最短路跑最短路后赋值 \(h_u\gets d_u\) 即可。

增广之后会加入反向边，考虑怎么更新势能。首先一条边的反向边被加入，其满足什么条件，然后推推式子（这里令 \(d'\) 为以 \(h\) 作为势能的最短路）：

\(d'_u+w_{u,v}+h_u-h_v=d'_v\)

将同类型的合并：

\((d'_u+h_u)-(d'_v+h_v)+w_{u,v}=0\)

有 \(w_{u,v}=-w_{v,u}\) 所以 \((d'_u+h_u)-(d'_v+h_v)=w_{v,u}\)

此时发现如果将新的势能 \(h'_u\gets d'_u+h_u\) 即满足反向边 \(=0\)，再来验证一下原先就存在的边：

根据三角形不等式 \(d'_u+w_{u,v}+h_u-h_v\geq d'_v\)

从而也有

\(w_{u,v}+(d'_u+h_u)-(d'_v+h_v)\geq 0\)

希望写得没问题/yun 不知道哪里常数写大了一点可能 至少比之前写的费用流要快。

#include<cstdio>
#include<vector>
#include<queue>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<ctime>
#include<random>
#include<assert.h>
#define pb emplace_back
#define mp make_pair
#define fi first
#define se second
#define dbg(x) cerr<<"In Line "<< __LINE__<<" the "<<#x<<" = "<<x<<'\n'
#define dpi(x,y) cerr<<"In Line "<<__LINE__<<" the "<<#x<<" = "<<x<<" ; "<<"the "<<#y<<" = "<<y<<'\n'
#define DE(fmt,...) fprintf(stderr, "Line %d : " fmt "\n",__LINE__,##__VA_ARGS__)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int>pii;
typedef pair<ll,int>pli;
typedef pair<ll,ll>pll;
typedef pair<int,ll>pil;
typedef vector<int>vi;
typedef vector<ll>vll;
typedef vector<pii>vpii;
typedef vector<pll>vpll;
template<typename T>T cmax(T &x, T y){return x=x>y?x:y;}
template<typename T>T cmin(T &x, T y){return x=x<y?x:y;}
template<typename T>
T &read(T &r){
	r=0;bool w=0;char ch=getchar();
	while(ch<'0'||ch>'9')w=ch=='-'?1:0,ch=getchar();
	while(ch>='0'&&ch<='9')r=r*10+(ch^48),ch=getchar();
	return r=w?-r:r;
}
template<typename T1,typename... T2>
void read(T1 &x,T2& ...y){read(x);read(y...);}
const int N=1000010;
const int inf=0x7fffffff;
const ll INF=0x7fffffffffffffff;
int n,m,S,T,tot;
int head[N],lst[N],pre[N],vis[N],ent=1;
ll flow[N],dis[N],h[N];
struct Edge{
	int fr,to,nxt;
	ll fl,co;
}e[N<<1];
inline void add(int x,int y,ll w,ll c){
	e[++ent]={x,y,head[x],w,c};head[x]=ent;
	e[++ent]={y,x,head[y],0,-c};head[y]=ent;
}
void spfa(){
	for(int i=1;i<=tot;i++)dis[i]=INF;
	queue<int>q;
	q.push(S);dis[S]=0;
	while(!q.empty()){
		int x=q.front();q.pop();
		vis[x]=0;
		for(int i=head[x];i;i=e[i].nxt){
			int v=e[i].to;
			if(dis[v]>dis[x]+e[i].co+h[x]-h[v]&&e[i].fl){
				dis[v]=dis[x]+e[i].co+h[x]-h[v];
				if(!vis[v]){
					q.push(v);
					vis[v]=1;
				}
			}
		}
	}
}
bool dij(){
	for(int i=1;i<=tot;i++)
		if(dis[i]!=INF)
			h[i]+=dis[i];
	for(int i=1;i<=tot;i++)lst[i]=pre[i]=vis[i]=0,dis[i]=flow[i]=INF;
	priority_queue<pll,vpll,greater<pll>>q;
	q.push(mp(0,S));dis[S]=0;
	while(!q.empty()){
		int x=q.top().se;q.pop();
		if(vis[x])continue;
		vis[x]=1;
		for(int i=head[x];i;i=e[i].nxt){
			int v=e[i].to;
			if(dis[v]>dis[x]+e[i].co+h[x]-h[v]&&e[i].fl){
				dis[v]=dis[x]+e[i].co+h[x]-h[v];
				pre[v]=x;lst[v]=i;
				flow[v]=min(flow[x],e[i].fl);
				q.push(mp(dis[v],v));
			}
		}
	}
	return pre[T]>0;
}
void MCMF(){
	spfa();
	ll mxfl=0,mnco=0;
	while(dij()){
		mxfl+=flow[T];
		mnco+=(dis[T]-h[S]+h[T])*flow[T]; 
		int now=T;
		while(now!=S){
			e[lst[now]].fl-=flow[T];
			e[lst[now]^1].fl+=flow[T];
			now=pre[now];
		}
	}
	cout << mxfl << ' ' << mnco << '\n';
}
signed main(){
	read(n,m,S,T);tot=n;
	for(int i=1;i<=m;i++){
		int u,v,w,c;read(u,v,w,c);
		add(u,v,w,c);
	}
	MCMF();
    #ifdef do_while_true
		cerr<<'\n'<<"Time:"<<clock()<<" ms"<<'\n';
	#endif
	return 0;
}