ISAP模板分析

概述：算法基于这样的一个事实：每次增广之后，任意结点到汇点（在残余网络中）的最短距离都不会减小。这样，我们可以利用d[i[表示结点i到汇点的距离的下界。然后再增广过程当中不断地修改这个下界。增广的时候和Dinic算法类似，只允许沿着d[i]==d[j]+1的弧(i,j)走。

不难证明，d[i[满足两个条件：(1)d[t]=0;(2)对任意的弧(i,j) d[i]≤d[j]+1。因为最坏的情况就是s到t是一条链，此时等号成立。因此，当d[s]≥n时，残余网络中不存在s-t路。

那么与Dinic算法类似，事先逆向bfs，找增广的过程就是沿着“允许弧”（即满足f< c且d[i]==d[j]+1的弧）往前走。（称为“前进”）。如果向前走不动了，那么就要考虑原路返回（称为“撤退”）。此时把d[i]修改为min{d[j]}+1即可。因为要满足d函数的条件(2)。修改后，原来的i值的个数就减少一个，而新i值的个数多一个。在程序中，用num数组来保存所有距离的个数，当把距离值从x修改为y时，num[x]–,num[y]++即可，然后检查num[x]是否为0，如果是0，那么s-t不连通，算法终止。原因显而易见：比如s-t的距离是3，如果距离为2的情况都已经没了，更别提走到距离为1的点了。这就是所谓的“gap优化”。

通过之前的分析，在数据结构方面，该算法只比Dinic算法的数据结构多了两个数组：用于记录父边以便于撤退的数组p，以及标记距离个数的数组num。增广的时候分为两步，第一步逆推求出可改进量a（即残余量的最小值）；第二步再逆推一遍，进行增广。主过程中，x走到汇点时增广。

最快的模板：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<string>
#include<vector>
#include<stack>
#include<bitset>
#include<cstdlib>
#include<cmath>
#include<set>
#include<list>
#include<deque>
#include<map>
#include<queue>
#define MAXN 100010  //点
#define MAXM 800010//边
#define inf 0x3f3f3f
using namespace std;
struct Node
{
    int from,to,next;
    int cap;
}edge[MAXM];
int tol;
int head[MAXN];
int dep[MAXN];
int gap[MAXN];//gap[x]=y :说明残留网络中dep[i]==x的个数为y

int n;//n是总的点的个数，包括源点和汇点

void init()
{
    tol=0;
    memset(head,-1,sizeof(head));
}

void addedge(int u,int v,int w)
{
    edge[tol].from=u;
    edge[tol].to=v;
    edge[tol].cap=w;
    edge[tol].next=head[u];
    head[u]=tol++;
    edge[tol].from=v;
    edge[tol].to=u;
    edge[tol].cap=0;
    edge[tol].next=head[v];
    head[v]=tol++;
}
void BFS(int start,int end)
{
    memset(dep,-1,sizeof(dep));
    memset(gap,0,sizeof(gap));
    gap[0]=1;
    int que[MAXN];
    int front,rear;
    front=rear=0;
    dep[end]=0;
    que[rear++]=end;
    while(front!=rear)
    {
        int u=que[front++];
        if(front==MAXN)front=0;
        for(int i=head[u];i!=-1;i=edge[i].next)
        {
            int v=edge[i].to;
            if(dep[v]!=-1)continue;
            que[rear++]=v;
            if(rear==MAXN)rear=0;
            dep[v]=dep[u]+1;
            ++gap[dep[v]];
        }
    }
}
int ISAP(int start,int end)
{
    int res=0;
    BFS(start,end);
    int cur[MAXN];
    int S[MAXN];
    int top=0;
    memcpy(cur,head,sizeof(head));
    int u=start;
    int i;
    while(dep[start]<n)
    {
        if(u==end)
        {
            int temp=inf;
            int inser;
            for(i=0;i<top;i++)
               if(temp>edge[S[i]].cap)
               {
                   temp=edge[S[i]].cap;
                   inser=i;
               }
            for(i=0;i<top;i++)
            {
                edge[S[i]].cap-=temp;
                edge[S[i]^1].cap+=temp;
            }
            res+=temp;
            top=inser;
            u=edge[S[top]].from;
        }
        if(u!=end&&gap[dep[u]-1]==0)//出现断层，无增广路
          break;
        for(i=cur[u];i!=-1;i=edge[i].next)
           if(edge[i].cap!=0&&dep[u]==dep[edge[i].to]+1)
             break;
        if(i!=-1)
        {
            cur[u]=i;
            S[top++]=i;
            u=edge[i].to;
        }
        else
        {
            int min=n;
            for(i=head[u];i!=-1;i=edge[i].next)
            {
                if(edge[i].cap==0)continue;
                if(min>dep[edge[i].to])
                {
                    min=dep[edge[i].to];
                    cur[u]=i;
                }
            }
            --gap[dep[u]];
            dep[u]=min+1;
            ++gap[dep[u]];
            if(u!=start)u=edge[S[--top]].from;
        }
    }
    return res;
}

int source;         // 源点
int sink;           // 汇点
int p[max_nodes];   // 可增广路上的上一条弧的编号
int num[max_nodes]; // 和 t 的最短距离等于 i 的节点数量
int cur[max_nodes]; // 当前弧下标
int d[max_nodes];   // 残量网络中节点 i 到汇点 t 的最短距离
bool visited[max_nodes];

// 预处理, 反向 BFS 构造 d 数组
bool bfs()
{
    memset(visited, 0, sizeof(visited));
    queue<int> Q;
    Q.push(sink);
    visited[sink] = 1;
    d[sink] = 0;
    while (!Q.empty()) {
        int u = Q.front();
        Q.pop();
        for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix) {
            Edge &e = edges[(*ix)^1];
            if (!visited[e.from] && e.capacity > e.flow) {
                visited[e.from] = true;
                d[e.from] = d[u] + 1;
                Q.push(e.from);
            }
        }
    }
    return visited[source];
}

// 增广
int augment()
{
    int u = sink, df = __inf;
    // 从汇点到源点通过 p 追踪增广路径, df 为一路上最小的残量
    while (u != source) {
        Edge &e = edges[p[u]];
        df = min(df, e.capacity - e.flow);
        u = edges[p[u]].from;
    }
    u = sink;
    // 从汇点到源点更新流量
    while (u != source) {
        edges[p[u]].flow += df;
        edges[p[u]^1].flow -= df;
        u = edges[p[u]].from;
    }
    return df;
}

int max_flow()
{
    int flow = 0;
    bfs();
    memset(num, 0, sizeof(num));
    for (int i = 0; i < num_nodes; i++) num[d[i]]++;
    int u = source;
    memset(cur, 0, sizeof(cur));
    while (d[source] < num_nodes) {
        if (u == sink) {
            flow += augment();
            u = source;
        }
        bool advanced = false;
        for (int i = cur[u]; i < G[u].size(); i++) { 
            Edge& e = edges[G[u][i]];
            if (e.capacity > e.flow && d[u] == d[e.to] + 1) {
                advanced = true;
                p[e.to] = G[u][i];
                cur[u] = i;
                u = e.to;
                break;
            }
        }
        if (!advanced) { // retreat
            int m = num_nodes - 1;
            for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix)
                if (edges[*ix].capacity > edges[*ix].flow)
                    m = min(m, d[edges[*ix].to]);
            if (--num[d[u]] == 0) break; // gap 优化
            num[d[u] = m+1]++;
            cur[u] = 0;
            if (u != source)
                u = edges[p[u]].from;
        }
    }
    return flow;
}

第二个代码：

struct Edge{  
    int v, c, n;  
};  
Edge edge[maxE];  
int adj[maxN], cntE;  
int Q[maxQ], head, tail, inq[maxN];  
int d[maxN], num[maxN], cur[maxN], pre[maxN];  
int s, t, nv;  
int n, m, dis[maxN], G[maxN][maxN];  
void addedge(int u, int v, int c){  
    edge[cntE].v = v; edge[cntE].c = c; edge[cntE].n = adj[u]; adj[u] = cntE++;  
    edge[cntE].v = u; edge[cntE].c = 0; edge[cntE].n = adj[v]; adj[v] = cntE++;  
}  
void REV_BFS(){  
    clear(d, -1);  
    clear(num, 0);  
    head = tail = 0;  
    d[t] = 0;  
    num[0] = 1;  
    Q[tail++] = t;  
    while(head != tail){  
        int u = Q[head++];  
        head %= maxQ;  
        for(int i = adj[u]; ~i; i = edge[i].n){  
            int v = edge[i].v;  
            if(~d[v]) continue;  
            d[v] = d[u] + 1;  
            num[d[v]]++;  
            Q[tail++] = v;  
            tail %= maxQ;  
        }  
    }  
}  
int ISAP(){  
    copy(cur, adj);  
    REV_BFS();  
    int flow = 0, u = pre[s] = s, i;  
    while(d[s] < nv){  
        if(u == t){  
            int f = oo, neck;  
            for(i = s; i != t; i = edge[cur[i]].v){  
                if(f > edge[cur[i]].c){  
                    f = edge[cur[i]].c;  
                    neck = i;  
                }  
            }  
            for(i = s; i != t; i = edge[cur[i]].v){  
                edge[cur[i]].c -= f;  
                edge[cur[i] ^ 1].c += f;  
            }  
            flow += f;  
            u = neck;  
        }  
        for(i = cur[u]; ~i; i = edge[i].n) if(edge[i].c && d[u] == d[edge[i].v] + 1) break;  
        if(~i){  
            cur[u] = i;  
            pre[edge[i].v] = u;  
            u = edge[i].v;  
        }  
        else{  
            if(0 == (--num[d[u]])) break;  
            int mind = nv;  
            for(i = adj[u]; ~i; i = edge[i].n){  
                if(edge[i].c && mind > d[edge[i].v]){  
                    mind = d[edge[i].v];  
                    cur[u] = i;  
                }  
            }  
            d[u] = mind + 1;  
            num[d[u]]++;  
            u = pre[u];  
        }  
    }  
    return flow;  
}