[ Median-of-Medians Algorithm ] 中位數演算法

這是一個可以在保證線性時間(c++ std::nth_element是隨機演算法)找出一個序列中第k大元素的演算法，網路上已經有不少教學，但是很多人都認為常數太大因此缺乏實作。

教學文在此：http://tmt514-blog.logdown.com/posts/484313-divide-and-conquer-method-iii-find-the-median

其實我高中時就想要試著去時做看看，但是因為那時的程式能力太差的關係，做出來的東西一直有bug，後來去忙其他事情後就被我忘掉了。最近因為學長面試有被問到一樣的問題跑來問我，才慢慢想起來有一份沒寫完的code，於是今天抱著不管怎樣都要寫出來的精神把他寫完了：

	#include<algorithm>
	template<typename IT,typename CMP>
	inline IT find_mid(const IT first,const IT last,CMP cmp){
	std::sort(first,last,cmp);
	return first+(last-first)/2;
	}
	template<typename IT,typename CMP>
	inline IT partition(IT first,IT last,IT pivot,CMP cmp){
	if(last-first<1)return pivot;
	std::iter_swap(--last,pivot);
	pivot=last;
	for(;;){
	while(cmp(*first,*pivot))++first;
	--last;
	while(cmp(*pivot,*last))--last;
	if(!(first<last)){
	std::iter_swap(first,pivot);
	return first;
	}
	std::iter_swap(first,last);
	++first;
	}
	}
	template<typename IT,typename CMP>
	IT kth(const IT first,const IT last,int k,CMP cmp){
	if(k==0||k>last-first) return last;
	IT mid_e=first,i=first;
	for(;i+5<last;i+=5){
	std::iter_swap(mid_e++,find_mid(i,i+5,cmp));
	}
	std::iter_swap(mid_e++,find_mid(i,last,cmp));
	int five_mid=mid_e-first;
	mid_e=(five_mid==1)?first:kth(first,mid_e,five_mid/2,cmp);
	IT pos=partition(first,last,mid_e,cmp);
	if(pos-first==k-1) return pos;
	if(pos-first<k) return kth(pos+1,last,k-(pos-first+1),cmp);
	return kth(first,pos,k,cmp);
	}

view raw kth.cpp hosted with ❤ by GitHub

	template<typename IT,typename CMP>
	void SM_sort(IT first,IT last,CMP cmp){ // 幫他寫了一個sort做測試
	if(last-first<=1)return;
	IT mid=kth(first,last,(last-first)/2,cmp);
	SM_sort(first,mid,cmp);
	SM_sort(mid+1,last,cmp);
	}

view raw SM_sort.cpp hosted with ❤ by GitHub

[ C++ std::sort python implement ] C++std::sort Python實作

最近有一個很靠北的課需要用python寫一些C++很容易做到的東西。
這次我們被要求寫一個quick sort，但是他要求的做法實在太糟糕了，於是我就參考了C++ std::sort來實做了一個quick sort(不包含heap sort的部分)：

	import random
	def quick_sort_implace(s):
	'''C++ std::sort implement without heap_sort'''
	first=0
	last=len(s)
	if(first!=last):
	__introsort_loop(s,first,last)
	__final_insertion_sort(s,first,last)
	def __final_insertion_sort(s,first,last):
	if(last-first>16):
	__insertion_sort(s,first,first+16)
	__unguarded_insertion_sort(s,first+16,last);
	else:
	__insertion_sort(s,first,last)
	def __insertion_sort(s,first,last):
	if(first==last):
	return
	for i in range(first+1,last):
	if(s[i]<s[first]):
	val=s[i]
	copy_backward(s,first,i,i+1)
	s[first]=val
	else:
	__unguarded_linear_insert(s,i)
	def copy_backward(s,first,last,result):
	while(first!=last):
	result-=1
	last-=1
	s[result]=s[last]
	return result
	def __unguarded_insertion_sort(s,first,last):
	for i in range(first,last):
	__unguarded_linear_insert(s,i)
	def __unguarded_linear_insert(s,last):
	val=s[last]
	__next=last
	__next-=1
	while(val<s[__next]):
	s[last]=s[__next]
	last=__next
	__next-=1
	s[last]=val
	def __introsort_loop(s,first,last):
	while(last-first>16):
	cut=__unguarded_partition_pivot(s,first,last)
	__introsort_loop(s,cut,last)
	last=cut
	def __unguarded_partition_pivot(s,first,last):
	mid=random.randint(first,last-1)
	__move_median_first(s,first,mid,last-1)
	return __unguarded_partition(s,first+1,last,s[first])
	def __move_median_first(s,a,b,c):
	if(s[a]<s[b]):
	if(s[b]<s[c]):
	s[a],s[b]=s[b],s[a]
	elif(s[a]<s[c]):
	s[a],s[c]=s[c],s[a]
	elif(s[a]<s[c]):
	return
	elif(s[b]<s[c]):
	s[a],s[c]=s[c],s[a]
	else:
	s[a],s[b]=s[b],s[a]
	def __unguarded_partition(s,first,last,pivot):
	while(True):
	while(s[first]<pivot):
	first+=1
	last-=1
	while(pivot<s[last]):
	last-=1
	if(not(first<last)):
	return first
	s[first],s[last]=s[last],s[first]
	first+=1
	# 以下為測試程式的部分
	def is_sorted(s):
	'''檢測序列s是否已排序'''
	n=len(s)
	i=1
	while(i<n):
	if(s[i]<s[i-1]):
	return False
	i+=1
	return True
	if __name__ == '__main__':
	import time
	s=[3 for n in range(1000000)]
	start=time.time()
	quick_sort_implace(s)
	end=time.time()
	print(is_sorted(s),end-start)
	s=[random.randint(0,1000000) for n in range(1000000)]
	start=time.time()
	quick_sort_implace(s)
	end=time.time()
	print(is_sorted(s),end-start)
	s=[random.randint(0,100) for n in range(1000000)]
	start=time.time()
	quick_sort_implace(s)
	end=time.time()
	print(is_sorted(s),end-start)
	s=[n for n in range(1000000)]
	start=time.time()
	quick_sort_implace(s)
	end=time.time()
	print(is_sorted(s),end-start)

view raw quick_sort.py hosted with ❤ by GitHub

[ inorder postorder construct tree ] 用中序、後序建樹保證$\ord{N}$的方法

昨天我學長因為要面試，所以努力地刷leetcode的題目，寫到了一題:

• 106. Construct Binary Tree from Inorder and Postorder Traversal  

他雖然AC了但是用的方法不太好，因此跑來問我有沒有看起來很帥速度又快的方法。

因為網路上的code都用一些很慢的方法來建樹，很多都$\ord{N^2}$，雖然也有看似$\ord{N}$的方法，但是用到像是unordered_map之類常數很大也不保證複雜度是$\ord{1}$(codeforces會卡內建hash)。

因為我查到的code都太糟糕了，因此我決定自己寫一個複雜度保證為$\ord{N}$又好寫的方法。

首先是找$root$的部分，因為有給後序的關係很容易就是到是誰了，但是找左右子樹就會出問題。經過觀察，只需要在dfs的過程用stack紀錄一些特別點就可以好好地維護左右子樹。

假設現在dfs在點$u$，stack紀錄的點就是從$root$到$u$的路徑所有點中，除了那些左小孩也在該路徑中的點之外的所有點，有點複雜看個圖會比較明白，紅色的點就是記錄在stack中的點。

至於為什麼記錄這些點就可以在dfs時判斷現在是不是NULL的節點，以及如果給的是preorder的情況就交給讀者思考了
以下附上程式碼:

	class Solution{
	int i,j;
	stack<int,vector<int>> st;
	const vector<int> *in,*po;
	TreeNode *dfs(){
	if(j==-1||(st.size()&&st.top()==in->at(i)))return 0;
	st.emplace(po->at(j));
	TreeNode *o=new TreeNode(po->at(j--));
	o->right=dfs();
	st.pop(),--i;
	o->left=dfs();
	return o;
	}
	public:
	TreeNode* buildTree(vector<int>& inorder, vector<int>& postorder) {
	i=j=int(inorder.size())-1;
	stack<int,vector<int>>().swap(st);
	in=&inorder,po=&postorder;
	return dfs();
	}
	};

view raw construct_tree.cpp hosted with ❤ by GitHub

[ Source code beautifier / syntax highlighter ] 在網頁/blog中插入彩色程式碼

先看看結果吧:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

#define x first #define y second #include<bits/stdc++.h> using namespace std; #define X(){\ sdfgsdfg;\ sdfgsdfg;\ } int main(){ //asdfasdfdfsdfghd\ asdfasdfasdf\ asdfasdfsdf wchar_t wc; cout<<"Jinkela"<<'\n'; cout<<R"jinkela( 7122)jinkela"; cout<<L"adsfasdf"<<endl; return 0; /* asdf */ }

這是使用http://hilite.me/ Style=monokai的結果，個人覺得效果不錯，只是raw string和一些比較難實作的東西沒有支援而已，其他的都還算可以

我已經把這個網址加到我的學習連結裡面了，C/C++ syntax highlighter (Style選monokai)那個。用法就是貼上程式碼，設定好語言和style，把產生的html貼在你想貼的位置，蠻簡單的

我在blog中如果程式碼量比較少我覺得沒必要加入模板也會用這個方法來貼code

[ Steiner tree problem in graphs ] 斯坦納樹

斯坦納樹問題是一個世界知名的NP-hard問題。在圖論上的斯坦納樹是對於一張無向圖$G=(V,E)$以及一個點集合$P \subseteq V$，通常會稱$P$集合為$terminal \; set$，對於每條邊$e=(u,v) \in E$，令$w(e)$表示它的權重。我們的目標是要從$G$的所有子圖中找出一棵生成樹$T=(V',E')$，使得$P \subseteq V'$且$\sum_{e \in E'} w(e)$最小。

簡單來說就是在圖$G$上找一棵子樹，可以把$P$中的點連通起來，且邊權總和最小

如果我們枚舉所有子圖，對每個子圖做最小生成樹演算法，就一定可以找到斯坦納樹，但是複雜度是$\ord{(\abs E + \abs V log \abs V ) \times 2^{\abs V}}$，非常糟糕。

如果$w(e)>0 ,e \in E$，且$\abs P \ll \abs V$，我們可以找到一個動態規劃的方法：
令$dp[S][i]$表示以點$i$為根，以$S \subseteq P$為$terminal \; set$構造出來的斯坦納樹，這樣我們最後的答案就會是$dp[P][u \in P]$

狀態轉移式可以寫成這樣

• $dp[S][i]=min(dp[T][j]+dp[S-T][j]+dis(i,j)\; : \; j \in V,T \subset S)$
  $dis(i,j)$表示$i \sim j$的最短路徑

 

任兩點間的最短路徑可以用floyd在$\ord {\abs V^3}$預先算出來，狀態有$2^{\abs P}\times \abs V$個，狀態轉移為$\ord{\abs V \times 枚舉子集合的時間}$，因此總複雜度為$\ord{\abs V^3+2^{\abs P} \times \abs V^2 \times 枚舉子集合的時間 }$

其中 $2^{\abs P} \times 枚舉子集合的時間$ 只是粗略的計算，實際上它是
$$\sum_{i=1}^{\abs P} \binom{\abs P}{i} \times (2^i -1) \simeq \sum_{i=0}^{\abs P} \binom{\abs P}{i} \times 2^i = (1+2)^{\abs P} = 3^{\abs P}$$ 因此總複雜度可以表示為$\ord{V^3+3^{\abs P} \times \abs V^2}$，但是這其實還可以優化，令$H[j] = min(dp[T][j]+dp[S-T][j] \; : \;T \subset S)$
則$dp[S][i]=min(H[j]+dis(i,j)\; : \;j \in \abs V)$
$H$是可以被預先算出來的，因此總複雜度就降為$\ord{\abs V^3 + \abs V 3^{\abs P}+\abs V^2 2^{\abs P}}$
以下附上程式碼:

	//n個點，其中r個要構成斯坦納樹
	//p表示要構成斯坦納樹的點集
	#define REP(i,n) for(int i=0;i<(int)n;++i)
	const int MAXN=30,MAXM=8;// 0-base
	const int INF=0x3f3f3f3f;
	int dp[1<<MAXM][MAXN];
	int g[MAXN][MAXN];//圖
	void init(){memset(g,0x3f,sizeof(g));}
	void add_edge(int u,int v,int w){
	g[u][v]=g[v][u]=min(g[v][u],w);
	}
	void steiner(int n,int r,int *p){
	REP(k,n)REP(i,n)REP(j,n)
	g[i][j]=min(g[i][j],g[i][k]+g[k][j]);
	REP(i,n)g[i][i]=0;
	REP(i,r)REP(j,n)dp[1<<i][j]=g[p[i]][j];
	for(int i=1;i<(1<<r);++i){
	if(!(i&(i-1)))continue;
	REP(j,n)dp[i][j]=INF;
	REP(j,n){
	int tmp=INF;
	for(int s=i&(i-1);s;s=i&(s-1))
	tmp=min(tmp,dp[s][j]+dp[i^s][j]);
	REP(k,n)dp[i][k]=min(dp[i][k],g[j][k]+tmp);
	}
	}
	}

view raw Steiner_tree_floyd.cpp hosted with ❤ by GitHub

有的時候圖是稀疏圖，也就是$\ord V=\ord E$，這種時候用這種算法效率其實不是很好，我們可以在dp的過程才用一些單源最短路徑算法算出最短路徑，這樣複雜度可以變成 $$\ord{\abs V 3^{\abs P}+ShortestPath(G) 2^{\abs P}}$$ 其中$ShortestPath(G)$是在圖$G$中計算最短路徑的時間，用dijkstra的話是$\ord{\abs E+\abs V log \abs V}$，這裡我用SPFA實作:

	//n個點，其中r個要構成斯坦納樹
	//p表示要構成斯坦納樹的點集
	#define x first
	#define y second
	#define REP(i,n) for(int i=0;i<n;++i)
	const int MAXN=30,MAXM=8;// 0-base
	const int INF=0x3f3f3f3f;
	int dp[1<<MAXM][MAXN];
	vector<pair<int,int> > g[MAXN];//圖
	void init(int n){REP(i,n)g[i].clear();}
	void add_edge(int u,int v,int w){
	g[u].push_back({v,w});
	g[v].push_back({u,w});
	}
	queue<int> Q;
	bool inq[MAXN];
	void SPFA(int *d){
	while(Q.size()){
	int u=Q.front();Q.pop(),inq[u]=0;
	for(auto &i:g[u]){
	if(d[i.x]>d[u]+i.y){
	d[i.x]=d[u]+i.y;
	if(inq[i.x])continue;
	Q.push(i.x),inq[i.x]=1;
	}
	}
	}
	}
	void steiner(int n,int r,int *p){
	memset(inq,0,sizeof(inq));
	memset(dp,0x3f,sizeof(dp));
	REP(i,r)dp[1<<i][p[i]]=0;
	for(int i=1;i<(1<<r);++i){
	REP(j,n){
	for(int s=i&(i-1);s;s=i&(s-1))
	dp[i][j]=min(dp[i][j],dp[s][j]+dp[i^s][j]);
	if(dp[i][j]<INF)Q.push(j),inq[j]=1;
	}
	SPFA(dp[i]);
	}
	}

view raw Steiner_tree_SPFA.cpp hosted with ❤ by GitHub

[ gcc Built-in Functions for binary ] gcc內建處理二進位函數

以下介紹的函數都是非標準的函數，他只能在GCC底下使用，不過一般的比賽環境都是支援的，所以熟記起來可以增加寫位元運算的效率

1. int __builtin_ffs (unsigned int x)
   int __builtin_ffsl (unsigned long)
   int __builtin_ffsll (unsigned long long)
   
   • 返回右起第一個1的位置
   • Returns one plus the index of the least significant 1-bit of x, or if x is zero, returns zero.
2. int __builtin_clz (unsigned int x)
   int __builtin_clzl (unsigned long)
   int __builtin_clzll (unsigned long long)
   
   • 返回左起第一個1之前0的個數
   • Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
    
3. int __builtin_ctz (unsigned int x)
   int __builtin_ctzl (unsigned long)
   int __builtin_ctzll (unsigned long long)
   
   • 返回右起第一個1之後的0的個數
   • Returns the number of trailing 0-bits in x, starting at the least significant bit position. If x is 0, the result is undefined.
    
4. int __builtin_popcount (unsigned int x)
   int __builtin_popcountl (unsigned long)
   int __builtin_popcountll (unsigned long long)
   
   • 返回1的個數
   • Returns the number of 1-bits in x.
    
5. int __builtin_parity (unsigned int x)
   int __builtin_parityl (unsigned long)
   int __builtin_parityll (unsigned long long)
   
   • 返回1的個數的奇偶性(1的個數 mod 2的值)
   •  Returns the parity of x, i.e. the number of 1-bits in x modulo 2. 

這種內建函數其實非常多，這邊有附上一個GCC內建函數的列表，有興趣的朋友可以參考
https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html

當你在GCC環境下，想直接用01寫二進位的東西，其實有簡單的方法:

• cout<<0b1010101;
• cout<<0b1010101LL;

這樣你的編譯器應該會警告你說這是GCC內建的東西(C++14以後的版本會支援它)，但是還是會正確執行，都是85，個人覺得蠻方便的

如果你是用C++11，可以用stoi,stol,stoll,stoul,stoull等函數來把二進位字串轉成int,long,long long,unsigned long,unsigned long long等，可以設定他用二進位轉換，像是這樣:

• cout<<stoi("1010101",NULL,2);
• cout<<stoll("1010101",NULL,2);

 另外有些編譯器會在<algorithm>實作std::__lg(n)這個函數，他會回傳$\floor{log_2n}$的值，可以參考這個
http://stackoverflow.com/questions/40434664/what-is-std-lg

[ Amortized analysis - Potential method ] 平攤分析 - 勢能方法

對於一個stack操作，pop和push都是$\ord{1}$的，這很好理解，現在定義了一個新的操作，pop_all，表示pop stack內的所有元素，顯然這是一個$\ord{N}$的操作。那我們進行一連串的push、pop、pop_all的複雜度上界是多少呢？

根據big-O的特性，因為pop_all是$\ord{N}$的，我們把每個操作都當作$\ord{N}$來看，執行$N$次操作的總複雜度就會是$\ord{N^2}$。這沒有錯，但怎麼想都覺得怪怪的，怎麼可能一直執行pop_all呢？執行一次pop_all之後stack就沒有元素了耶！

這三種操作不是平行的，而是互相影響的。換言之，你每次的push創造“機會”給pop和pop_all。pop和pop_all才能“消費”著這些機會，不存在無限制的消費。

因此這個複雜度是高估的，那要究竟怎麼真正去計算$N$次操作的總複雜度呢？

勢能方法

對於一個資料結構$DS$，我們定義$\Phi(DS)$表示$DS$的“勢能”，而且要保證在任何情況下$\Phi(DS) \geq 0$。

通常$\Phi(DS)$我們會定義他是$DS$的某個性質，像是元素個數啦，如果是二元樹的話可能是所有節點左子樹大小減右子樹大小的絕對值的總和啊、或是葉節點的個數啊...各種東西都可以當作是$DS$的勢能。

 那$\Phi(DS)$能用來幹甚麼？我們定義$\Phi_0$表示$DS$在還沒有進行任何操作時的勢能，假設有$N$個操作，第$i$個操作的時間花費為$t_i$，這個操作結束之後的勢能為$\Phi_i$，$i>0$，我們定義第$i$個操作的均攤時間$a_i=t_i+C(\Phi_i-\Phi_{i-1})$，這裡$C>0$是一個常數

可以發現總均攤花費時間：
$$\sum^{N}_{i=1}a_i=\sum^{N}_{i=1}(t_i+C(\Phi_i-\Phi_{i-1})) \\ =(t_1+t_2+...+t_n)+C(-\Phi_0+(\Phi_1-\Phi_1)+(\Phi_2-\Phi_2)+...+(\Phi_{n-1}-\Phi_{n-1})+\Phi_n) \\ =\sum^{N}_{i=1} t_i +C(\Phi_N-\Phi_0)$$
最後得到：
$$\sum^{N}_{i=1}t_i=\sum^{N}_{i=1}a_i-C(\Phi_N-\Phi_0)$$
這個特性告訴我們，只要好好選擇$\Phi(DS)$函數，就可以在$\ord{\sum^{N}_{i=1}t_i}$很難直接求的情況下透過$\ord{\sum^{N}_{i=1}a_i-C(\Phi_N-\Phi_0)}$求出$\ord{\sum^{N}_{i=1}t_i}$。

證明範例

有了勢能方法，可以很輕鬆的用它來計算一些資料結構操作的均攤複雜度。

以剛剛stack的例子來說，我們設定stack $S$它的勢能函數$\Phi(S)$為$S$的元素個數，可以確定$\Phi_0=0$且$\Phi(S) \geq 0$。

我們設一次的push、一次的pop花費一單位的時間，並設常數$C=1$，在第$i$次操作：

• 如果是push操作的話$t_i=1$，stack的元素個數增加$1$
  因此$\Phi_i-\Phi_{i-1}=1$
  $a_i=t_i+\Phi_i-\Phi_{i-1}=1+1=2$
• 如果是pop操作的話$t_i=1$，stack的元素個數減少$1$
  因此$\Phi_i-\Phi_{i-1}=-1$
  $a_i=t_i+\Phi_i-\Phi_{i-1}=1-1=0$
•  如果是pop_all操作的話$t_i=n$，stack的元素個數減少$n$
  因此$\Phi_i-\Phi_{i-1}=-n$
  $a_i=t_i+\Phi_i-\Phi_{i-1}=n-n=0$

$a_i$的最大值是$2$，經過$N$次操作之後$\Phi_N-\Phi_0$的最小值為$0$
可以知道：
$$\ord{\sum^{N}_{i=0}t_i}=\ord{\sum^{N}_{i=1}a_i-(\Phi_N-\Phi_0)}=\ord{2N+0}=\ord{N}$$
因此經過$N$次stack的任何有效操作之後，總花費的時間為$\ord{N}$，這才是我們滿意的結果。

對了，通常來說$\Phi_0$都會是0，因此在大部分的情況下：
$$\ord{\sum^{N}_{i=0}t_i}=\ord{\sum^{N}_{i=1}a_i}$$
所以大部分的證明都會忽略掉$\Phi_N-\Phi_0$的部分

[ 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