①什么是斯坦納點？
假設原來已經給定了個點，庫朗等指出需要引進的點數至多為，此種點稱為斯坦納點。過每一斯坦納點，至多有三條邊通過。若為三條邊，則它們兩兩交成120°角；若為兩條邊，則此斯坦納點必為某一已給定的點，且此兩條邊交成的角必大於或等於120°。其中最小的網絡稱為已給定點的集合的最小斯坦納樹，記作SMT。若此SMT的斯坦納點中有等於給定點的點，則稱此SMT為退化的，此給定點稱為退化點。
構造方法：
已知B，C，D，E可知B點為轉軸線段BC繞B順時針旋轉60度得到正三角形，再以頂點F為轉軸，FD構成的線段逆時針旋轉得到新的正三角形頂點G，劣弧DF上任意一點都能和D，F構成三個，相同的，劣弧CB上的點也是。
故將第四點E與G相連接在劣弧上得到一個交點，再由交點與F連接交劣弧CB於一點，即構成了非退化情況下的兩斯坦納點，枚舉得到斯坦納最小生成樹，當與頂點連線不與劣弧有交點時則為該種結構的退化點情況.
②什么是斯坦納樹？
斯坦納樹問題是組合優化學科中的一個問題。將指定點集合中的所有點連通，且邊權總和最小的生成樹稱為最小斯坦納樹（Minimal Steiner Tree），其實最小生成樹是最小斯坦納樹的一種特殊情況。而斯坦納樹可以理解為使得指定集合中的點連通的樹，但不一定最小。
③如何求解最小斯坦納樹？
可以用DP求解，dp[i][state]表示以i為根，指定集合中的點的連通狀態為state的生成樹的最小總權值。
轉移方程有兩重：
第一重，先通過連通狀態的子集進行轉移。
dp[i][state]=min{ dp[i][subset1]+dp[i][subset2] }
枚舉子集的技巧可以用 for(sub=(state-1)&state;sub;sub=(sub-1)&state)。
第二重，在當前枚舉的連通狀態下，對該連通狀態進行松弛操作。
dp[i][state]=min{ dp[i][state], dp[j][state]+e[i][j] }
為什么只需對該連通狀態進行松弛？因為更后面的連通狀態會由先前的連通狀態通過第一重轉移得到，所以無需對別的連通狀態松弛。松弛操作用SPFA即可。
復雜度 O(n*3^k+cE*2^k)
c為SPFA復雜度中的常數，E為邊的數量，但幾乎達不到全部邊的數量，甚至非常小。3^k來自於子集的轉移sum{C(i,n)*2^i} (1<=i<=n)，用二項式展開求一下和。
模版如下：
1 /* 2 * Steiner Tree：求，使得指定K個點連通的生成樹的最小總權值 3 * st[i] 表示頂點i的標記值，如果i是指定集合內第m(0<=m<K)個點，則st[i]=1<<m 4 * endSt=1<<K 5 * dptree[i][state] 表示以i為根，連通狀態為state的生成樹值 6 */ 7 #define CLR(x,a) memset(x,a,sizeof(x)) 8 9 int dptree[N][1<<K],st[N],endSt; 10 bool vis[N][1<<K]; 11 queue<int> que; 12 13 int input() 14 { 15 /* 16 * 輸入，並且返回指定集合元素個數K 17 * 因為有時候元素個數需要通過輸入數據處理出來，所以單獨開個輸入函數。 18 */ 19 } 20 21 void initSteinerTree() 22 { 23 CLR(dptree,-1); 24 CLR(st,0); 25 for(int i=1;i<=n;i++) CLR(vis[i],0); 26 endSt=1<<input(); 27 for(int i=1;i<=n;i++) 28 dptree[i][st[i]]=0; 29 } 30 31 void update(int &a,int x) 32 { 33 a=(a>x || a==-1)? x : a; 34 } 35 36 void SPFA(int state) 37 { 38 while(!que.empty()){ 39 int u=que.front(); 40 que.pop(); 41 vis[u][state]=false; 42 for(int i=p[u];i!=-1;i=e[i].next){ 43 int v=e[i].vid; 44 if(dptree[v][st[v]|state]==-1 || 45 dptree[v][st[v]|state]>dptree[u][state]+e[i].w){ 46 47 dptree[v][st[v]|state]=dptree[u][state]+e[i].w; 48 if(st[v]|state!=state || vis[v][state]) 49 continue; //只更新當前連通狀態 50 vis[v][state]=true; 51 que.push(v); 52 } 53 } 54 } 55 } 56 57 void steinerTree() 58 { 59 for(int j=1;j<endSt;j++){ 60 for(int i=1;i<=n;i++){ 61 if(st[i] && (st[i]&j)==0) continue; 62 for(int sub=(j-1)&j;sub;sub=(sub-1)&j){ 63 int x=st[i]|sub,y=st[i]|(j-sub); 64 if(dptree[i][x]!=-1 && dptree[i][y]!=-1) 65 update(dptree[i][j],dptree[i][x]+dptree[i][y]); 66 } 67 if(dptree[i][j]!=-1) 68 que.push(i),vis[i][j]=true; 69 } 70 SPFA(j); 71 } 72 }
④學習心得
參考09年姜碧野神牛寫的論文《SPFA的優化與應用》，里面提到了一道題——[WC2008]游覽計划。這題讓我立刻聯想到了去年北京賽區的E題，差不多的模型，大概就是在一個圖中求給定的k個點的斯坦納生成樹，給定點的個數k<=10。
首先我們知道，最優解必然是一棵樹，然后這棵樹又是由若干棵子樹合並成的，於是我們可以狀態壓縮，把k個節點的連通狀態用一個二進制數j表示，dp[i][j]表示以i為根和對應狀態為j的節點連通的子樹的最小權值。有兩種轉移方法：
枚舉子樹的形態：dp[ i ][ j ]=min{ dp[ i ][ j ]，dp[ i ][ k ]+dp[ i ][ l ] }，其中k和l是對j的一個划分。
按照邊進行松弛：dp[ i ][ j ]=
min{ dp[ i ][ j ]，dp[ i' ][ j ]+w[ i ][ i' ] }，其中i和i'之間有邊相連。
對於第一種轉移，我們直接枚舉子集就行了。對於第二種轉移，我們仔細觀察可以發現這個方程和最短路的約束條件是很類似的，於是我們可以用spfa或者dij來進行狀態轉移。枚舉子集的復雜度=n*sum{C(k，i)*2^i，0<i=k}=n*3^k，spfa的復雜度為n*2^k。所以總復雜度為O(n*3^k)。
具體實現的時候我試了好幾種不同的方法，一開始是直接把兩種轉移都看成圖中的邊，一遍spfa得出結果，大概如下所示：
1 void spfa(){ 2 while(!Q.empty()){ 3 int x=Q.front()/10000,y=Q.front()%10000; 4 in[x][y]=0; 5 Q.pop(); 6 for(edge *i=Adj[x];i;i=i->nxt) //對當前節點的每條邊都進行松弛操作 7 update(i->v,s[i->v]|y,d[x][y]+i->w); 8 int t=nn-1-y; 9 for(int i=t;i;i=(i-1)&t) //枚舉補集的所有子集，進行松弛操作 10 update(x,y|i,d[x][y]+d[x][i|s[x]]); 11 } 12 }
這么做的復雜度是沒有變的，但是常數非常大，hdu上跑了2500ms才過，仔細一想，我們發現第二松弛操作其實做了很多無用功，考慮能不能進行優化。
第二種松弛操作非常的耗時間，所以我們就不把它加到spfa里面進行轉移，直接在外面進行枚舉，實現更新，避免大量的重復計算。先枚舉連通性j，對於所有的1<=i<=n，我們先進行第一種轉移，既枚舉子集進行更新。如果dp[i][j]被更新了，我們就把它加到隊列里，最后再進行spfa()，這樣按j分層的進行轉移，大概如下：
1 for(int y=0;y<nn;y++) //枚舉連通性 2 for(int x=1;x<=n;x++){ 3 bool flag=0; 4 for(int i=(y-1)&y;i;i=(i-1)&y) //枚舉所有子集，進行第一種轉移 5 flag|=update(x,y,d[x][i|s[x]]+d[x][(y-i)|s[x]]); 6 if(flag) Q.push(x*10000+y); //如果節點被更新則加入隊列 7 spfa(); //spfa進行第二種轉移 8 }
我本來以為這樣會更快一些，結果跑了4700ms = =！頓時吐槽無力。
為啥這樣會更慢呢？我覺的大概是由於spfa()的次數過多，所以導致很多節點被重復的更新了很多次，又產生了大量了重復計算，所以反而更慢了。那么就沒有什么好辦法嗎？仔細一想，我發現進行spfa的時候只需要對當前層的節點進行spfa就行了，不需要整個圖完全松弛一遍，因為更高的層都可以通過枚舉子集而變成若干個更低的層，這樣一次spfa的復雜度一下就降了下來，變成了O(n)級別，大概如下：
1 for(edge *i=Adj[x];i;i=i->nxt) 2 if(update(i->v,y|s[i->v],d[x][y]+i->w)&&y==(y|s[i->v])&&!in[i->v][y]) //只把處於相同層的節點加到隊列中 3 in[i->v][y]=1,Q.push(i->v*10000+y);
這樣修改以后效果果然非常明顯，1000ms就AC了。但還是不夠快，別人最快的能夠達到500ms。於是我baidu了一下，發現他們沒有用spfa！大概就是把第二種轉移表示成了另外一種形式：
dp[ i ][ j ]=min{ dp[ i ][ j ] , dp[ k ][ j ]+d[ k ][ i ] }，其中d[ k ][ i ]表示k到i的最短路。
很容易就能證明這樣寫方程也是對的，於是我們就可以先用floyed預處理出任意兩點間的最短路，然后直接DP。這樣做的總復雜度為O(n^3+n^2*2^k+n*3^k)，這個復雜度並不比上面的方法低，但由於hdu4085的n比較小，所以這樣寫反而比上一種方法要快上不少。但對於
[WC2008]游覽計划、ZOJ 3613 Wormhole Transport這兩道題就不行了，n都達到了100甚至200的大小，這種方法要比前面一種慢。所以最后得出結論，還是前一種方法最穩定 = ^ =
⑤更多習題分享
HDU 4085 Peach Blossom Spring
11年北京賽區的E題，這題有點不同的地方在於，最后的答案可能是一個森林，所以我們要先求出斯坦納樹后進行DP。轉移的時候要注意一點，只有人的個數和房子的個數相等的時候才算合法狀態，所以我們要加一個check()函數進行檢查。
1 #include<cstdio> 2 #include<cstring> 3 #include<vector> 4 #include<queue> 5 #include<algorithm> 6 #define N 60 7 #define INF 2000000 8 using namespace std; 9 struct edge{ 10 int v,w; 11 edge *nxt; 12 }E[2009],*Adj[N],*cur; 13 int n,m,K,nn; 14 int s[N],in[N][1<<10]; 15 int d[N][1<<10],dp[1<<10]; 16 queue<int> Q; 17 void addedge(int u,int v,int w){cur->v=v,cur->w=w,cur->nxt=Adj[u],Adj[u]=cur++;} 18 bool check(int x){ 19 int r=0; 20 for(int i=0;x;i++,x>>=1) 21 r+=(x&1)*(i<K?1:-1); 22 return r==0; 23 } 24 inline bool update(int x,int y,int w){ 25 if(w<d[x][y]) return d[x][y]=w,true; 26 return false; 27 } 28 void spfa(){ 29 while(!Q.empty()){ 30 int x=Q.front()/10000,y=Q.front()%10000; 31 in[x][y]=0; 32 Q.pop(); 33 for(edge *i=Adj[x];i;i=i->nxt) 34 if(update(i->v,y|s[i->v],d[x][y]+i->w)&&y==(y|s[i->v])&&!in[i->v][y]) 35 in[i->v][y]=1,Q.push(i->v*10000+y); 36 37 } 38 } 39 void init(){ 40 cur=E; 41 memset(Adj,0,sizeof(Adj)); 42 memset(s,0,sizeof(s)); 43 scanf("%d%d%d",&n,&m,&K); 44 nn=1<<(2*K); 45 for(int i=1;i<=n;i++) 46 for(int j=0;j<nn;j++) 47 d[i][j]=INF; 48 while(m--){ 49 int u,v,w; 50 scanf("%d%d%d",&u,&v,&w); 51 addedge(u,v,w); 52 addedge(v,u,w); 53 } 54 for(int i=1;i<=K;i++){ 55 s[i]=1<<(i-1),d[i][s[i]]=0; 56 s[n-i+1]=1<<(K+i-1),d[n-i+1][s[n-i+1]]=0; 57 } 58 } 59 int main(){ 60 int T; 61 scanf("%d",&T); 62 while(T--){ 63 init(); 64 for(int y=0;y<nn;y++){ 65 for(int x=1;x<=n;x++){ 66 for(int i=(y-1)&y;i;i=(i-1)&y) 67 d[x][y]=min(d[x][y],d[x][i|s[x]]+d[x][(y-i)|s[x]]); 68 if(d[x][y]<INF) Q.push(x*10000+y),in[x][y]=1; 69 } 70 spfa(); 71 } 72 for(int j=0;j<nn;j++){ 73 dp[j]=INF; 74 for(int i=1;i<=n;i++) dp[j]=min(dp[j],d[i][j]); 75 } 76 for(int i=1;i<nn;i++) 77 if(check(i)) 78 for(int j=i&(i-1);j;j=(j-1)&i) 79 if(check(j)) 80 dp[i]=min(dp[i],dp[j]+dp[i-j]); 81 if(dp[nn-1]>=INF) puts("No solution"); 82 else printf("%d\n",dp[nn-1]); 83 } 84 }
[WC2008]游覽計划
這題要求一棵滿足要求的斯坦納樹，基本上按照上面的做法寫就行了，不過有一點惡心的就是要輸出一組可行方案，所以DP的時候還要記錄一下路徑。
1 #include<cstdio> 2 #include<cstring> 3 #include<vector> 4 #include<queue> 5 #include<algorithm> 6 #define INF 2000000 7 #define N 10 8 using namespace std; 9 int dx[]={0,1,0,-1}, 10 dy[]={1,0,-1,0}; 11 int max_s,n,m; 12 int mat[N][N],st[N][N],vis[N][N],cnt; 13 int d[N][N][1<<N],pre[N][N][1<<N]; 14 bool in[N][N][1<<N]; 15 queue<int> Q; 16 void spfa(){ 17 int x,y,s,tx,ty,ts; 18 while(!Q.empty()){ 19 x=Q.front()/100000; 20 y=(Q.front()-x*100000)/10000; 21 s=Q.front()-x*100000-y*10000; 22 Q.pop(); 23 in[x][y][s]=0; 24 for(int i=0;i<4;i++){ 25 tx=x+dx[i],ty=y+dy[i]; 26 if(tx>=n||ty>=m||tx<0||ty<0) continue; 27 ts=s|st[tx][ty]; 28 if(d[x][y][s]+mat[tx][ty]<d[tx][ty][ts]){ 29 d[tx][ty][ts]=d[x][y][s]+mat[tx][ty]; 30 pre[tx][ty][ts]=x*100000+y*10000+s; 31 if(!in[tx][ty][ts]&&s==ts) in[tx][ty][ts]=1,Q.push(tx*100000+ty*10000+ts); 32 } 33 } 34 } 35 } 36 void go(int x,int y,int s){ 37 vis[x][y]=1; 38 int t=pre[x][y][s],tx,ty,ts; 39 if(!t) return; 40 tx=t/100000; 41 ty=(t-tx*100000)/10000; 42 ts=t-tx*100000-ty*10000; 43 go(tx,ty,ts); 44 if(x==tx&&y==ty) go(x,y,(s-ts)|st[x][y]); 45 } 46 int main(){ 47 //freopen("in.in","r",stdin); 48 scanf("%d%d",&n,&m); 49 for(int i=0;i<n;i++) 50 for(int j=0;j<m;j++){ 51 scanf("%d",&mat[i][j]); 52 if(!mat[i][j]) st[i][j]=1<<(cnt++); 53 } 54 max_s=1<<cnt; 55 for(int i=0;i<n;i++) 56 for(int j=0;j<m;j++){ 57 for(int k=0;k<max_s;k++) 58 d[i][j][k]=INF; 59 if(st[i][j]) d[i][j][st[i][j]]=0; 60 } 61 for(int k=1;k<max_s;k++){ 62 for(int i=0;i<n;i++) 63 for(int j=0;j<m;j++){ 64 if(st[i][j]&&!(st[i][j]&k)) continue; 65 for(int x=(k-1)&k;x;x=(x-1)&k){ 66 int t=d[i][j][x|st[i][j]]+d[i][j][(k-x)|st[i][j]]-mat[i][j]; 67 if(t<d[i][j][k]) d[i][j][k]=t,pre[i][j][k]=i*100000+j*10000+(x|st[i][j]); 68 } 69 if(d[i][j][k]<INF) Q.push(i*100000+j*10000+k),in[i][j][k]=1; 70 } 71 spfa(); 72 } 73 for(int i=0;i<n;i++) 74 for(int j=0;j<m;j++) 75 if(st[i][j]){ 76 printf("%d\n",d[i][j][max_s-1]); 77 go(i,j,max_s-1); 78 for(int x=0;x<n;x++){ 79 for(int y=0;y<m;y++){ 80 if(st[x][y]) putchar('x'); 81 else if(vis[x][y]) putchar('o'); 82 else putchar('_'); 83 } 84 puts(""); 85 } 86 return 0; 87 } 88 }
ZOJ 3613 Wormhole Transport
ZOJ Monthly, June 2012的C題。和HDU 4085差不多，有一點不同的是一個星球可能有很多個工廠，但是含有資源和含有工廠的星球個數都不超過4。還是先狀態壓縮，然后DP求出斯坦納樹。最優的方案有可能是森林，所以我們還要DP，dp[ i ]表示對應的工廠節點和資源節點組成的斯坦樹森林的最優值。那么：
dp[ i ]=min{ dp[ i ]，dp[ j ]+dp[ k ] }，其中j和k為i的一個划分。
這里要注意一點，所有的狀態i、j、k都要滿足一個條件，就是連通的星球上工廠的個數要大於等於資源的個數，這樣才是一個合法的狀態，所以要加一個check()函數。最后再找到所含資源最多，花費最小的合法方案就是答案。
1 #include<cstdio> 2 #include<cstring> 3 #include<queue> 4 #include<vector> 5 #include<algorithm> 6 #define N 209 7 using namespace std; 8 9 struct edge{int v,w;edge *nxt;}E[10009],*Adj[N],*cur; 10 int n,m,nn; 11 int d[N][1<<8],dp[1<<8]; 12 bool in[N][1<<8]; 13 int S[N],P[N],st[N],fac[4],cf,cs; 14 queue<int> Q; 15 void addedge(int u,int v,int w){cur->v=v,cur->w=w,cur->nxt=Adj[u],Adj[u]=cur++;} 16 void up(int &a,int b){if(a==-1||a>b) a=b;} 17 void spfa(){ 18 while(!Q.empty()){ 19 int x=Q.front()/1000,y=Q.front()%1000; 20 Q.pop(); 21 in[x][y]=0; 22 for(edge *i=Adj[x];i;i=i->nxt) 23 if(d[i->v][y|st[i->v]]==-1||d[x][y]+i->w<d[i->v][y|st[i->v]]){ 24 d[i->v][y|st[i->v]]=d[x][y]+i->w; 25 if(y==(y|st[i->v])&&!in[i->v][y]) in[i->v][y]=1,Q.push(i->v*1000+y); 26 } 27 } 28 } 29 bool check(int x){ 30 int t=0; 31 for(int i=0;x;i++,x>>=1) 32 t+=(x&1)*(i<cf?fac[i]:-1); 33 return t>=0; 34 } 35 int cnt(int x){ 36 int r=0; 37 for(int i=0;x;i++,x>>=1) 38 r+=(x&1)*(i<cf?0:1); 39 return r; 40 } 41 int main(){ 42 while(scanf("%d",&n)+1){ 43 cur=E; 44 cf=cs=0; 45 memset(Adj,0,sizeof(Adj)); 46 memset(st,0,sizeof(st)); 47 memset(d,-1,sizeof(d)); 48 memset(dp,-1,sizeof(dp)); 49 int ans=0; 50 for(int i=1;i<=n;i++){ 51 scanf("%d%d",P+i,S+i); 52 if(S[i]&&P[i]) P[i]--,S[i]=0,ans++; 53 if(P[i]) st[i]=1<<cf,fac[cf++]=P[i],d[i][st[i]]=0; 54 } 55 for(int i=1;i<=n;i++) 56 if(S[i]) 57 st[i]=1<<(cf+cs++),d[i][st[i]]=0; 58 nn=1<<(cf+cs); 59 60 scanf("%d",&m); 61 while(m--){ 62 int u,v,w; 63 scanf("%d%d%d",&u,&v,&w); 64 addedge(u,v,w); 65 addedge(v,u,w); 66 } 67 68 for(int y=1;y<nn;y++){ 69 for(int x=1;x<=n;x++){ 70 if(st[x]&&!(st[x]&y)) continue; 71 for(int i=(y-1)&y;i;i=(i-1)&y) 72 if(d[x][i|st[x]]!=-1&&d[x][(y-i)|st[x]]!=-1) 73 up(d[x][y],d[x][i|st[x]]+d[x][(y-i)|st[x]]); 74 if(d[x][y]!=-1) Q.push(x*1000+y),in[x][y]=1; 75 } 76 spfa(); 77 } 78 for(int i=1;i<=n;i++) 79 for(int j=0;j<nn;j++) 80 if(d[i][j]!=-1) 81 up(dp[j],d[i][j]); 82 int num=0,cost=0; 83 for(int i=1;i<nn;i++) 84 if(check(i)){ 85 for(int j=(i-1)&i;j;j=(j-1)&i) 86 if(check(j)&&check(i-j)&&dp[j]!=-1&&dp[i-j]!=-1) 87 up(dp[i],dp[j]+dp[i-j]); 88 int t=cnt(i); 89 if(dp[i]!=-1&&(t>num||(t==num&&dp[i]<cost))) 90 num=t,cost=dp[i]; 91 } 92 printf("%d %d\n",num+ans,cost); 93 } 94 }