[ 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