##网络流##

HLPP参考：http://blog.sina.com.cn/s/blog_60707c0f01010oqt.html

ISAP参考：http://blog.sina.com.cn/s/blog_6165e6920100ea08.html

参考：http://blog.sina.com.cn/s/blog_476a25110100mag7.html

http://wenku.baidu.com/link?url=Q3y_0x9BaO1L1W5Kr3AJzgCUVJVKNjG8UGw2nt44d6HDWTT3QLWygZ8IT_lb-V6GNjYoV4Fh74bsB88uTMeJA4xWBnbAyNfUr9VRCo-rCFS

//algo1://Ford-Fulkerson//O(n*m*U)#define N 115struct arc{    int c,f;}e[N][N];int flag[N],al[N],pre[N];int n,m,np,nc;int ford(int s,int t,int n){   while (1)        {   memset(flag,255,sizeof(flag));            memset(pre,255,sizeof(pre));            memset(al,255,sizeof(al));            queue<int > q;            flag[s]=0; pre[s]=-1; al[s]=inf;            q.push(s);            while (!q.empty()&&flag[t]==-1)                {                    int v=q.front(); q.pop();                    for (int i=0;i<n;i++)                    if (flag[i]==-1)                        {                            if (e[v][i].f<e[v][i].c)                                {                                    flag[i]=0; pre[i]=v;                                    al[i]=min(al[v],e[v][i].c-e[v][i].f);                                    q.push(i);                                }                            else if (e[i][v].f>0)                                {   flag[i]=0; pre[i]=-v;                                    al[i]=min(al[v],e[i][v].f);                                    q.push(i);                                }                        }                    flag[v]=1;                }            if (flag[t]==-1||al[t]==0) break;            int k1=t,k2=abs(pre[k1]),alph=al[t];            while (1)                {   if (e[k2][k1].c>0)                         e[k2][k1].f+=alph;                    else e[k1][k2].f-=alph;                    if (k2==s) break;                    k1=k2; k2=abs(pre[k1]);                }        }    int ans=0;    for (int i=0;i<n;i++)        ans+=e[s][i].f;    return ans;}void doit(){   for (int i=0;i<N;i++)    for (int j=0;j<N;j++)        e[i][j].c=e[i][j].f=0;    //add edges    printf("%d\n",ford(s,t,N));}

//algo2:  //Dinic//O(n*n*m)struct dinic{    int g[N],d[N],G[N];    int np[M],to[M],cap[M];    int n,s,t,e;    void AddEdge(int x,int y,int z,int z2=0)         {  to[e]=y;cap[e]=z;            np[e]=g[x];g[x]=e++;            to[e]=x;cap[e]=z2;            np[e]=g[y];g[y]=e++;         }    void init(int n)         {  this->n=n;e=0;            memset(g,255,sizeof(g));         }    int maxflow(int s,int t)         {  this->s=s; this->t=t;            int ans=0;            while (bfs())                {for(int i=0;i<n;i++)G[i]=g[i];                 ans+=dfs(s,inf);                }            return ans;         }    bool bfs()//建立层次图    {        memset(d,-1,sizeof(d));        queue<int > q;        q.push(s);        d[s]=0;        while(!q.empty())        {            int x=q.front(); q.pop();            for(int i=g[x];~i;i=np[i])                 {                     int y=to[i];                     if(cap[i]&&d[y]==-1)                     {  d[y]=d[x]+1;                        if(y==t)return 1;                        q.push(y);                     }                 }        }        return 0;    }   int dfs(int x,int sum)    {        if (x==t|| sum==0) return sum;        int ss=sum,ret;        for(int &i=G[x];~i;i=np[i]) //当前弧优化        {            int y=to[i];            if(d[y]==d[x]+1&&(ret=dfs(y,min(sum,cap[i]))))                {                    cap[i]-=ret;                    cap[i^1]+=ret;                    sum-=ret;                    if (!sum) break;                }        }        return ss-sum;    }}it;

SAP的两个优化：

GAP优化——如果一次重标号时，出现距离断层，则可以证明ST无可行流，此时则可以直接退出算法。
当前弧优化——为了使每次找增广路的时间变成均摊O(V)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧。

//algo3_1://SAP//O(n*n*m)struct SAP{    int g[N],d[N],G[N],pre[N],numd[N];    int np[M],to[M],cap[M];    int n,e,s,t;    void AddEdge(int x,int y,int z,int z2=0)        {  to[e]=y;cap[e]=z;           np[e]=g[x];g[x]=e++;           to[e]=x;cap[e]=z2;           np[e]=g[y];g[y]=e++;        }    void init(int n)        {   this->n=n;            e=0;memset(g,255,sizeof(g));        }    void bfs()//先用广度优先算出每个点的d值        {   memset(numd,0,sizeof(numd));            for(int i=0; i<n; i++)                numd[d[i]=n]++;            d[t]=0; numd[n]--; numd[0]++;            queue<int> Q;Q.push(t);            while(!Q.empty())            {   int v=Q.front();Q.pop();                int i=g[v];                while(i!=-1)                    {                        int u=to[i];                        if(d[u]<n){i=np[i];continue;}                        d[u]=d[v]+1;                        numd[n]--;numd[d[u]]++;                        Q.push(u);                        i=np[i];                    }            }        }    int maxflow(int s,int t)    {   this->s=s; this->t=t;        for(int i=0;i<n;i++)            G[i]=g[i];        int max_flow=0;        bfs();        int u=s ;//从源点搜一条到汇点的增广路        while(d[s]<n)//就算所有的点连成一条线源点的d值也是最多是n-1        {            if(u==t)//如果找到一条增广路径            {                    int cur_flow=inf,neck;//找到那条瓶颈边                    for(int from=s; from!=t; from=to[G[from]])                        if(cur_flow>cap[G[from]])                        {   neck=from;cur_flow=cap[G[from]];}                for(int from=s; from!=t; from=to[G[from]]) //修改增广路上的边的容量                    {                    int tmp=G[from];                    cap[tmp]-=cur_flow;                    cap[tmp^1]+=cur_flow;                    }                max_flow+=cur_flow;//累加计算最大流                u=neck;//下一次搜索直接从瓶颈边的前一个节点搜起            }            int i;            for(i=G[u];~i; i=np[i]) //从当前点开始找一条允许弧            if(cap[i]&&d[u]==d[to[i]]+1)//如果找到跳出循环                break;            if(i!=-1)//找到一条允许弧            {                G[u]=i;//从点u出发的允许弧的地址                pre[to[i]]=u;//允许弧上下一个点的前驱为u                u=to[i];//u变成下一个点继续搜直到搜出一条增广路            }            else  //如果没有搜到允许弧            {                numd[d[u]]--; //d[u]将被修改所以numd[d[u]]减一                if(!numd[d[u]]) break;  //如果没有点的d值为d[u]则不可能再搜到增广路结束搜索                G[u]=g[u];  //当前点的允许弧为第一条边                int tmp=n;                for(int j=g[u]; ~j; j=np[j]) //搜与u相连的点中d值最小的                    if(cap[j]&&tmp>d[to[j]])tmp=d[to[j]];                d[u]=tmp+1; //修改d[u]                numd[d[u]]++;                if(u!=s)u=pre[u];//从u的前驱搜，因为从u没有搜到允许弧            }        }    return max_flow;    }}it;

//algo3_2://SAP//O(n*n*m)#define N 105int n,m,np,nc,st,ed;struct SAP{    int g[N][N];//直接在原图上操作了    int d[N],gap[N];    int n,s,t;    void AddEdge(int x,int y,int z)        {            g[x][y]=z;        }    void init(int n)        {   this->n=n;            memset(g,0,sizeof(g));            memset(gap,0,sizeof(gap));            memset(d,0,sizeof(d));        }    int maxflow(int s,int t)        {  this->s=s;this->t=t;           int ans=0;           for(gap[0]=n;d[s]<n;)                ans+=sap(s,inf);           return ans;        }    int sap(int u,int flow)        {            if(u==t) return flow;            int res=flow;            for(int i=0;i<n;++i)                if(g[u][i] && d[u]==d[i]+1)                {                    int tmp=sap(i,min(res,g[u][i]));                    g[u][i]-=tmp,g[i][u]+=tmp;                    if(!(res-=tmp))return flow;                }            if(!--gap[d[u]])d[s]=n;            ++gap[++d[u]];            return flow-res;        }}it;

//algo4:  //ISAP  //O(n*n*m)  //ISAP 的一个所谓 GAP 优化:  //由于从 s 到 t 的一条最短路径的顶点距离标号单调递减，且相邻顶点标号差严格等于1，  //因此可以预见如果在当前网络中距离标号为 k (0 <= k < n) 的顶点数为 0，  //那么可以知道一定不存在一条从 s 到 t 的增广路径，此时可直接跳出主循环。  struct ISAP  {   struct Edge      {          int from,to,cap,flow;          Edge(){}          Edge(int a,int b,int c,int d):from(a),to(b),cap(c),flow(d){}      };      int n,m,s,t;//结点数，边数(含反向弧)，源点，汇点      vector<Edge> edges;//边表，edges[e]&edges[e^1]互为反向弧      vector<int> G[N];//邻接表，G[i][j]表示结点i的第j条边在e数组中的序号      bool vis[N];//BFS使用      int d[N];//从起点到i的距离      int cur[N];//当前弧下标      int p[N];//可增广路上的上一条弧      int num[N];//距离标号计数        void AddEdge(int from,int to,int cap)//重边不影响      {          edges.push_back(Edge(from,to,cap,0));          edges.push_back(Edge(to,from,0,0));//容量为0,表示反向弧          m=edges.size();          G[from].push_back(m-2);          G[to].push_back(m-1);      }        void init(int n)      {          this->n=n;          for(int i=0;i<n;++i) G[i].clear();          edges.clear();      }        void BFS()//反向      {   for (int i=0;i<n;i++)            d[i]=n+10; //!!!        memset(vis,0,sizeof(vis));        queue<int> Q;        Q.push(t);        d[t]=0;        vis[t]=1;        while(!Q.empty())        {            int x=Q.front();            Q.pop();            for(int i=0; i<G[x].size(); ++i)            {                Edge& e=edges[G[x][i]^1];                if(!vis[e.from]&&e.cap>e.flow)                {                    vis[e.from]=1;                    d[e.from]=d[x]+1;                    Q.push(e.from);                }            }        }    }        int Augment()      {          int x=t,a=INF;          while(x!=s)          {              Edge& e=edges[p[x]];              a=min(a,e.cap-e.flow);              x=edges[p[x]].from;          }          x=t;          while(x!=s)          {              edges[p[x]].flow+=a;              edges[p[x]^1].flow-=a;              x=edges[p[x]].from;          }          return a;      }        int maxflow(int s,int t)//结点数      {          this->s=s,this->t=t;          int flow=0;          BFS();          memset(num,0,sizeof(num));          for(int i=0;i<n;++i) if (d[i]!=n+10)++num[d[i]];//!!!         int x=s;          memset(cur,0,sizeof(cur));          while(d[s]<n)          {              if(x==t)              {                  flow+=Augment();                  x=s;              }              int ok=0;              for(int i=cur[x];i<G[x].size();++i)              {                  Edge& e=edges[G[x][i]];                  if(e.cap>e.flow&&d[x]==d[e.to]+1)//Advance                  {                      ok=1;                      p[e.to]=G[x][i];                      cur[x]=i;                      x=e.to;                      break;                  }              }              if(!ok)//Retreat              {                  int m=n-1;                  for(int i=0;i<G[x].size();++i)                  {                      Edge& e=edges[G[x][i]];                      if(e.cap>e.flow) m=min(m,d[e.to]);                  }                  if(--num[d[x]]==0) break;//gap优化                  num[d[x]=m+1]++;                  cur[x]=0;                  if(x!=s) x=edges[p[x]].from;              }          }          return flow;      }  }it;  

algorithm Improved-Shortest-Augmenting-Path  f <-- 0  从终点 t 开始进行一遍反向 BFS 求得所有顶点的起始距离标号 d(i)  i <-- s  while d(s) < n do      if i = t then    // 找到增广路          Augment          i <-- s      // 从源点 s 开始下次寻找      if Gf 包含从 i 出发的一条允许弧 (i,j)          then Advance(i)         else Retreat(i)    // 没有从 i 出发的允许弧则回退  return fprocedure Advance(i)  设 (i,j) 为从 i 出发的一条允许弧  pi(j) <-- i    // 保存一条反向路径，为回退时准备  i <-- j        // 前进一步，使 j 成为当前结点procedure Retreat(i)  d(i) <-- 1 + min{d(j):(i,j)属于残量网络Gf}    // 称为重标号的操作  if i != s     then i <-- pi(i)    // 回退一步procedure Augment  pi 中记录为当前找到的增广路 P  delta <-- min{rij:(i,j)属于P}  沿路径 P 增广 delta 的流量  更新残量网络 Gf

~~~~~~~~~~~~~~~~~~~~~~~~~~~增广路与预流推进的分割线~~~~~~~~~~~~~~~~~~~~~~~~~~~~

定义e(u)为节点u上的余流，实际上我们并不需要计算各个h(v)，只需将h(s)标为n，其余标为0，算法也能逐步计算出各个点的h。网络流中的节点(除s, t)仅需要满足流入量不小于流出量。其中流入量大于流出量的节点，我们称之为活动节点。预流推进算法就是不断地将活动结点，变为非活动结点，使得预流成为可行流。具体分为两个方法：
①压入操作，选出一个活动顶点u。并依次判断残余网络G中每条边(u, v)，若h(u) = h(v) + 1，则顺着这里条边推流，流量变化量df(u,v)=min(e[u],cf(u,v))，同时更新相应的e(u)，e(v)；
②重标记操作，选出一个活动结点。如果对于残余网络G中每条边(u, v)，都有h(v)大于或等于h(u)，则需要对u进行重新标号，h(u) = min{h(v) + 1}；
具体实现又可以根据对活跃点的推进过程分为：
①一般的预流推进算法（Generic Preflow-Push Algorithm），在残余网络中，维护一个预流，不断对活跃点执行push操作，或者relable操作来重新调整这个预流，直到不能操作。时间复杂度为O(V*E*E)；
②先进先出预流推进算法，在残余网络中，以先进先出队列维护活跃点，时间复杂度为O(V*V*V)；
③最高标号预流推进算法（Highest Preflow-Push Algorithm），在残余网络图中，每次检查最高标号的活跃点，需要用到优先队列，时间复杂度为O(V*V*E^0.5)；

//algo5://HLPP//O(n*n*m^(1/2))#define N 105int d[N];struct HLPP{    int e[N];    int g[N][N];    struct node{         int index;         friend bool operator < (const node& a, const node& b)            {return d[a.index] > d[b.index];}    };    int n,s,t;    void Relabel(int u)        {            int m1=inf;            for(int i=0; i<n; ++i)                if(g[u][i] && d[i] < m1)m1 =d[i];            d[u]=m1+1;        }    void AddEdge(int x,int y,int z)        {            g[x][y]+=z;        }    void init(int n)        {   this->n=n;            memset(g,0,sizeof(g));        }    int maxflow(int s,int t)        {  node no;           priority_queue<node> Q;           memset(e,0,sizeof(e));           memset(d,0,sizeof(d));           d[s]=n;           for(int i=0;i<n;++i)           if(g[s][i])                {                    int df=g[s][i];                    e[s] -= df;                    e[i] += df;                    g[s][i] -= df;                    g[i][s] += df;                    no.index = i;                    if(i!=s&& i!=t) Q.push(no);                }           while(!Q.empty())                {                 no=Q.top();Q.pop();                 int u=no.index;                 while(e[u])                     {  int i;                        for(i = 0; i < n; ++i)                            if(g[u][i] && d[u] > d[i])break;                        if(i == n) Relabel(u);                        for(int i=0; i<n; ++i)                             if(d[u]==d[i] + 1 && g[u][i])                             {                                 int df = min(e[u], g[u][i]);                                 g[u][i] -= df;                                 g[i][u] += df;                                 e[u] -= df;                                 e[i] += df;;                                 no.index = i;                                 if(e[i] > 0&&i!=s &&i!=t) Q.push(no);                             }                     }                }            return e[t];        }}it;

~~~~~~~~~~~~~~~~~~~~~~~~~~~最大流与最小费用最大流间的分割线~~~~~~~~~~~~~~~~~~~~~~~~~~~~

struct MCMF{   struct Edge    {        int from,to,cap,flow,cost;        Edge(){}        Edge(int a,int b,int c,int d,int e):from(a),to(b),cap(c),flow(d),cost(e){}    };    int n,m,s,t;    vector<Edge> edges;    vector<int> G[N];    int inq[N],d[N],p[N],a[N];    queue<int> Q;    void init(int n)    {        this->n=n;        for(int i=0;i<n;++i) G[i].clear();        edges.clear();    }    void AddEdge(int from,int to,int cap,int cost)    {   //printf("%d %d %d %d\n",from,to,cap,cost);        edges.push_back(Edge(from,to,cap,0,cost));        edges.push_back(Edge(to,from,0,0,-cost));        m=edges.size();        G[from].push_back(m-2);        G[to].push_back(m-1);    }    bool SPFA(int s,int t,int & flow,int & cost)    {        memset(d,0x3f,sizeof(d));        memset(inq,0,sizeof(inq));        d[s]=0;inq[s]=1;p[s]=0;a[s]=inf;        Q.push(s);        while(!Q.empty())        {            int u=Q.front();Q.pop();            inq[u]=0;            for(int i=0;i<G[u].size();++i)            {                Edge& e=edges[G[u][i]];                if(e.cap>e.flow&&d[e.to]>d[u]+e.cost)                {                    d[e.to]=d[u]+e.cost;                    p[e.to]=G[u][i];                    a[e.to]=min(a[u],e.cap-e.flow);                    if(!inq[e.to]) {Q.push(e.to);inq[e.to]=1;}                }            }        }        if(d[t]==inf) return false;        flow+=a[t];//固定流量k，检查，若flow+a>=k，只增加k-flow个单位        cost+=d[t]*a[t];        int u=t;        while(u!=s)        {            edges[p[u]].flow+=a[t];            edges[p[u]^1].flow-=a[t];            u=edges[p[u]].from;        }        return true;    }    int mincost(int s,int t)    {        int flow=0,cost=0;        while(SPFA(s,t,flow,cost));        return cost;    }}it;