二分图相关与KM算法
前言
呜，好久没写博客了，DDL 也有好多，一不留神就轮到我了呜。
看了一眼其它同学写的博客，什么数模啊，什么 CTF 啊，什么 Python 爬虫啊，感觉自己真是越来越菜了呜。
然后在我一愁莫展之际，我突然发现，我去年开学初还有一篇写了一半的博客没有发。
于是我就把它补全重新排了下版然后扔上去敷衍部长了
本帖为蓝书中 5.5二分图匹配 相关内容的自学笔记
对书中大多数难懂的内容如最佳完美匹配做了详细的解释并加上了自己的理解
（毕竟原文是真的谜语人）
可能大概似乎可以为图论 网络流的学习作参考
还含有本人尝试使用dfs实现n^3的心路历程
、
很菜，轻喷。如果你发现下面有的地方没有用markdown写公式 那是因为我实在懒得改了。
（加粗看心情，加粗的东西不一定有用，只是方便阅读）
一.何为二分图
用黑白两种颜色给图上所有节点染色 若能在染色后使每一条边两端的点颜色不同 则为二分图
二.何为匹配
选出边 使取出的任意两条边间无公共点 则选出的边构成图的一个匹配
我们将按上述原则取出来的边集中的边称为匹配边 每条匹配边两端的点称为匹配点
特别的：
在二分图中 能选出的边数最大的边集 称为二分图最大匹配
当所有点都是匹配点时 我们称这个匹配为完美匹配
显然 完美匹配一定是最大匹配
三.增广路与最大匹配
从某一个未被匹配的点出发
走一条非匹配边到某个匹配点 然后再走匹配边走到另一个匹配点 再走非匹配边走到第三个匹配点 重复到最后一个匹配点后 若还能沿非匹配边走到另一个非匹配点
则此时，若将路上经过的所有非匹配边变为匹配边、匹配边变为非匹配边 形成的匹配会比原先的匹配多一条边
我们称找到的这一条非匹配边、匹配边、非匹配边...匹配边、非匹配边的交替路为增广路
当无法找到下一条增广路时 已达到最大匹配
四.从无向图上的移子删点游戏理解完美匹配
Alice在无向联通图上选一个点 然后放下棋子
接下来Bob与Alice轮流移动棋子 并在移动后删去原先站着的点
若轮到某人时已经无法移动 则这个人输掉游戏
在这个游戏中Alice与Bob都使绝对聪明的
我们可以发现
若图存在完美匹配
Alice放下棋子的点必然是一个匹配点 Bob必然能沿着这个匹配点对应的匹配边移动
而在这之后 由于在匹配中 每个匹配点连接的匹配边只有一条 所以Alice只能沿一条非匹配边前往下一个点
而由于完美匹配中图上所有的点都是匹配点 所以这个点又是一个匹配点 又有一条对应的匹配边...
当Bob走完最后一条匹配边后 此时场上只剩下一个点 因而Alice无法移动棋子 会输掉比赛
若图不存在完美匹配
只要Alice只要先求出最大匹配 然后将棋子放在一个非匹配点上
若此非匹配点还连接有另一个非匹配点 那么另一个非匹配点一定可以连到这个匹配中 从而使这个匹配更大 也就是说会存在更大的匹配
但是Alice求出的已经是最大匹配 所以Bob的移动必然往一个匹配点上走
此后Alice只要不断沿着匹配边走就行了
只要Bob能动 Bob必然是沿一条非匹配边到另一个匹配点
而之后Alice就总能沿着匹配边往下走
因为如果Bob可以找到一条非匹配边连接到另一个非匹配点
那么从这个多出来的非匹配点连回某个匹配点 就可以形成一条增广路(毕竟至少比原有的匹配对应的交替路多了两条边)
这会使匹配变得更大(边数更多)
然而Alice求出的已经是最大匹配了
顺带一提 无向联通图中若无法构成完美匹配 最大匹配外余出来的点最多只有一个 否则必然存在增广路使这个匹配更大
五.如何寻找一条增广路
上面介绍增广路时讲了 选择一个未匹配点 不断走 尝试找到下一个非匹配点 找到下一个非匹配点时 便获得了一条增广路
具体实现以luogu P3386为例子
这道题直接分出了二分图的两个点集 所以省去了染色找点的步骤
当通过 dfs 找到某条增广路后 该点已经匹配 没有必要再搜一次了 (找增广路时首尾是非匹配点)
如果某个点此时搜不到增广路 那么之后也不可能搜到从这个点出发的增广路 也就是说每个点只需要搜索一次
每个点搜索增广路时 dfs 过一次的点都没有再跑一遍的必要 因为这个点出发的所有路径都已经被搜索过了
就算之后再从其它路径抵达 既然之前没有找到 这次也不可能找到
这就是单次搜索中 vis 标记在回溯时不需要清0的原因
至于为什么 vis 标记打在另一个点集里 是因为 dfs 找增广路时 我们会让它朝着匹配路前进
所以对于每个 vis 过的点 实际上下一次扫到会再跑的路径也是一样的 即沿增广路 所以不需要再跑一次
当然 vis 打在前一个点集里也是可以的 细节上有小变动 就是判断的不是vis[i]而是vis[to[i]]
#include <cstdio>
#include <cstring>
#define REG register
#define rep(i,a,b) for(REG int i=a;i<=b;i++)
#define Rep(i,a,b) for(REG int i=a;i>=b;i--)
inline char getc(){
static char buf[1<<14],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<14,stdin),p1==p2)?EOF:*p1++;
}
inline int scan(){
REG int x=0; REG char ch=0;
while(ch<48) ch=getc();
while(ch>=48) x=x*10+ch-48,ch=getc();
return x;
}
const int N=502;
int n,m,c,to[N],link[N],res; bool e[N][N],vis[N];
inline bool dfs(int u){
rep(i,1,m){
if(!e[u][i]||vis[i]) continue;
vis[i]=1;
if(!to[i]||dfs(to[i])){
to[i]=u; link[u]=i;
return 1;
}
}
return 0;
}
int main(){
n=scan(); m=scan(); c=scan();
rep(i,1,c){
REG int u,v;
u=scan(); v=scan();
e[u][v]=1;
if(!link[u]&&!to[v]) res++,to[v]=u,link[u]=v;
}
rep(i,1,n){
if(link[i]) continue;
memset(vis,0,sizeof(vis));
if(dfs(i)) res++;
}
printf("%d",res);
return 0;
}
六.何为匈牙利树
匈牙利树就是：用 dfs 找增广路且没能找到时 遍历完所有的点后 走过的边所形成的树
简单的说 就是沿交替路遍历出来的一棵 dfs树
我们记 dfs 起点也就是未匹配点所在点集为A 另一个点集则记为B
由于我们走的是交替路 我们对 dfs 做出如下要求
当 dfs 进行到B点集时 必须沿着匹配边往下走
(当然你也可以把A和B倒过来 这里为了下文讨论方便如此规定)
七.二分图最佳完美匹配
就是在完全二分图中求出边权和最大的完美匹配
(个人觉得 这个理解起来比较麻烦 大概是我比较菜吧)
这里书中引入了一个函数 \(l(u)\) 我们将它叫做顶标
这个函数的值是人为规定的 只要这个值满足在图中总是有 \(l(u)+l(v)>=edge(u,v)\) 就行
比如说我们把任意点 \(x\) 的 \(l(x)\) 的值全部定为 \(+\infty\) 那么上述条件同样满足
只要满足条件 我们就可以称 \(l(x)\) 为可行顶标
而如果我们将 \(l(x)\) 的值确定后 却发现某条边两端点的顶标之和比这条边的边权还小 此时 这个顶标是不可行的
接下来书中介绍了一个极为重要的定理 ：
从这个完全二分图G中取出一定数量的边生成一个子图M 若子图含有所有点 且总有 \(l(u)+r(v)=edge(u,v)\) 我们称这个子图为相等子图
那么 若此子图有完全匹配 则此时的完美匹配便是最佳匹配
证明书里很清楚：
记相等子图M的完美匹配为M' 则显然M'的边权和为 \(l(x)\) 的总和
而我们从二分图中取出一个完美匹配 由于总有 \(l(u)+l(v)>=edge(u,v)\) 我们求出的任何一个完美匹配的边权和总会小于或等于 \(l(x)\) 的总和
因此 利用相等子图求出的完美匹配一定是最佳完美匹配
那么问题就变为 如何调整顶标的值 使相等子图有完美匹配
（下面是我边看边想边写整出来的）
先来理论一波：
现在假设我们以某种方式规定好了各个点的 \(l(x)\)
在 dfs 时 我们无视 \(l(u)+l(v)>w(u,v)\) 的边 只走 \(l(u)+l(v)=w(u,v)\) 的
对每个点都 dfs 后 我们扫出了一个基于相等子图的最大匹配
之后再从相等子图中某个非匹配点进行 dfs 可以找到一棵基于相等子图的匈牙利树
(这里之所以说一棵 是因为相等子图并不一定是连通的)
我们记这棵树上分别含于二分图两个点集的点集为 \(S\) 和 \(T\) 而二分图两个点集中不在树上的点的集合记为 \(S'\) 和 \(T'\)
(即 \(S,S'\subset A\)，\(T,T'\subset B\) )
P.S.我在下面的描述中加上了"合法"两个字 意思是指这条边满足 \(l(u)+l(v)=w(u,v)\) 的要求
书中给出的性质如下：
-
\(S\) 中的点到 \(T'\) 没有连合法边
我们知道，在寻找增广路也就是 dfs 的过程中，我们会把扫到的 \(A\) 点集中的点的所有连到 \(B\) 点集的边都走一遍。
(这一点在上面找最大匹配的代码中有体现)
因此，如果 \(S\) 到 \(T'\) 有连边，dfs 过程中一定会被扫到，那么这个点就会被归到T。
(这就是书中证明部分所谓的"匈牙利树可以继续生长") -
\(S'\) 到 \(T\) 的所有合法边都是非匹配边
这个很明显，如果是匹配边 ，dfs 的时候必然会沿着匹配边跑到 \(S'\) 这个点就必然在匈牙利树上，会被归到 \(S\) 集合里
接下来我们可以试着执行某种操作，使这个子图仍然是相等子图
我们不难发现 对于S中点 \(l(x)\) 同时加上某个值 再对T中的对应的 \(l(y)\) 同时减去相同的值 \(l(x)+l(y)=w(x,y)\) 依旧成立
我们知道，由于性质一，所以 \(S\) 到 \(T'\) 的边总满足 \(l(u)+l(v)>w(u,v)\) 那么如果我们将 \(S\) 的点顶标降低
因为 \(l(x)\) 变小了 就有可能使某条边变为为 \(l(x)+l(y)=w(x,y)\) 进而将这条边加入相等子图
下面我们记 \(x\in S\)，\(y\in T'\)。
原本有 \(l(x)+l(y)>w(x,y)\)
则只要让改变后的 \(l'(x)= l(x) -( l(x) +l(y)-w(x,y))\)
便能使边(x,y)进入相等子图
但是由于我们要维持 \(S\) 到 \(T\) 的原有边合法（即要维持原有的顶标相加等于边权的关系），所以 \(S\) 中所有点的顶标必须同时减去一个数 \(a\)（同时 \(T\) 中的点顶标相应加上 \(a\)），而不能为了使某条特定边进入而单独改动某个点。
我们可以看出
\(a=min\{ l(x) +l(y)-w(x,y)|x \in S,y \in T'\}\)
如果减的比这个大 会出现 \(l(x)+l(y)<w(x,y)\) 的边 这不满足顶标的定义 你的顶标会变得不可行
如果减的比这个小 一条边也进不去
当有新边进入时
这条新边的两端中x点在原有相等子图中
(因为右端点的顶标增加了 不会有新边从那边进来的)
- 若y是非匹配点
我们必然能找到一条增广路：
若x也是非匹配点 则形成一条只有两个点的增广路
若x是匹配点 则形成一条有匈牙利树起点(某个非匹配点)经过x到y的增广路 - 若y是匹配点 意味着我们连通了相等子图两个匹配 (或是说联通了两个相等子图)
那么 我们可以直接将这个y点加入T集合 把y点匹配的点z加入S集合
这条匹配边(z,y)与非匹配边(x,y)则并入这棵匈牙利树中
总之 这一波操作后 我们得到两种结果:
- 找到增广路
- 一条原本独立的匹配边经由一条新加入的非匹配边进入匈牙利树
(就是树上多了两条边 但整个合法二分图的匹配数并没有增加 因为匹配是本来就有的)
这样的操作重复几次后 我们最终会找到一个含有所有点的最佳完美匹配
(前提是有解)
接下来我们考虑实现
(以luogu P6577 【模板】二分图最大权完美匹配 为例)
朴素算法
首先 我们把每个点的顶标初始化为其所连接的边中最大的边权 显然 这样的顶标是可行的
接下来我们枚举每一个A点集中的点
对每个点不断进行 dfs 直到找到相等子图上以之为起点的增广路
根据前面讲增广路算法时的原理 找到增广路后此点已匹配 便不再 dfs 了 这时我们再去 dfs 下一个点
每次 dfs 中
我们将扫到的A点集中的点归入S B点集中的点归入T
每次 dfs 时 如果未找到增广路（相等子图中的这条增广路已无法再扩展） 我们会扩展出一棵匈牙利树
接下来我们枚举S集合中的点 算出a也就是这个需要同时减去/增加的值 然后同时更改S与T集合的顶标
如果找到增广路 我们就找下一个点...
下面的代码，我是自己按着理解打的。看起来跟书中的有很大的区别。
不过本质思想上是一样的，只不过是我在存边啊之类的细节上使用了不同的手段
复杂度可以简单的这样理解：
每次修改顶标 我们需要枚举两个点集的点来计算a 总共是O(n^2)
每次修改后 至少有一条边进入相等子图
我们最多需要把\(n^2\)条边(或者说所有的m条边)都加入进去 总复杂度就是\(O(n^4)\)或是\(O(n^2*m)\)
然而事实上\(n^4\)是不可能出现的 其本身也不可能卡到满 朴素算法在随机数据下可以有优秀的表现
#include <cstdio>
#include <cstring>
typedef long long ll;
#define REG register
#define rep(i,a,b) for(REG int i=a;i<=b;i++)
#define Rep(i,a,b) for(REG int i=a;i>=b;i--)
inline char getc(){
static char buf[1<<14],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<14,stdin),p1==p2)?EOF:*p1++;
}
inline int scan(){
REG int x=0,f=1; REG char ch=0;
while(ch<48) ch=='-'&&(f=-1),ch=getc();
while(ch>=48) x=x*10+ch-48,ch=getc();
return x*f;
}
template<typename T> inline void print(const T&x){
if(x>9) print(x/10),putchar(x%10+48);
else putchar(x+48);
}
template<typename T> inline T max(const T&x,const T&y){
return x<y?y:x;
}
template<typename T> inline T min(const T&x,const T&y){
return x<y?x:y;
}
const int N=503,M=3e5+1,INF=1e9;
int n,m,lA[N],lB[N],left[N]; bool link[N];
int head[N],nex[M],to[M],wei[M];
int S[N],cs; bool vis[N];
//e为边的标号 w为边权 l为A/B顶标 left为B点集的匹配点 link为A点集是否已匹配
//邻接表 事实证明我们只需要连条A到B的有向边就行了 反正dfs的规则就在那 你又不会往回扫
//S用于存储进入集合的各个点 vis用于标记B集合中的点是否已被扫过
inline void addEdge(int u,int v,int w){
static int cc=0;
nex[++cc]=head[u]; to[cc]=v; wei[cc]=w; head[u]=cc;
}
inline bool dfs(int u){
S[++cs]=u;
for(REG int i=head[u];i;i=nex[i]){
REG int &v=to[i];
if(lA[u]+lB[v]!=wei[i]||vis[v]) continue;
vis[v]=1;
if(!left[v]||dfs(left[v])){
left[v]=u; link[u]=1;
return 1;
}
}
return 0;
}
inline void update_lx(){
REG int d=INF;
rep(p,1,cs){
REG int &u=S[p];
for(REG int i=head[u];i;i=nex[i]){
if(vis[to[i]]) continue;
d=min(d,lA[u]+lB[to[i]]-wei[i]);
}
}
rep(p,1,cs) lA[S[p]]-=d;
rep(i,1,n) if(vis[i]) lB[i]+=d;
}
int main(){
n=scan(); m=scan();
rep(i,1,m){
REG int u=scan(),v=scan(),w=scan();
addEdge(u,v,w);
lA[u]=max(lA[u],w); lB[v]=max(lB[v],w);
}
Rep(i,n,1){
if(link[i]) continue;
while(1){
cs=0; memset(vis,0,sizeof(vis));
if(dfs(i)) break;
update_lx();
}
}
REG ll ans=0;
rep(i,1,n) ans+=lA[i]+lB[i];
print(ans); putchar('\n');
rep(i,1,n) print(left[i]),putchar(' ');
return 0;
}
\(O(n^3)\) 算法尝试1
我们考虑降低顶标修改时计算 \(a\) 的复杂度
我们可以考虑某种预处理 提前算出每个点y对应的 \(l(x)+l(y)-w(x,y)\) 的最小值 \(slack(y)\)
然后我试着在 dfs 的时候顺带跑一跑，然后T飞了
其实我这么写依然是 \(O(n^4)\) 算法
因为我无法随着 dfs 同时进行 slack 的更改，必须另开一个循环来计算，否则会漏掉点
但是这样一来，最坏情况下，我就不得不遍历所有的边，一样会被卡成 \(n^4\)
只不过被卡的位置从顶标修改操作变成了 dfs 操作
#include <cstdio>
#include <cstring>
typedef long long ll;
#define REG register
#define rep(i,a,b) for(REG int i=a;i<=b;i++)
#define Rep(i,a,b) for(REG int i=a;i>=b;i--)
inline char getc(){
static char buf[1<<14],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<14,stdin),p1==p2)?EOF:*p1++;
}
inline int scan(){
REG int x=0,f=1; REG char ch=0;
while(ch<48) ch=='-'&&(f=-1),ch=getc();
while(ch>=48) x=x*10+ch-48,ch=getc();
return x*f;
}
template<typename T> inline void print(const T&x){
if(x>9) print(x/10),putchar(x%10+48);
else putchar(x+48);
}
template<typename T> inline T max(const T&x,const T&y){
return x<y?y:x;
}
template<typename T> inline T min(const T&x,const T&y){
return x<y?x:y;
}
const int N=503,M=3e5+1,INF=1e9;
int n,m,lA[N],lB[N],left[N]; bool link[N];
int head[N],nex[M],to[M],wei[M];
int S[N],cs;
int slack[N]; bool vis[N];
//slack用于存储T'中每个y对应的min(l(x)+l(y)-w(x,y))
inline void addEdge(int u,int v,int w){
static int cc=0;
nex[++cc]=head[u]; to[cc]=v; wei[cc]=w; head[u]=cc;
}
inline bool dfs(int u){
S[++cs]=u;
for(REG int i=head[u];i;i=nex[i]){
REG int &v=to[i];
slack[v]=min(slack[v],lA[u]+lB[v]-wei[i]);
}
for(REG int i=head[u];i;i=nex[i]){
REG int &v=to[i];
if(lA[u]+lB[v]!=wei[i]||vis[v]) continue;
vis[v]=1;
if(!left[v]||dfs(left[v])){
left[v]=u; link[u]=1;
return 1;
}
}
return 0;
}
inline void update_lx(){
REG int d=INF;
// rep(p,1,cs){
// REG int &u=S[p];
// for(REG int i=head[u];i;i=nex[i]){
// if(vis[to[i]]) continue;
// d=min(d,lA[u]+lB[to[i]]-wei[i]);
// }
// }
rep(i,1,n) if(!vis[i]) d=min(d,slack[i]);
rep(p,1,cs) lA[S[p]]-=d;
rep(i,1,n) if(vis[i]) lB[i]+=d;
// printf("now:");
// rep(i,1,n) printf("%d ",slack[i]);
// putchar('\n');
// printf("SS:%d\n",d);
}
int main(){
n=scan(); m=scan();
rep(i,1,m){
REG int u=scan(),v=scan(),w=scan();
addEdge(u,v,w);
lA[u]=max(lA[u],w); lB[v]=max(lB[v],w);
}
Rep(i,n,1){
if(link[i]) continue;
// static int cc; cc=14;
while(1){
cs=0;
memset(vis,0,sizeof(vis));
rep(i,1,n) slack[i]=INF;
if(dfs(i)) break;
update_lx();
// if(!(--cc)) break;
}
}
REG ll ans=0;
rep(i,1,n) ans+=lA[i]+lB[i];
print(ans); putchar('\n');
rep(i,1,n) print(left[i]),putchar(' ');
return 0;
}
继续尝试 \(n^3\) 算法
虽然有点小问题 但是我们打开书 发现至少原理是没有错的
只不过 书中的描述是在完成每一次增广后初始化slack
再利用匈牙利树扩展过程中S或T中顶标的变化是相同的这个性质 实现对slack的 \(O(n)\) 修改
而不是在 dfs 的时候改 毕竟 dfs 过程中并没有什么更改顶标的操作
果然我还是太菜了orz
然而我仔细想了想 在 dfs 过程中改其实没有问题
刚才打炸了主要是我没有考虑到这个slack可以 \(O(n)\) 修改
既然书中没有具体实现的代码 我就照着书上的文字自己摸索吧
首先是slack的初始化
如果找到增广路 则不需要修改顶标 那就直接pass
如果没找到增广路 直接按原思路 边 dfs 边修改顶标就行了 在没有增广路的情况下是会遍历完所有的边的
其实到这里就已经奠定了本次尝试爆炸的必然，只是我实在太菜了没看出来。刚才我们另开了一个循环修改顶标，这次我们没有开，而是顺着 dfs 走 ，但这样仍存在需要遍历所有边的恶劣情况，就算撇开修改顶标的操作，仅仅只是dfs就已经T飞了
无论如何 我们先研究** slack 的更改**
我们知道 slack 算的是 \(S\) 到 \(T'\)，如果不在这里头的点 slack 是不能算的
所以我们需要判断，\(T'\) 中的某个点是否与 \(S\) 联通
这里有一点很显然，就是已经计算出了slack的点显然与 \(S\) 联通
考虑匈牙利树的扩展
要么找到增广路 要么并入一条边
当有新的点并入S时 我们需要计算由这个点引出的slack值
这其实直接在 dfs 中修改就行了 不会出现边被跳过的情况
因为会跳过边的原因是找到了增广路 找到增广路的话直接下一轮再处理就行了
然后又T飞了 原因上文有提到
#include <cstdio>
#include <cstring>
typedef long long ll;
#define REG register
#define rep(i,a,b) for(REG int i=a;i<=b;i++)
#define Rep(i,a,b) for(REG int i=a;i>=b;i--)
inline char getc(){
static char buf[1<<14],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<14,stdin),p1==p2)?EOF:*p1++;
}
inline int scan(){
REG int x=0,f=1; REG char ch=0;
while(ch<48) ch=='-'&&(f=-1),ch=getc();
while(ch>=48) x=x*10+ch-48,ch=getc();
return x*f;
}
template<typename T> inline void print(const T&x){
if(x>9) print(x/10),putchar(x%10+48);
else putchar(x+48);
}
template<typename T> inline T max(const T&x,const T&y){
return x<y?y:x;
}
template<typename T> inline T min(const T&x,const T&y){
return x<y?x:y;
}
const int N=503,M=3e5+1,INF=1e9;
int n,m,lA[N],lB[N],left[N]; bool link[N];
int head[N],nex[M],to[M],wei[M];
int S[N],cs; bool vis[N];
int slack[N];
inline void addEdge(int u,int v,int w){
static int cc=0;
nex[++cc]=head[u]; to[cc]=v; wei[cc]=w; head[u]=cc;
}
inline bool dfs(int u){
S[++cs]=u;
for(REG int i=head[u];i;i=nex[i]){
REG int &v=to[i];
slack[v]=min(slack[v],lA[u]+lB[v]-wei[i]);
if(lA[u]+lB[v]!=wei[i]||vis[v]) continue;
vis[v]=1;
if(!left[v]||dfs(left[v])){
left[v]=u; link[u]=1;
return 1;
}
}
return 0;
}
inline void update_lx(){
REG int d=INF;
rep(i,1,n) if(!vis[i]) d=min(d,slack[i]);
rep(p,1,cs) lA[S[p]]-=d;
rep(i,1,n) if(vis[i]) lB[i]+=d;
rep(i,1,n) if(!vis[i]&&slack[i]!=INF) slack[i]-=d;
}
int main(){
n=scan(); m=scan();
rep(i,1,m){
REG int u=scan(),v=scan(),w=scan();
addEdge(u,v,w);
lA[u]=max(lA[u],w); lB[v]=max(lB[v],w);
}
Rep(i,n,1){
if(link[i]) continue;
rep(j,1,n) slack[j]=INF;
while(1){
cs=0; memset(vis,0,sizeof(vis));
if(dfs(i)) break;
update_lx();
}
}
REG ll ans=0;
rep(i,1,n) ans+=lA[i]+lB[i];
print(ans); putchar('\n');
rep(i,1,n) print(left[i]),putchar(' ');
return 0;
}
这次我们从 dfs 本身的性质想想？
(逐渐暴躁.jpg)
书中只讲到了如何降低顶标修改的复杂度
但是我们需要找n次增广路 每次寻找最多要 dfs 个n次 而每次 dfs 的最坏复杂度是n^2
这方面的优化书中只字不提
那我就从 dfs 过程本身来考虑吧
我们可以发现 在寻找同一条增广路时
多次的 dfs 中总会有一段是在做重复的工作
也就是我们每次 dfs 前都清空子图 然后又把整个子图再跑出来一次 只为了加一条边
这就很傻逼
我们希望之前 dfs 出来的匈牙利树可以保留
扩展时直接从S集合中的点往下扩展
这就需要保留住在匈牙利树最后一层的所有节点和根节点
这时我们发现 由于 dfs 是基于栈的性质
我们保存的并不是一层节点 而是一条从最底层开往根节点的路径
(同一深度但路径不同的点会在回溯时从栈中离开)
所以 要执行这个优化 我们只能使用 bfs
类似的 对slack的操作我们边 bfs 边做修改
等等 仔细想想 就算真的存下S集合中的所有点加入队列
复杂度并不会少 我们照样还要把所有的边再扫一遍
或许我理解错了 bfs 复杂度变低的原理 或许其原理在于继承之前增广时扫出的子图 而不需要在扫描下一个点时重来
（已疯）
直接粗暴地上数据结构好了（？
我们知道 实际上每次加入的新边 都是从S出发的slack最小的边
考虑对每个点预处理出其每条边的slack 然后建立n棵splay
每次尝试寻找增广路时 我们取出当前S集合中最小slack的边
然后给这棵splay打上标记 整棵树减去slack
若这边连到一个匹配点i 就把点left[j]对应的splay与当前点的splay合并
可以借助并查集来合并
我们需要找n次增广路 每次最多合并n棵树 每次合并、取边、删边的最高复杂度是 \(O(\log n^2)\)
可知 复杂度大概为 \(O(n^2\log n^2)\)
接下考虑如何合并splay
首先把两个的最小值旋到root 然后比较二者的大小
？？？？
我在说什么鬼话？？
代码如下
// 实现个鬼
~~P.S. 突然发现，题目提示的最后一行如是写到：“善良的杨村花又提醒你，你的复杂度可能只是「看起来」很对哦”~
好的我终于还是受不了了开始百度 orz
跟我上面提到的一样，确实是将 dfs 改成 bfs，以避免重复的扩展。
在寻找增广路过程中，将发现的新加入相等子图的边但本身已经在相等子图中的点加入队列。
同样的，每次只要找到一条增广路就可以了。对每个点 bfs 找增广路。
在对每个点的 bfs 中，有两种可能的结果：
- 找到增广路
- 没找到，扩展出一棵匈牙利树
若为后者，则接下来需进行顶标修改。
与 dfs 不同，dfs 在修改后需要重新从头开始走一遍以找到新加入的边，而 bfs 则将顶标修改后满足要求的合法边能到达的点直接加入队列，然后直接从这些点开始下一轮的增广路寻找。
这样做避免了重复的跑图过程。这样的话实际上需要把所有边都遍历一遍的点就少了很多，bfs 的平均复杂度与 \(O(n^2)\) 或 \(O(m)\) 的距离更远了。
（因为我无法证明它的复杂度确实达到了 \(O(n)\) ，所以我这里用了“远不及 \(O(n^2)\)” 这样模糊的表述）
顺带一提，看完题解我才发现，自己之前写的代码有很多不必要的部分。
比如初始化 \(l(x)\) 时，其实只需要对其中一个点集初始化，另一个设为 0 就行了。
还有就是每次 dfs 或 bfs 只是找出一条增广路，所以每个点无论如何都是需要 dfs 一次的，并不需要专门用 link 这个数组来判断。
还有就是 pre 数组的设置也有很多细节需要注意，具体可以看代码实现。
在重构和理解后，为 KM 算法的学习划上久违了一个学期的句号吧。
#include <cstdio>
#include <algorithm>
#include <queue>
using namespace std;
typedef long long ll;
const int N=503,M=250005;
const ll INF=1e11;
int n,m,head[N],to[M],nex[M]; ll wei[M];
ll lA[N],lB[N],slack[N];
inline void addEdge(int u,int v,ll w){
static int cc=0;
to[++cc]=v; nex[cc]=head[u]; head[u]=cc; wei[cc]=w;
}
bool vis1[N],vis2[N]; // A,B 中点是否已经被访问过
int Left[N],Right[N]; // 表示 B 中点 i 对应的匹配点是 Left[i], A 中点 i 则对应 Right[i]
int pre[N]; // 存遍历到这个点的上一个点 我们是根据这个点出去的非匹配边将 B 中的这个点加入队列中的
inline void reset(int u){
// 找到增广路后 更新匹配点
// 因为没有了 dfs 的自然回溯，所以只能通过多存下相关信息来实现这一过程
while(u){
Left[u]=pre[u];
swap(u,Right[Left[u]]);
}
}
inline void update_lx(){
ll d=INF;
for(int i=1;i<=n;i++) if(!vis2[i]) d=min(d,slack[i]);
for(int i=1;i<=n;i++){
if(vis1[i]) lA[i]-=d;
if(vis2[i]) lB[i]+=d;
else if(slack[i]!=INF) slack[i]-=d;
}
}
inline void KM(int u){
queue<int> Q; Q.push(u);
while(1){
while(!Q.empty()){
u=Q.front(); vis1[u]=1; Q.pop();
for(int i=head[u];i;i=nex[i]){
if(vis2[to[i]]) continue; // 不要污染了 pre 的设置
if(lA[u]+lB[to[i]]-wei[i]<slack[to[i]]){
slack[to[i]]=lA[u]+lB[to[i]]-wei[i];
// 如果这条边不符合要求 因为它可能是 slack 最小的对应边
// 可能会在后面的循环中直接作为增广路返回 所以同样需要设置 pre
pre[to[i]]=u;
}
if(lA[u]+lB[to[i]]!=wei[i]) continue;
// 如果这条边符合要求 则 slack 一定为 0 一定不会在上面设置 pre
// 所以我们在这里再设置一次以保证正确性
vis2[to[i]]=1; pre[to[i]]=u;
// 这个点已经有匹配点 则继续走下去
// 这个点还没匹配点 则已经找到了增广路
if(!Left[to[i]]) return reset(to[i]);
else Q.push(Left[to[i]]);
}
}
update_lx();
// bfs 的关键部分, 把这次更新后可以新扩展出的点直接加入队列中
for(int i=1;i<=n;i++){
if(slack[i]==0&&!vis2[i]){
// 讨论新出现的相等子图中的点
if(!Left[i]) return reset(i); // 完成增广
else Q.push(Left[i]),vis2[i]=1; // 这个 vis2 让我调了两个小时 吐了
}
}
}
}
int main(){
scanf("%d%d",&n,&m);
while(m--){
int u,v; ll w; scanf("%d%d%lld",&u,&v,&w);
addEdge(u,v,w); lA[u]=max(lA[u],w);
}
for(int i=1;i<=n;i++){
// pre[j] 不用初始化 反正都会被覆盖掉
for(int j=1;j<=n;j++) slack[j]=INF,vis1[j]=vis2[j]=0;
KM(i);
}
ll ans=0;
for(int i=1;i<=n;i++) ans+=lA[i]+lB[i];
printf("%lld\n",ans);
for(int i=1;i<=n;i++) printf("%d ",Left[i]);
return 0;
}
