标签：
0. 引子
基础的算法和数据结构已经学习的差不多了，上学期期末就打算重点研究研究STOC和FOCS上面的论文。
做这件事情的初衷是了解别人是如何改进原有算法的，搞清楚目前比较热的算法问题有哪些，更重要的是acm的很多算法或者
书里的算法都是别人整理的，很多年以前的了，学习新东西总会有很多收获的。
关于算法，很多人认为不需要了解太多。大二以前吧，我也是这么认为的，大二以后我就不这么想了。
真的，算法是一件很神奇的事情。不了解的人永远不懂，你写的代码没用到你学习的算法只能说明一个问题——你做的东西太太太简单了。
1. 简介
SSSP: Single Source Shortest Paths.
APSP: All Pairs Shortest Paths.
这篇论文质量挺高的，基本没有错误，同时算法也确实可以提高效率。两位作者也一直致力于这方面的研究。
forward-backward(以下简称FB)算法主要用来解单源最短路径问题。
对于SSSP问题，熟知的算法主要有Dijkstrea(我认为文中所说的Dijkstra并不是O(n^2)的，而是类似于优先级队列实现的SPFA的)，使用Fibo堆（见算法导论）时间复杂度O(m+nlgn)。
而文中提到的Spira则是我们不熟知的算法，效率也很不错。这个算法最早在1950-1960年提出。Spira算法成立的前提是，对于每个结点u的邻接链表，按照边的权重单调非递减顺序排序。
新的SSSP算法基于这样一个前提：对于K_n的图，有很高的概率仅需要检测Omega(nlgn)条边。
简单说，由于排序后的有序关系，我们可以建立某些限制条件，使得有些边不需要检测。这边论文的精华就是限制条件建立的巧妙，而且恰恰可以保证最短路径的性质。
2. Spira算法
Spira算法每次while循环中从队列P取出的边(u->v)，都是当前s起始的最短路径。如果v不在S内，则说明该路径就是s到v的最短路径。
每当从队列释放一条路径时，就对u进行一次forward；每当把v加入S中，就对v进行一次forward。
forward操作可以理解成对SPT（Shortest Paths Tree）的拓展操作，这棵树的根节点为s。
理解清楚forward操作，基本上也就理解了Spira算法了。算法源代码如下
1 /* Spira */ 2 #include <iostream> 3 #include <string> 4 #include <map> 5 #include <queue> 6 #include <set> 7 #include <stack> 8 #include <vector> 9 #include <deque> 10 #include <algorithm> 11 #include <cstdio> 12 #include <cmath> 13 #include <ctime> 14 #include <cstring> 15 #include <climits> 16 #include <cctype> 17 #include <cassert> 18 #include <functional> 19 #include <iterator> 20 #include <iomanip> 21 using namespace std; 22 //#pragma comment(linker,"/STACK:102400000,1024000") 23 24 #define sti set<int> 25 #define stpii set<pair<int, int> > 26 #define mpii map<int,int> 27 #define vi vector<int> 28 #define pii pair<int,int> 29 #define vpii vector<pair<int,int> > 30 #define rep(i, a, n) for (int i=a;i<n;++i) 31 #define per(i, a, n) for (int i=n-1;i>=a;--i) 32 #define clr clear 33 #define pb push_back 34 #define mp make_pair 35 #define fir first 36 #define sec second 37 #define all(x) (x).begin(),(x).end() 38 #define SZ(x) ((int)(x).size()) 39 #define lson l, mid, rt<<1 40 #define rson mid+1, r, rt<<1|1 41 42 typedef struct node_t { 43 int u, v, w; 44 node_t() {} 45 node_t(int u_, int v_, int w_): 46 u(u_), v(v_), w(w_) {} 47 friend bool operator< (const node_t& a, const node_t& b) { 48 return a.w > b.w; 49 } 50 } node_t; 51 52 const int INF = 0x1f1f1f1f; 53 const int maxv = 205; 54 const int maxe = maxv * maxv + 5; 55 node_t ND[maxe]; 56 int dis[maxv]; 57 bool inS[maxv]; 58 int nv, ne; 59 int V[maxe], W[maxe], nxt[maxe]; 60 int head[maxv], head_[maxv]; 61 62 void addEdge(int u, int v, int w) { 63 V[ne] = v; 64 W[ne] = w; 65 nxt[ne] = head_[u]; 66 head_[u] = ne++; 67 } 68 69 void forward(priority_queue<node_t>& Q, int u) { 70 int& k = head[u]; 71 72 if (k != -1) { 73 Q.push(node_t(u, V[k], dis[u]+W[k])); 74 k = nxt[k]; 75 } 76 } 77 78 void spira(int s = 1) { 79 priority_queue<node_t> Q; 80 node_t nd; 81 82 memset(inS, false, sizeof(inS)); 83 memset(dis, INF, sizeof(dis)); 84 memcpy(head, head_, sizeof(head)); 85 inS[s] = true; 86 dis[s] = 0; 87 forward(Q, s); 88 89 while (!Q.empty()) { 90 nd = Q.top(); 91 Q.pop(); 92 forward(Q, nd.u); 93 if (!inS[nd.v]) { 94 inS[nd.v] = true; 95 dis[nd.v] = nd.w; 96 forward(Q, nd.v); 97 } 98 } 99 } 100 101 void input() { 102 int m; 103 104 scanf("%d %d", &nv, &m); 105 rep(i, 0, m) 106 scanf("%d %d %d", &ND[i].u, &ND[i].v, &ND[i].w); 107 108 ne = 0; 109 memset(head_, -1, sizeof(head_)); 110 111 // pre sort 112 sort(ND, ND+m); 113 rep(i, 0, m) 114 addEdge(ND[i].u, ND[i].v, ND[i].w); 115 } 116 117 void solve() { 118 spira(); 119 120 rep(i, 1, nv+1) 121 printf("%d: %d\n", i, dis[i]); 122 } 123 124 int main() { 125 ios::sync_with_stdio(false); 126 #ifndef ONLINE_JUDGE 127 freopen("data.in", "r", stdin); 128 freopen("data.out", "w", stdout); 129 #endif 130 131 input(); 132 solve(); 133 134 #ifndef ONLINE_JUDGE 135 printf("time = %d.\n", (int)clock()); 136 #endif 137 138 return 0; 139 }
3. 验证SPT
验证SPT其实就是验证算法正确与否的过程。
最短路径满足的条件是dis[u]+E[u->v]>=dis[v], 因此，我们可以遍历每条边u->v, 检测权重w是否大于等于dis[v]-dis[u]。
这里面有几个有趣的定理和定义。SPT的验证也引发了FB算法。
定理1: forward-only验证算法对于边的期望验证数量是(1+O(1))nlgn。
定义1：
对于阈值M(任意值)，边u->v为out-pertinent需要满足E[u->v] <= 2*(M - dis[u])，边u->v为in-pertinent需要满足E[u->v] < 2*(dis[v] - M)。而in-pertinent和out-pertinent都称为pertinent。
有了pertinent的定义后，我们可以得到一个有趣的结论，SPT中的任何一条边都必须为in-pertinent或者out-pertinent，并且不能同时满足两者。
这个条件非常有意思，可以简单证明：
假设E[u->v]同时满足两者。那么，将两个不等式相加可得
2 * E[u->v] < 2 * (dis[v] - dis[u]), 即E[u->v] < dis[v]-dis[u]。
显然与SPT成立条件矛盾。
那么假设E[u->v]两者都不满足，即E[u->v]>2*(M-dis[u]), E[u->v]>=2*(dis[v]-M)。将两者相加
2 * E[u->v] > 2 * (dis[v] - dis[u])。同时也与SPT成立条件矛盾。
那么可以得证，SPT中的边满足in-pertinent和out-pertinent。
4. FB算法
在FB算法中M取dis数组的中位数（不得不想到今天CF的C题啊，全是泪水啊）。
由上述有趣的定理。我们可以先通过Spira算出一半结点，然后得到M。然后，由M定界把out-pertinent中的边也加入候选边。
FB算法可以仔细想想M为什么选中位数，大概估计一下这样选择后有多少条会被扫到。
其中队列Q的while循环，如果P为空则min(P)=INF。算法代码如下：
1 /* foward-backward spira */ 2 #include <iostream> 3 #include <string> 4 #include <map> 5 #include <queue> 6 #include <set> 7 #include <stack> 8 #include <vector> 9 #include <deque> 10 #include <algorithm> 11 #include <cstdio> 12 #include <cmath> 13 #include <ctime> 14 #include <cstring> 15 #include <climits> 16 #include <cctype> 17 #include <cassert> 18 #include <functional> 19 #include <iterator> 20 #include <iomanip> 21 using namespace std; 22 //#pragma comment(linker,"/STACK:102400000,1024000") 23 24 #define sti set<int> 25 #define stpii set<pair<int, int> > 26 #define mpii map<int,int> 27 #define vi vector<int> 28 #define pii pair<int,int> 29 #define vpii vector<pair<int,int> > 30 #define rep(i, a, n) for (int i=a;i<n;++i) 31 #define per(i, a, n) for (int i=n-1;i>=a;--i) 32 #define clr clear 33 #define pb push_back 34 #define mp make_pair 35 #define fir first 36 #define sec second 37 #define all(x) (x).begin(),(x).end() 38 #define SZ(x) ((int)(x).size()) 39 #define lson l, mid, rt<<1 40 #define rson mid+1, r, rt<<1|1 41 42 typedef struct node_t { 43 int u, v, w; 44 node_t() {} 45 node_t(int u_, int v_, int w_): 46 u(u_), v(v_), w(w_) {} 47 friend bool operator< (const node_t& a, const node_t& b) { 48 return a.w > b.w; 49 } 50 } node_t; 51 52 const int INF = 0x1f1f1f1f; 53 const int maxv = 205; 54 const int maxe = maxv * maxv + 5; 55 node_t ND[maxe]; 56 int dis[maxv]; 57 bool inS[maxv], active[maxv], out[maxv]; 58 int U[maxe], V[maxe], W[maxe]; 59 int RV[maxe], RW[maxe]; 60 int inxt[maxe], onxt[maxe], rnxt[maxe]; 61 int ihead[maxv], ohead[maxv], rhead[maxv]; 62 int nv, ne = 0, rne = 0, MID; 63 64 65 void addEdge(int u, int v, int w) { 66 U[ne] = u; 67 V[ne] = v; 68 W[ne] = w; 69 70 onxt[ne] = ohead[u]; 71 ohead[u] = ne; 72 73 inxt[ne] = ihead[v]; 74 ihead[v] = ne; 75 ++ne; 76 } 77 78 void addReqEdge(int u, int v, int w) { 79 RV[rne] = v; 80 RW[rne] = w; 81 rnxt[rne] = rhead[u]; 82 rhead[u] = rne++; 83 } 84 85 void forward(priority_queue<node_t>& Q, int u) { 86 int v, w; 87 88 if (out[u]) { 89 int& k = ohead[u]; 90 91 if (k == -1) { 92 out[u] = false; 93 } else { 94 v = V[k]; 95 w = W[k]; 96 if (w > 2*(MID - dis[u])) 97 out[u] = false; 98 k = onxt[k]; 99 active[u] = true; 100 } 101 } 102 103 if (!out[u]) { 104 int& k = rhead[u]; 105 106 if (k == -1) { 107 active[u] = false; 108 } else { 109 v = RV[k]; 110 w = RW[k]; 111 k = rnxt[k]; 112 active[u] = true; 113 } 114 } 115 116 if (active[u]) { 117 Q.push(node_t(u, v, dis[u]+w)); 118 } 119 } 120 121 void backward(priority_queue<node_t>& Q, int v) { 122 int& k = ihead[v]; 123 124 if (k != -1) { 125 Q.push(node_t(k, v, W[k])); 126 k = inxt[k]; 127 } 128 } 129 130 void request(priority_queue<node_t>& Q, int k) { 131 int u = U[k], v = V[k], w = W[k]; 132 133 if (w <= 2*(MID - dis[u])) 134 return ; 135 136 // append(Req[u], v) 137 addReqEdge(u, v, w); 138 139 if (inS[u] && !active[u]) { 140 forward(Q, u); 141 } 142 } 143 144 void sssp(int s = 1) { 145 priority_queue<node_t> P, Q; 146 node_t nd; 147 int u, v, k, n = 1; 148 int mnp, mnq; 149 int nth = (nv+1) >> 1; 150 151 MID = INF; 152 memset(dis, INF, sizeof(dis)); 153 memset(inS, false, sizeof(inS)); 154 memset(out, true, sizeof(out)); 155 dis[s] = 0; 156 inS[s] = true; 157 forward(P, s); 158 159 while (n<nv && !P.empty()) { 160 nd = P.top(); 161 P.pop(); 162 forward(P, nd.u); 163 if (!inS[nd.v]) { 164 ++n; 165 inS[nd.v] = true; 166 dis[nd.v] = nd.w; 167 forward(P, nd.v); 168 if (n == nth) { 169 MID = nd.w; 170 rep(i, 1, nv+1) 171 if (!inS[i]) 172 backward(Q, i); 173 } 174 } 175 176 while (!Q.empty()) { 177 mnp = P.empty() ? INF:P.top().w; 178 mnq = Q.top().w; 179 if (mnq >= 2*(mnp - MID)) 180 break; 181 nd = Q.top(); 182 k = nd.u; 183 Q.pop(); 184 if (!inS[nd.v]) { 185 backward(Q, nd.v); 186 request(P, k); 187 } 188 } 189 } 190 } 191 192 void input() { 193 int m; 194 195 scanf("%d %d", &nv, &m); 196 rep(i, 0, m) 197 scanf("%d %d %d", &ND[i].u, &ND[i].v, &ND[i].w); 198 199 // init 200 ne = rne = 0; 201 memset(ihead, -1, sizeof(ihead)); 202 memset(ohead, -1, sizeof(ohead)); 203 memset(rhead, -1, sizeof(rhead)); 204 sort(ND, ND+m); 205 rep(i, 0, m) 206 addEdge(ND[i].u, ND[i].v, ND[i].w); 207 } 208 209 void solve() { 210 sssp(); 211 212 rep(i, 1, nv+1) 213 printf("%d: %d\n", i, dis[i]); 214 } 215 216 int main() { 217 ios::sync_with_stdio(false); 218 #ifndef ONLINE_JUDGE 219 freopen("data.in", "r", stdin); 220 freopen("data.out", "w", stdout); 221 #endif 222 223 input(); 224 solve(); 225 226 #ifndef ONLINE_JUDGE 227 printf("time = %d.\n", (int)clock()); 228 #endif 229 230 return 0; 231 }
FB要比Spira提高了很多效率。但是，不得不说的是需要维护P、Q两个优先级队列，而且，实际上空间大小是Spira的几倍。
相当于重新建了个request的图。
FB和Spira也比SPFA会快一些。
FB算法后面的定理，我觉得没什么意思。关于每条SPT的边为in-pertinent或out-pertinent这个定理，我觉得挺精彩的。
最后，思考一下为什么Spira没人用，而都选择Dijkstra(SPFA)。
我认为主要还是预处理的条件，其实Spira这个算法还是挺容易写的。而且Spira很适合解决APSP问题（要比floyd好）。
【SSSP】A forward-backward single-source paths algorithm
标签：
原文地址：http://www.cnblogs.com/bombe1013/p/4909700.html
