题解 P3381 【【模板】最小费用最大流】
远航之曲
·
·
题解
不知道luogu的LaTeX支持如何。。。要是崩了的话可以去我的blog看看
这里提供一种新的思想,用Dij求解费用流
有人可能会说,Dijkstra只能在正权图上求最短路呀,而费用流的网络构图中w[x][y] = -w[y][x],势必会产生负权边,怎么能用Dijkstra呢？(下文用w[i][j]表示费用,p[i]表示点,方便讲解)
我们求解费用流一般用spfa或zkw流,这两种方法在稀疏图上运行效率还不错,但是在稠密图上的表现却无法令人满意,而且其运行效率与图的形态关系太大,无法令人安心。
反观Dijkstra,它在稠密图O(n^2)和稀疏图O(mlogm)上的表现都不错,而且其运行时间与图的形态关系不大,较为稳定,唯一的缺点是不能运行在负权图上。
事实上在对网络进行一些处理后,Dijkstra是可以用来求解费用流的。
我们对网络G中的每个点i设置一个势函数h[i],在算法运行的每一时刻,对于每条残余网络中的边(x, y),都要满足h[x] + w[x][y] - h[y] >= 0(三角不等式)。维护h[i]的方法下文再讨论。
构建新网络G',V' = V,w'[x][y] = h[x] + w[x][y] - h[y],这样在新网络中所有边的权值均非负,可以使用Dijkstra求解最短路。
$$\begin{aligned}
w'[path]&=w'[p[1]][p[2]] + w'[p[2]][p[3]] + \cdots + w'[p[m-1]][p[m]]\\
&=h[p[1]]-h[p[2]] + w[p[1]][p[2]]\\
&+h[p[2]]-h[p[3]] + w[p[2]][p[3]]\\
&+\cdots\\
&+h[p[m-1]]-h[p[m]] + w[p[m-1]][p[m]]\\
&=w[p[1]][p[2]] + w[p[2]][p[3]] + \cdots + w[p[m-1]][p[m]] + h[p[1]]-h[p[m]] \\
&=w[path] + h[p[1]]-h[p[m]]
\end{aligned}$$
而当p[1]和p[m]确定时,最优化$w[path]$与最优化$w'[path]$是等价的,所以我们可以在G'中求解S到T的最短路$w'[path]$,G中S到T的最短路$w[path] = w'[path] - h[S] + h[T]$。
下面我们来讨论如何维护势函数$h[]$。
初始情况下,如果所有的权值都为正,那么可以简单地将所有$h[i]$设置为0；如果有负权值,那么我们做一遍spfa,让h[]等于距离函数。如果初始网络为DAG,也可以采用递推$h[i] = min(h[j] + w[j][i], g[j][i] = true)$。
让$d'[i]$表示G'中S到点i的距离,当某次增广结束后,G'中会新加入某些边$(j, i)$,而这些(j, i)必定满足$d'[i] + w'[i][j] = d'[j]$(否则(i, j)就不会在增广路中)。对上式进行一些代数变换：
$$\begin{aligned}
d'[i] + w'[i][j] = d'[j]
\\d'[i] + w[i][j] + h[i]-h[j] = d'[j]
\\(d'[j] + h[j])-(d'[i] + h[i])-w[i][j] = 0
\\(d'[j] + h[j])-(d'[i] + h[i]) + w[j][i] = 0
\end{aligned}$$
(因为是费用流,所以有w[i][j] = -w[j][i])
因此让所有$h[i] += d'[i]$后,新加入的边(j, i)也会满足势函数的性质。
同时我们有：
$$\begin{aligned}
d'[i] + w'[i][j] - d'[j] >= 0\\
d'[i] + h[i] - h[j] + w[i][j] - d'[j] >= 0\\
(d'[i] + h[i]) - (d'[j] + h[j]) + w[i][j] >= 0\\
\end{aligned}$$
因此修改$h[i]$后,$(i, j)$依然会满足势函数的性质。
算法过程如下：
>S1 初始化h[]
>S2 在残留网络中做Dijkstra
>S3 若S到T有可行路径,则修改增广路上的边的容量并所有h[i] += d'[i],转S2,否则退出
算法时间复杂度为：$O(spfa() + K * Dijkstra())$或$O(K * Dijkstra())$,这取决于初始化h[]时是否调用spfa,K表示增广次数。
```cpp
#include <cstdio>
#include <algorithm>
#include <queue>
#include <cstring>
#define M 200100
#define N 10010
#define INF 0x33333333
#define min(x,y) ((x<y)?(x):(y))
using namespace std;
typedef pair<int,int> Pair;
struct node
{
int from,to,next,flow,cost;
}e[M];
int tot=-1,st[M];
int n,m,x,y,z;
void add(int x,int y,int z,int zz)
{
e[++tot].to=y;
e[tot].from=x;
e[tot].flow=z;
e[tot].cost=zz;
e[tot].next=st[x];
st[x]=tot;
}
Pair main_pro(int s,int t)
{
static int h[N];
int flow=0,cost=0;
while(1)
{
static int dis[N],pv[N],pe[N];
memset(dis,0x33,sizeof dis);
dis[s]=0;
priority_queue<Pair,vector<Pair>,greater<Pair> >que;
que.push(Pair(0,s));
while(!que.empty())
{
Pair now=que.top();
que.pop();
if (now.first!=dis[now.second]) continue;
if (now.second==t) break;
for (int i=st[now.second];~i;i=e[i].next)
{
int nowcost=e[i].cost+h[now.second]-h[e[i].to];
if (e[i].flow>0 && dis[e[i].to]>dis[now.second]+nowcost)
{
dis[e[i].to]=dis[now.second]+nowcost;
que.push(Pair(dis[e[i].to],e[i].to));
pv[e[i].to]=now.second;
pe[e[i].to]=i;
}
}
}
if (dis[t]==INF)
break;
for (int i=0;i<n;i++)
h[i]=min(h[i]+dis[i],INF);
int newflow=INF;
for (int x=t;x!=s;x=pv[x])
newflow=min(newflow,e[pe[x]].flow);
flow+=newflow;
cost+=newflow*h[t];
for (int x=t;x!=s;x=pv[x])
e[pe[x]].flow-=newflow,e[pe[x]^1].flow+=newflow;
}
return make_pair(flow,cost);
}
int main()
{
int m, from, to;
scanf("%d%d%d%d",&n,&m,&from,&to);
memset(e,-1,sizeof e);
memset(st,-1,sizeof st);
for(int i = 0; i < m; ++i)
{
int u,v,flow,cost;
scanf("%d%d%d%d",&u,&v,&flow,&cost);
add(u,v,flow,cost);
add(v,u,0,-cost);
}
Pair ans = main_pro(from, to);
printf("%d %d\n",ans.first,ans.second);
}
```