ISAP算法

为什么叫ISAP

　　ISAP(Improved Shortest Augment Path)：改进的最短增广路，属于增广路算法

算法

　　Dinic算法中，我们每次都需要BFS出层次图，而在ISAP中，我们只需要初始化时BFS出层次图（从$T$向$S$进行），然后在増广的过程中维护标号$d$（就是到汇点$T$的距离）。

　　増广的过程和Dinic类似，搜索时沿着$d_v=d_u-1$的边走，走出来的肯定是最短路。不同的是，増广完后，并不会马上对标号$d$进行更新，而是増广到没办法再増广时再更新。具体来说，就是不断BFS找增广路进行増广。如果増广过程中发现没有从$S$到$u$再到$T$的增广路，就要对$d_u$进行更新：令$d_u=\min_{(u,v)\in E'}d_v+1$，其中$E'$为残留网络的边集。这样可以保证下一次找到的増广路还是最短路。如果$u$没有出边，就令$d_u=n$，这样以后増广时就不会经过这个点。更新标号后要往回退一步。

　　那什么时候停止呢？当然是不存在$S$到$T$的最短路的时候。因为最短路是无环的，所以如果$S$到$T$有路径那么一定有$d_S\leq n-1$。

　　时间复杂度是$O(n^2m)$，然而一般不会到这个上界。这就是为什么増广路算法的时间复杂度被称为玄学。

步骤

• 初始化，BFS出层次图
• 每次沿着$d_v=d_u-1$的边进行搜索，找到一条増广路就増广
• 如果没有从$u$出发到$T$的路径，就更新$d_u$，往回退一步
• 当$d_S\geq n$时停止

一些优化

• 多次増广。每次找到一条増广路后可以往回退，尝试着找第二条増广路。

• DFS优化。这样每次找到増广路后就不用把$S$到$T$路径的信息全部重新更新一遍了。

• GAP优化。如果不存在标号为$x$的点，说明残留网络出现了断层。修改标号时如果发现当前点是最后一个$h_u=x$的点，就直接退出。可以维护一个$num$数组，$num_i$表示标号为$i$的点有多少个。

• 当前弧优化。如果一个点的标号没有被修改，那么这个点已经遍历过的边就不用再遍历了，因为不可能通过这些边再找到増广路。修改一个点的标号时顺便把遍历过的边设为空。实现时维护没遍历过得第一条边。

• 标号$d$的修改有连续性，即我们不需要找到最小的$d_v$，而是直接把$d_u+1$。

• 如果流量用完就退出。

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
int h[1000010];
int v[8000010];
ll c[8000010];
int t[8000010];
int cnt=0;
void add(int x,int y,ll z)
{
    cnt++;
    v[cnt]=y;
    c[cnt]=z;
    t[cnt]=h[x];
    h[x]=cnt;
    cnt++;
    v[cnt]=x;
    c[cnt]=0;
    t[cnt]=h[y];
    h[y]=cnt;
}
int op(int x)
{
    return ((x-1)^1)+1;
}
int S,T;
int d[1000010];
int to[1000010];
int q[1000010];
int num[1000010];
int cur[1000010];
int n,m;
void bfs()
{
    memset(d,-1,sizeof d);
    d[T]=0;
    int head=1,tail=0;
    q[++tail]=T;
    while(tail>=head)
    {
        int i,x=q[head++];
        num[d[x]]++;
        for(i=h[x];i;i=t[i])
            if(c[op(i)]&&d[v[i]]==-1)
            {
                d[v[i]]=d[x]+1;
                q[++tail]=v[i];
            }
    }
}
int stop;
ll dfs(int x,ll flow)
{
    if(x==T)
        return flow;
    ll s=0,u;
    int &i=cur[x];
    for(;i;i=t[i])
        if(d[v[i]]==d[x]-1&&c[i])
        {
            u=dfs(v[i],min(c[i],flow));
            s+=u;
            flow-=u;
            c[i]-=u;
            c[op(i)]+=u;
            if(stop)
                return s;
            if(!flow)
                return s;
        }
    num[d[x]]--;
    if(!num[d[x]])
    {
        d[S]=n;
        stop=1;
        return s;
    }
    d[x]++;
    num[d[x]]++;
    cur[x]=h[x];
    return s;
}
ll maxflow()
{
    stop=0;
    bfs();
    ll ans=0;
    memcpy(cur,h,sizeof h);
    while(d[S]<=n-1)
        ans+=dfs(S,0x7fffffffffffffffll);
    return ans;
}
int main()
{
#ifdef DEBUG
    freopen("isap.in","r",stdin);
    freopen("isap.out","w",stdout);
#endif
    scanf("%d%d%d%d",&n,&m,&S,&T);
    int i,x,y;
    ll z;
    for(i=1;i<=m;i++)
    {
        scanf("%d%d%lld",&x,&y,&z);
        add(x,y,z);
    }
    ll ans=maxflow();
    printf("%lld\n",ans);
    return 0;
}