[ Minimum Arborescence / zhu_liu ] 朱劉算法 - 最小樹形圖

在有向圖中，給定一個點$r$作為生成樹的根，找出有向圖最小生成樹。

首先我們要能保證從$r$能夠走到圖上的所有點，這樣生成樹才會存在，這很簡單，一次DFS即可，再來是把圖上的所有自環移除，因為一顆樹裡面很明顯是不會有自環的。

之後就是算法的主要步驟了，可以先想一下，除了$r$以外的每一個點都有一條儘可能小的邊指向自己，最好的情況就是我們枚舉每一個點（除了根節點）並找到最小的一條指向這個點的邊，如果這些邊不構成有向環，就形成了一個所求的最小樹形圖。

但是實際上會出現環啊，但是這些環一定是獨立的，為甚麼呢？因為只有$|V|-1$條邊啊，只有是一棵樹的時候才會是連通的狀態。換句話說，如果圖連通了，就一定是最小樹形圖。

我們嘗試去換一些邊，使圖連通，在換的過程中我們總是選擇較小的邊，那麼得到的就是最小樹形圖。你可能會去枚舉一些邊把有向環拆掉，但是這樣的話可能會產生新的有向環，不是一個好做法。

朱劉算法就不直接去換邊，它也不去拆掉環，而是在不增加邊的情況下讓圖連通，怎麼做呢？就是用一個新的點代替原來圖的一個環（也就是所謂的「縮點」，和強連通分量有點像），並且修改跟這個環裡的點有關的邊的權值。

為何要修改邊的權重呢？當我們每更換一個點的入邊的時候我們就要去掉原來那個入邊，於是我們把這個點所有可能的入邊全部減少原來選取的那個入邊的權值，這樣每增加一條入邊無形中就刪去了原來那條邊。

上圖中紅色部分是要進行縮點的有向環

每個環上的點所有可能的入邊全部減少原來選取的那個入邊的權值

接著把環縮成一個點就可以了

假設我們想要把原來縮環之前3那條邊換成4那條邊，那我們換完的結果如下:

可以發現修改邊權後，不需要把邊刪掉，直接去計算選取邊的權重和就會和換邊的結果一樣

朱劉算法主算法的過程就是：找最小入邊->判斷有沒有環(沒有環就退出，算法成功)->縮點，改權值，如此反覆，一般來說為了方便不會去記錄縮點後虛擬節點裡包含了那些點，如果需要找出最小樹形圖包含的邊，就必須要真的去記錄他。

時間複雜度來說的話，用當時論文提出的實作方式，修改邊權的部分為$\ord{|E|}$，縮點最多執行$|V|-1$次，所以總複雜度是$\ord{|V|*|E|}$。
我自己有想了一個$\ord{|E| \; log \; |E|}$的方法，需要用到一種可以合併和把裡面所有元素加上某個值 的heap，又因為每個點最多只會連出去$|V|-1$條邊，也就是heap裡面只有$|V|$個元素是有用的，所以可以在heap大小為$2|V|$時把後$|V|$個元素刪掉，用斐式堆可以做到$\ord{|E|+|V| \; log|V|}$。

以下為$\ord{|V|*|E|}$的code:

	#define INF 0x3f3f3f3f
	template<typename T>
	struct zhu_liu{
	static const int MAXN=110;
	struct edge{
	int u,v;
	T w;
	edge(int u=0,int v=0,T w=0):u(u),v(v),w(w){}
	};
	vector<edge>E;// 0-base
	int pe[MAXN],id[MAXN],vis[MAXN];
	T in[MAXN];
	void init(){E.clear();}
	void add_edge(int u,int v,T w){
	if(u!=v)E.push_back(edge(u,v,w));
	}
	T build(int root,int n){
	T ans=0;
	for(;;){
	for(int u=0;u<n;++u)in[u]=INF;
	for(size_t i=0;i<E.size();++i)
	if(E[i].u!=E[i].v&&E[i].w<in[E[i].v])
	pe[E[i].v]=i,in[E[i].v]=E[i].w;
	for(int u=0;u<n;++u)//無解
	if(u!=root&&in[u]==INF)return -INF;
	int cntnode=0;
	memset(id,-1,sizeof(int)*n);
	memset(vis,-1,sizeof(int)*n);
	for(int u=0;u<n;++u){
	if(u!=root)ans+=in[u];
	int v=u;
	for(;vis[v]!=u&&id[v]==-1&&v!=root;v=E[pe[v]].u)
	vis[v]=u;
	if(v!=root&&id[v]==-1){
	for(int x=E[pe[v]].u;x!=v;x=E[pe[x]].u)
	id[x]=cntnode;
	id[v]=cntnode++;
	}
	}
	if(!cntnode)break;//無環
	for(int u=0;u<n;++u)if(id[u]==-1)id[u]=cntnode++;
	for(size_t i=0;i<E.size();++i){
	int v=E[i].v;
	E[i].u=id[E[i].u];
	E[i].v=id[E[i].v];
	if(E[i].u!=E[i].v)E[i].w-=in[v];
	}
	n=cntnode;
	root=id[root];
	}
	return ans;
	}
	};

view raw zhu_liu.cpp hosted with ❤ by GitHub

接著是$\ord{|E| \; log^2|E|}$的code，使用啟發式合併(感謝59491、編譯器幫忙debug):

	template<typename T>
	struct zhu_liu{
	static const int MAXN=110;
	struct edge{
	int u,v;
	T w;
	edge(int u=0,int v=0,T w=0):u(u),v(v),w(w){}
	bool operator<(const edge &b)const{
	return w>b.w;
	}
	};
	typedef pair<priority_queue<edge>,T> heap;
	#define x first
	#define y second
	heap *pq[MAXN*2],mem[MAXN];//靜態記憶體
	int st[MAXN*2],id[MAXN*2];
	edge E[MAXN*2];
	void init(int n){
	for(int i=1;i<=n;++i){
	pq[i]=&mem[i],mem[i]=heap();
	E[i]=edge();
	st[i]=id[i]=i;
	}
	}
	void add_edge(int u,int v,T w){
	if(u!=v)pq[v]->x.push(edge(u,v,w));
	}
	heap *merge(heap *a,heap *b){//啟發式合併
	if(a->x.size()<b->x.size())swap(a,b);
	while(b->x.size()){
	edge e=b->x.top();e.w+=b->y-a->y;
	a->x.push(e),b->x.pop();
	}
	return a;
	}
	int find(int x,int *st){
	return st[x]==x?x:st[x]=find(st[x],st);
	}
	T build(int root,int n){
	T ans=0;int N=n,all=n;
	for(int i=1;i<=N;++i){
	if(i==root||pq[i]->x.empty())continue;
	while(pq[i]->x.size()){
	E[i]=pq[i]->x.top();pq[i]->x.pop();
	if(find(E[i].u,id)!=find(i,id))break;
	}
	if(find(E[i].u,id)==find(i,id))continue;
	ans+=(E[i].w+=pq[i]->y);
	if(find(E[i].u,st)==find(i,st)){
	pq[i]->y-=E[i].w;
	pq[++N]=pq[i],id[N]=N;
	for(int u=find(E[i].u,id);u!=i;u=find(E[u].u,id)){
	pq[u]->y-=E[u].w,id[find(u,id)]=N;
	pq[N]=merge(pq[N],pq[u]);
	}
	st[N]=find(i,st);
	id[find(i,id)]=N;
	}else st[find(i,st)]=find(E[i].u,st),--all;
	}
	return all==1?ans:-INT_MAX;//圖不連通就無解
	}
	};

view raw zhu_liu_logMlogM.cpp hosted with ❤ by GitHub

接著是$\ord{|E| \; log|E|}$的code，使用skew heap:

	template<typename T>
	struct zhu_liu{
	static const int MAXN=110,MAXM=10005;
	struct node{
	int u,v;
	T w,tag;
	node *l,*r;
	node(int u=0,int v=0,T w=0):u(u),v(v),w(w),tag(0),l(0),r(0){}
	void down(){
	w+=tag;
	if(l)l->tag+=tag;
	if(r)r->tag+=tag;
	tag=0;
	}
	}mem[MAXM];//靜態記憶體
	node *pq[MAXN*2],*E[MAXN*2];
	int st[MAXN*2],id[MAXN*2],m;
	void init(int n){
	for(int i=1;i<=n;++i){
	pq[i]=E[i]=0;
	st[i]=id[i]=i;
	}m=0;
	}
	node *merge(node *a,node *b){//skew heap
	if(!a||!b)return a?a:b;
	a->down(),b->down();
	if(b->w<a->w)return merge(b,a);
	swap(a->l,a->r);
	a->l=merge(b,a->l);
	return a;
	}
	void add_edge(int u,int v,T w){
	if(u!=v)pq[v]=merge(pq[v],&(mem[m++]=node(u,v,w)));
	}
	int find(int x,int *st){
	return st[x]==x?x:st[x]=find(st[x],st);
	}
	T build(int root,int n){
	T ans=0;int N=n,all=n;
	for(int i=1;i<=N;++i){
	if(i==root||!pq[i])continue;
	while(pq[i]){
	pq[i]->down(),E[i]=pq[i];
	pq[i]=merge(pq[i]->l,pq[i]->r);
	if(find(E[i]->u,id)!=find(i,id))break;
	}
	if(find(E[i]->u,id)==find(i,id))continue;
	ans+=E[i]->w;
	if(find(E[i]->u,st)==find(i,st)){
	if(pq[i])pq[i]->tag-=E[i]->w;
	pq[++N]=pq[i],id[N]=N;
	for(int u=find(E[i]->u,id);u!=i;u=find(E[u]->u,id)){
	if(pq[u])pq[u]->tag-=E[u]->w;
	id[find(u,id)]=N;
	pq[N]=merge(pq[N],pq[u]);
	}
	st[N]=find(i,st);
	id[find(i,id)]=N;
	}else st[find(i,st)]=find(E[i]->u,st),--all;
	}
	return all==1?ans:-INT_MAX;//圖不連通就無解
	}
	};

view raw zhu_liu_logM.cpp hosted with ❤ by GitHub