[理性愉悦]最大流算法学习

本文总结了一些OI中常用和不常用的最大流算法，并对运行效率进行简单测试。主要是闲着无聊没事干

Tip：除了Dinic以外其他算法我都是现学现写的，理解/代码上可能有各种错误，欢迎拍砖qwq

0. Edmonds-Karp

由于这种方法实在太慢了，这里不提。

不过这个做法告诉了我们一个重要道理：网络流找增广路要加反向边。

1. Dinic

这个经典网络流算法相信大家都耳熟能详。

我们在残余网络上bfs，把访问到的点按层标号，接着将跨层的边进行dfs增广。（当然一句话讲的也不是很清楚，不会的自行百度吧）

正确的dinic分层姿势是从T开始倒着分层，从S开始正着分层复杂度似乎不变，但是从T开始跑的快一点。

Dinic主要有两个优化，一个是当前弧优化（即不重复访问一条边），另一个是当一个点无法增广时就放弃这个点。还有一个优化比较玄学，因为bfs很慢，所以我们随便搜出一条路就返回，复杂度全看脸，详见代码中的注释。

#define SZ 666666
typedef long long ll;
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],ff[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ];
bool bfs()
{
    memset(d,-1,4*(N+1));
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h!=t)
    {
        int cur=q[h++];
        for(int e=fst[cur];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==-1&&cap[e^1])
            {
                d[b]=d[cur]+1, q[t++]=b;
                if(b==S) return 1; //这是玄学 
            }
        }
    }
    return d[S]!=-1;
}
int dfs(int x,int f)
{
    if(f<=0) return 0;
    if(x==T) return f;
    int ca=0;
    for(int& e=fst[x];e;e=nxt[e]) //当前弧优化 
    {
        int b=vb[e];
        if(d[b]+1!=d[x]) continue;
        int w=dfs(b,(cap[e]<f-ca)?cap[e]:(f-ca));
        cap[e]-=w; cap[e^1]+=w; ca+=w;
        if(ca==f) break;
    }
    if(!ca) d[x]=-1; //另一个优化 
    return ca;
}
#define inf 2e9+3
int dinic()
{
    int ans=0;
    memcpy(ff,fst,4*(N+1));
    while(bfs())
    {
        //当前弧优化
        ans+=dfs(S,inf);
        memcpy(fst,ff,4*(N+1));
    }
    return ans;
}

理论时间复杂度为$O(n^2m)$，据说寻找增广路时使用LCT可以优化到$O(nmlog(n))$（）。

2. Sap

考虑仍然像dinic一样把从T开始反着把图分层，这回我们增广时采用dfs，只进行单路增广，只走相邻层的边，找到一条增广路就增广返回起点。

这个做法听上去非常naive，但是可以加入两个优化。

Dinic的时候我们每次增广完一堆路径，这时候我们会重新跑一遍bfs，然后再增广一堆路径，而Sap不需要重新跑一遍bfs，而是当某个点找不到增广路时，把它到T的层数用它所有有流量的后继层数最小值+1来更新再返回，这样就避免了重新跑bfs。

此外还有一个叫gap的玄学优化，就是说如果从S到T中间有一个层一个点也没有，那么这个层就是一个“gap”，阻止了流量的传输，那么我们就可以直接退出bfs节省时间。

当然当前弧优化也是可以使用的。实际写的时候并不需要真的dfs，而是记录每个点的pre（反正增广也要用），直接回溯。

#define SZ 666666
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],ff[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ],gap[SZ];
int pre[SZ],low[SZ],pe[SZ];
int sap()
{
    for(int i=0;i<=N;++i) d[i]=N,ff[i]=fst[i];
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h^t)
    {
        int x=q[h++];
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==N&&cap[e^1])
                d[b]=d[x]+1, q[t++]=b;
        }
    }
    for(int i=1;i<=N;++i) ++gap[d[i]];
    int x=S,ans=0;low[S]=2e9;
    while(d[S]<N)
    {
        int ce=-1;
        for(int& e=fst[x];e;e=nxt[e]) //当前弧优化 
        {
            if(!cap[e]||d[vb[e]]+1!=d[x]) continue;
            ce=e; break;
        }
        if(ce==-1)
        {
            --gap[d[x]];
            if(!gap[d[x]]) break; //gap优化 
            d[x]=N; fst[x]=ff[x]; //回退当前弧 
            for(int e=fst[x];e;e=nxt[e])
                (cap[e]&&d[vb[e]]+1<d[x])?(d[x]=d[vb[e]]+1):0; //sap
            ++gap[d[x]];
            if(x!=S) x=pre[x]; //回溯 
        }
        else
        {
            int b=vb[ce]; pre[b]=x;pe[b]=ce;
            low[b]=(low[x]<cap[ce])?low[x]:cap[ce];
            x=b; if(x!=T) continue;
            //找到了增广路
            ans+=low[x];
            for(int p=x;p!=S;p=pre[p])
                cap[pe[p]]-=low[x],cap[pe[p]^1]+=low[x];
            x=S;
        }
    }
    return ans;
}

它的理论时间复杂度仍然为$O(n^2m)$。

3. HLPP

http://blog.****.net/hechenghai/article/details/42024683

https://en.wikipedia.org/wiki/Push%E2%80%93relabel_maximum_flow_algorithm

这个算法全称叫最高标号预流推进，和Sap类似，仍然从T开始给每个点一个距离标号，无法增广时就更新标号。开始每个点的距离标号可以赋为0，也可以正常bfs。

既然叫“最高标号”预流推进算法，这个算法按标号从大到小维护了一个优先队列，表示要增广的点。

那么什么叫“预流推进”呢？“预流”就表示想往外流，每次考虑要增广一条边(a,b)，那么a的预流就减掉这条边的流量，b的预流加上这条边的流量。

这个算法主体用bfs增广，增广一个点时，考虑每条通往下层的边，按上面的做法推流，推到了没进队的点就加入优先队列，当这个点的预流为0时它的历史使命就完成了。如果一个点想要推流却啥也没推出去，那么就更新它的标号。最后检查这个点的预流是否为0，没0继续加入队列。

注意S和T除了一开始S加入队列外不应该加入队列。开始将S的预流赋为∞，该推的流都推完以后T的预流就是答案。

当前弧优化在这个算法中似乎不适用。

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define VIZ {printf("digraph G{\n"); for(int i=1;i<=n;i++) for es(i,e) printf("%d->%d;\n",i,vb[e]); puts("}");}
#define VIZ2 {printf("graph G{\n"); for(int i=1;i<=n;i++) for es(i,e) if(vb[e]>=i)printf("%d--%d;\n",i,vb[e]); puts("}");}
#define SZ 666666
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ],ex[SZ];
bool inq[SZ]; priority_queue<pii> pq;
#define gq(x) ((x)!=S&&(x)!=T&&ex[x]&&!inq[x])
int hlpp()
{
    for(int i=0;i<=N;++i) d[i]=N+2;
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h^t)
    {
        int x=q[h++];
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==N+2&&cap[e^1])
                d[b]=d[x]+1, q[t++]=b;
        }
    }
    d[S]=N+2;
    pq.push(pii(d[S],S)); inq[S]=1; ex[S]=2e9;
    while(!pq.empty())
    {
        int x=pq.top().se; pq.pop(); inq[x]=0;
        bool p=0;
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(!cap[e]||(x!=S&&d[b]+1!=d[x])) continue;
            p=1; int f=(cap[e]<ex[x])?cap[e]:ex[x];
            ex[x]-=f; ex[b]+=f;
            cap[e]-=f; cap[e^1]+=f;
            if(gq(b)) inq[b]=1,pq.push(pii(d[b],b));
            if(!ex[x]) break;
        }
        if(!p)
        {
            d[x]=N+2;
            for(int e=fst[x];e;e=nxt[e])
                (cap[e]&&d[vb[e]]+1<d[x])?(d[x]=d[vb[e]]+1):0;
        }
        if(gq(x)) inq[x]=1,pq.push(pii(d[x],x));
    }
    return ex[T];
}
namespace FF
{
char ch,B[1<<15],*S=B,*T=B;
#define getc() (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
#define isd(c) (c>='0'&&c<='9')
int aa,bb;int F(){
    while(ch=getc(),!isd(ch)&&ch!='-');ch=='-'?aa=bb=0:(aa=ch-'0',bb=1);
    while(ch=getc(),isd(ch))aa=aa*10+ch-'0';return bb?aa:-aa;
}
}
#define gi FF::F()
int main()
{
    N=gi; int m=gi;
    S=1; T=N;
    for(int i=1;i<=m;i++)
    {
        int x=gi,y=gi,z=gi;
        ad_dl(x,y,z);
    }
    printf("%d\n",hlpp());
}

UPD：注意！这份代码是错误的，原因我现在还不清楚，反正是错的，下面有关HLPP的测试结果请忽略。

理论时间复杂度比上面两个优秀！为$O(n^2\sqrt{m})$，这个界是紧的。

4. Goldberg-Rao

http://www.cs.princeton.edu/courses/archive/fall06/cos528/handouts/Goldberg-Rao.pdf

http://www.cs.tau.ac.il/~haimk/adv-alg-2013/goldbergrao.pdf

https://people.orie.cornell.edu/dpw/orie633/LectureNotes/lecture10.pdf

首先我们需要知道，dinic的复杂度十分玄学，在容量为1的图上跑的飞快，经过一番严谨的证明可以证出上界是$O(min(n^{2/3},m^{1/2})m)$（详见维基）。

我们考虑将dinic算法改造一番将它在容量大的时候也跑的很快。

在开始之前，我们设$\Lambda=min(n^{2/3},m^{1/2})$，最多能增广的量不超过F，$\Delta=\lceil{F}/\Lambda\rceil$。

考虑形式化地定义dinic算法，我们定义一条边u->v是好的，当且仅当$d(v)=d(u)+l$，其中l表示这条边的长度。在原始dinic算法中，l恒为1，但在这个算法中，l可能为0或1。d的定义同理，形式化地，$d(T)=0$，对于每个点v，$d(v)\leq{d(w)+l(v~to~w)}$，并且等号能取到。

我们定义一条边为大的，当且仅当它的残余容量$\geq3\Delta$；一条边为特别的，这个条件比较多，有三个条件：残余容量在$[2\Delta,3\Delta)$，两端点距离标号相等，它的反向边残余容量$\geq3\Delta$。

那么我们现在可以定义一条边的l了，定义十分简单，如果这条边是大的或者特别的，l为0，否则为1。

仔细一想你会觉得有些不对，l依赖“特别”来定义，“特别”依赖距离标号来定义，而距离标号却又依赖l来计算...但是事实上我们可以先忽略特别的边，即只考虑小的边将其l设为0，这样计算出距离标号之后再重新计算一遍l，因为特别的边显然不会改变距离标号。

然后我们开始先定一个F的上界，为了简化我们设U为每条边的最大容量，那么就开始设$F=mU$。

算法流程大致如下：

初始化F

while F>=1:
  计算△
  重复5Λ次操作:
    初始化距离标号和边的l值
    在好的边上找一个堵塞的或者流量为△的流
    增广
  F/=2

“堵塞”的流就表示是增广出来满的一个流。是不是看起来十分简单呢？

我们还遗漏了一个重要细节。在正常dinic中，l都为1，增广轻松愉快。但在这个算法中l可能等于0，那么有可能出现流量循环的情况，因为同一层的点间可能会有流量。

流量循环？Tarjan！我们在好的边上跑一遍tarjan，把每个强连通分量缩成一个点，那么就不会有什么循环的情况出现了！

那么缩完点之后的流量怎么办？我们在强连通分量里随便找个点r当做根，同样在好边上dfs两遍搞出一棵In Tree，一棵Out Tree。盗图一张：

如果有了流量，我们通过In Tree上的边流给r，再从Out Tree上的边流出去，这样就搞完啦！

不过还有一个小优化，考虑每次F/=2显得非常蠢，我们考虑想办法乱搞求出残余流量的一个上界，进行估价，每次F=min(F/2,计算出的上界)。

我们考虑最小割=最大流，那么考虑算出一个比较小的割。考虑在原图上bfs求出距离标号，然后每次考虑对于每个k，把距离标号为k和k+1的点之间的边全部割掉，这显然是一个合法的割，虽然不是很小，我们就用这个进行估价。同时，开始的F也可以设为这个值。

当前弧优化还是可以继续使用的。另外如果找到的堵塞流流量为0，可以直接break。

正确性和复杂度证明留作习题

注意下面这份代码可能有各种各样的错误，数据中的边权都是随机的，所以增广一次就够了...更强的数据我并不会造，就先这样好了...

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define VIZ {printf("digraph G{\n"); for(int i=1;i<=n;i++) for es(i,e) printf("%d->%d;\n",i,vb[e]); puts("}");}
#define VIZ2 {printf("graph G{\n"); for(int i=1;i<=n;i++) for es(i,e) if(vb[e]>=i)printf("%d--%d;\n",i,vb[e]); puts("}");}
#define SZ 666666
#define FEdge \
int M=1,fst[SZ],vb[SZ],cap[SZ],nxt[SZ],va[SZ];\
void ad_de_(int a,int b,int c)\
{++M; nxt[M]=fst[a]; vb[M]=b; cap[M]=c; fst[a]=M; va[M]=a;}\
void ad_de(int a,int b,int c)\
{ad_de_(a,b,c); ad_de_(b,a,0);}
#define GetD(l) \
int qs[SZ];\
bool inq[SZ];\
void spfa()\
{\
    int h=0,t=1;qs[0]=T;\
    for(int i=1;i<=n;i++) d[i]=1e9;\
    d[T]=0; inq[T]=1;\
    while(h^t)\
    {\
        int x=qs[h++];\
        for(int e=fst[x];e;e=nxt[e])\
        {\
            if(cap[e^1]);else continue;\
            int y=d[x]+l(e^1),b=vb[e];\
            if(d[b]<=y) continue;\
            d[b]=y; if(inq[b]) continue;\
            inq[b]=1; qs[t++]=b;\
        }\
        inq[x]=0;\
    }\
}
int n,m,S,T;
namespace dn
{
int T,d[SZ],M=1,fst[SZ],vb[SZ],cap[SZ],nxt[SZ],l[SZ];
int addedge(int a,int b,int c,int g)
{++M; nxt[M]=fst[a]; vb[M]=b; cap[M]=c; fst[a]=M; l[M]=g; return M;}
int dfs(int x,int f)
{
    if(f<=0) return 0;
    if(x==T) return f;
    int ca=0;
    for(int& e=fst[x];e;e=nxt[e])
    {
        int b=vb[e];
        if(d[b]+l[e]!=d[x]) continue;
        int w=dfs(b,min(cap[e],f-ca));
        cap[e]-=w; cap[e^1]+=w; ca+=w;
        if(ca==f) break;
    }
    if(!ca) d[x]=-233;
    return ca;
}
}
namespace ng
{
FEdge int d[SZ];
namespace og
{
int d[SZ],sm[SZ]; GetD(1+0*)
int gj()
{
    spfa();
    int mx=d[S];
    if(mx==1e9) return 0;
    for(int i=0;i<=mx;++i) sm[i]=0;
    for(int i=1;i<=n;i++)
    {
        for(int e=fst[i];e;e=nxt[e])
        {
            int b=vb[e];
            if(!cap[e]||d[b]+1!=d[i]) continue;
            sm[d[i]]+=cap[e];
        }
    }
    int ans=2147483647;
    for(int i=1;i<=mx;++i)
        ans=min(ans,sm[i]);
    return ans;
}
}
ll F; int L,D,vs[SZ],ans=0;
bool bg(int e) {return cap[e]>=D+D+D;}
bool fb(int e) {return cap[e]<D+D+D;}
int val(int e)
{return !(bg(e)||(cap[e]>=D+D&&cap[e]<D+D+D
&&d[va[e]]==d[vb[e]]&&cap[e^1]>=D+D+D));}
GetD(fb)
int low[SZ],dfn[SZ],cur=0;
bool ins[SZ];
int ss[SZ],sn=0,bl[SZ],C=0;
void tarjan(int x)
{
    dfn[x]=low[x]=++cur; ss[++sn]=x; ins[x]=1;
    for esb(x,e,b)
    {
        if(vs[e]) continue;
        if(ins[b]) low[x]=min(low[x],dfn[b]);
        else if(!dfn[b])
            tarjan(b),low[x]=min(low[x],low[b]);
    }
    if(dfn[x]!=low[x]) return;
    int cs=0; ++C;
    while(1)
    {
        ++cs; int t=ss[sn--];
        bl[t]=C; ins[t]=0;
        if(t==x) break;
    }
}
int fa[SZ],tmp[SZ],oi[SZ]/*in-out*/;
void dfs(int x,int g)
{
    for esb(x,e,b)
    {
        if(bl[b]!=g||fa[b]) continue;
        fa[b]=x; dfs(b,g);
        cap[e]+=oi[b]; cap[e^1]-=oi[b];
        oi[x]+=oi[b]; oi[b]=0;
    }
}
//大概有错 
void doit()
{
    L=min(pow(n,2.0/3),pow(m,1.0/2));
    if(L<1) L=1; F=1e9;
    int tm=0;
    while(F>=1)
    {
        D=ceil(F/(ld)L);
        for(int i=1;i<=5*L;++i)
        {
            spfa();
            for(int i=2;i<=M;++i) vs[i]=val(i);
            sn=cur=C=0;
            for(int i=1;i<=n;i++)
                dfn[i]=low[i]=ins[i]=oi[i]=0;
            for(int i=1;i<=n;i++)
                if(!low[i]) tarjan(i);
            if(bl[S]==bl[T]) break;
            dn::M=1; dn::T=bl[T];
            for(int i=1;i<=C;i++) dn::fst[i]=0;
            for(int i=1;i<=n;i++) dn::d[bl[i]]=d[i];
            for(int i=2;i<=M;i++)
            {
                int a=va[i],b=vb[i];
                if(bl[a]!=bl[b]) tmp[i]=dn::addedge(bl[a],bl[b],cap[i],vs[i]);
            }
            ll cur=dn::dfs(bl[S],D); ans+=cur;
            if(!cur) break;
            for(int i=1;i<=n;i++)
                for esb(i,e,b)
                    if(bl[i]!=bl[b])
                    {
                        ll l=cap[e]-dn::cap[tmp[e]];
                        oi[i]-=l; cap[e]-=l;
                    }
            for(int i=1;i<=n;i++)
                if(!low[i]) fa[i]=-1,dfs(i,bl[i]);
        }
        F=min(int(F/2),og::gj());
    }
    printf("%d\n",ans);
}
}
namespace FF
{
char ch,B[1<<15],*S=B,*T=B;
#define getc() (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
#define isd(c) (c>='0'&&c<='9')
int aa,bb;int F(){
    while(ch=getc(),!isd(ch)&&ch!='-');ch=='-'?aa=bb=0:(aa=ch-'0',bb=1);
    while(ch=getc(),isd(ch))aa=aa*10+ch-'0';return bb?aa:-aa;
}
}
#define gi FF::F()
int main()
{
    n=gi,m=gi; S=1; T=n;
    for(int i=1;i<=m;++i)
    {
        int a,b,c;
        a=gi,b=gi,c=gi;
        ng::ad_de(a,b,c);
    }
    ng::doit();
}

理论复杂度为$O(min(n^{2/3},m^{1/2})mlog(n)log(mU))$，其中U表示边的最大容量。

测试数据中n=3000，m=200000，11组数据，时限2s，采用各种姿势生成，目测挺强的。如果按理论复杂度这些算法一个也过不了

sap：100分/492ms，dinic：100分/508ms，假的goldberg-rao：100分/1812ms，HLPP：45~54分（写挂了）。

终于填上了这个大坑好高兴啊，过几天开新坑。

[理性愉悦]最大流算法学习

相关推荐