【SSSP】A forward-backward single-source paths algorithm

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算法了。算法源代码如下

/* Spira */
#include <iostream>
#include <string>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <vector>
#include <deque>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cstring>
#include <climits>
#include <cctype>
#include <cassert>
#include <functional>
#include <iterator>
#include <iomanip>
using namespace std;
//#pragma comment(linker,"/STACK:102400000,1024000")

#define sti                set<int>
#define stpii            set<pair<int, int> >
#define mpii            map<int,int>
#define vi                vector<int>
#define pii                pair<int,int>
#define vpii            vector<pair<int,int> >
#define rep(i, a, n)     for (int i=a;i<n;++i)
#define per(i, a, n)     for (int i=n-1;i>=a;--i)
#define clr                clear
#define pb                 push_back
#define mp                 make_pair
#define fir                first
#define sec                second
#define all(x)             (x).begin(),(x).end()
#define SZ(x)             ((int)(x).size())
#define lson            l, mid, rt<<1
#define rson            mid+1, r, rt<<1|1

typedef struct node_t {
    int u, v, w;
    node_t() {}
    node_t(int u_, int v_, int w_):
        u(u_), v(v_), w(w_) {}
    friend bool operator< (const node_t& a, const node_t& b) {
        return a.w > b.w;
    }
} node_t;

const int INF = 0x1f1f1f1f;
const int maxv = 205;
const int maxe = maxv * maxv + 5;
node_t ND[maxe];
int dis[maxv];
bool inS[maxv];
int nv, ne;
int V[maxe], W[maxe], nxt[maxe];
int head[maxv], head_[maxv];

void addEdge(int u, int v, int w) {
    V[ne] = v;
    W[ne] = w;
    nxt[ne] = head_[u];
    head_[u] = ne++;
}

void forward(priority_queue<node_t>& Q, int u) {
    int& k = head[u];

    if (k != -1) {
        Q.push(node_t(u, V[k], dis[u]+W[k]));
        k = nxt[k];
    }
}

void spira(int s = 1) {
    priority_queue<node_t> Q;
    node_t nd;

    memset(inS, false, sizeof(inS));
    memset(dis, INF, sizeof(dis));
    memcpy(head, head_, sizeof(head));
    inS[s] = true;
    dis[s] = 0;
    forward(Q, s);

    while (!Q.empty()) {
        nd = Q.top();
        Q.pop();
        forward(Q, nd.u);
        if (!inS[nd.v]) {
            inS[nd.v] = true;
            dis[nd.v] = nd.w;
            forward(Q, nd.v);
        }
    }
}

void input() {
    int m;

    scanf("%d %d", &nv, &m);
    rep(i, 0, m)
        scanf("%d %d %d", &ND[i].u, &ND[i].v, &ND[i].w);

    ne = 0;
    memset(head_, -1, sizeof(head_));
    
    // pre sort
    sort(ND, ND+m);
    rep(i, 0, m)
        addEdge(ND[i].u, ND[i].v, ND[i].w);
}

void solve() {
    spira();

    rep(i, 1, nv+1)
        printf("%d: %d
", i, dis[i]);
}

int main() {
    ios::sync_with_stdio(false);
    #ifndef ONLINE_JUDGE
        freopen("data.in", "r", stdin);
        freopen("data.out", "w", stdout);
    #endif

    input();
    solve();

    #ifndef ONLINE_JUDGE
        printf("time = %d.
", (int)clock());
    #endif

    return 0;
}

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。算法代码如下：

/* foward-backward spira */
#include <iostream>
#include <string>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <vector>
#include <deque>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <ctime>
#include <cstring>
#include <climits>
#include <cctype>
#include <cassert>
#include <functional>
#include <iterator>
#include <iomanip>
using namespace std;
//#pragma comment(linker,"/STACK:102400000,1024000")

#define sti                set<int>
#define stpii            set<pair<int, int> >
#define mpii            map<int,int>
#define vi                vector<int>
#define pii                pair<int,int>
#define vpii            vector<pair<int,int> >
#define rep(i, a, n)     for (int i=a;i<n;++i)
#define per(i, a, n)     for (int i=n-1;i>=a;--i)
#define clr                clear
#define pb                 push_back
#define mp                 make_pair
#define fir                first
#define sec                second
#define all(x)             (x).begin(),(x).end()
#define SZ(x)             ((int)(x).size())
#define lson            l, mid, rt<<1
#define rson            mid+1, r, rt<<1|1

typedef struct node_t {
    int u, v, w;
    node_t() {}
    node_t(int u_, int v_, int w_):
        u(u_), v(v_), w(w_) {}
    friend bool operator< (const node_t& a, const node_t& b) {
        return a.w > b.w;
    }
} node_t;

const int INF = 0x1f1f1f1f;
const int maxv = 205;
const int maxe = maxv * maxv + 5;
node_t ND[maxe];
int dis[maxv];
bool inS[maxv], active[maxv], out[maxv];
int U[maxe], V[maxe], W[maxe];
int RV[maxe], RW[maxe];
int inxt[maxe], onxt[maxe], rnxt[maxe];
int ihead[maxv], ohead[maxv], rhead[maxv];
int nv, ne = 0, rne = 0, MID;


void addEdge(int u, int v, int w) {
    U[ne] = u;
    V[ne] = v;
    W[ne] = w;
    
    onxt[ne] = ohead[u];
    ohead[u] = ne;
    
    inxt[ne] = ihead[v];
    ihead[v] = ne;
    ++ne;
}

void addReqEdge(int u, int v, int w) {
    RV[rne] = v;
    RW[rne] = w;
    rnxt[rne] = rhead[u];
    rhead[u] = rne++;
}

void forward(priority_queue<node_t>& Q, int u) {
    int v, w;
    
    if (out[u]) {
        int& k = ohead[u];
        
        if (k == -1) {
            out[u] = false;
        } else {
            v = V[k];
            w = W[k];
            if (w > 2*(MID - dis[u]))
                out[u] = false;
            k = onxt[k];
            active[u] = true;
        }
    }
    
    if (!out[u]) {
        int& k = rhead[u];
        
        if (k == -1) {
            active[u] = false;
        } else {
            v = RV[k];
            w = RW[k];
            k = rnxt[k];
            active[u] = true;
        }
    }
    
    if (active[u]) {
        Q.push(node_t(u, v, dis[u]+w));
    }
}

void backward(priority_queue<node_t>& Q, int v) {
    int& k = ihead[v];
    
    if (k != -1) {
        Q.push(node_t(k, v, W[k]));
        k = inxt[k];
    }
}

void request(priority_queue<node_t>& Q, int k) {
    int u = U[k], v = V[k], w = W[k];
    
    if (w <= 2*(MID - dis[u]))
        return ;
    
    // append(Req[u], v)
    addReqEdge(u, v, w);
    
    if (inS[u] && !active[u]) {
        forward(Q, u);
    }
}

void sssp(int s = 1) {
    priority_queue<node_t> P, Q;
    node_t nd;
    int u, v, k, n = 1;
    int mnp, mnq;
    int nth = (nv+1) >> 1;
    
    MID = INF;
    memset(dis, INF, sizeof(dis));
    memset(inS, false, sizeof(inS));
    memset(out, true, sizeof(out));
    dis[s] = 0;
    inS[s] = true;
    forward(P, s);
    
    while (n<nv && !P.empty()) {
        nd = P.top();
        P.pop();
        forward(P, nd.u);
        if (!inS[nd.v]) {
            ++n;
            inS[nd.v] = true;
            dis[nd.v] = nd.w;
            forward(P, nd.v);
            if (n == nth) {
                MID = nd.w;
                rep(i, 1, nv+1)
                    if (!inS[i])
                        backward(Q, i);
            }
        }
        
        while (!Q.empty()) {
            mnp = P.empty() ? INF:P.top().w;
            mnq = Q.top().w;
            if (mnq >= 2*(mnp - MID))
                break;
            nd = Q.top();
            k = nd.u;
            Q.pop();
            if (!inS[nd.v]) {
                backward(Q, nd.v);
                request(P, k);
            }
        }
    }
}

void input() {
    int m;
    
    scanf("%d %d", &nv, &m);
    rep(i, 0, m)
        scanf("%d %d %d", &ND[i].u, &ND[i].v, &ND[i].w);
        
    // init
    ne = rne = 0;
    memset(ihead, -1, sizeof(ihead));
    memset(ohead, -1, sizeof(ohead));
    memset(rhead, -1, sizeof(rhead));
    sort(ND, ND+m);
    rep(i, 0, m)
        addEdge(ND[i].u, ND[i].v, ND[i].w);
}

void solve() {
    sssp();
    
    rep(i, 1, nv+1)
        printf("%d: %d
", i, dis[i]);
}

int main() {
    ios::sync_with_stdio(false);
    #ifndef ONLINE_JUDGE
        freopen("data.in", "r", stdin);
        freopen("data.out", "w", stdout);
    #endif
    
    input();
    solve();
    
    #ifndef ONLINE_JUDGE
        printf("time = %d.
", (int)clock());
    #endif
    
    return 0;
}

FB要比Spira提高了很多效率。但是，不得不说的是需要维护P、Q两个优先级队列，而且，实际上空间大小是Spira的几倍。
相当于重新建了个request的图。
FB和Spira也比SPFA会快一些。
FB算法后面的定理，我觉得没什么意思。关于每条SPT的边为in-pertinent或out-pertinent这个定理，我觉得挺精彩的。

最后，思考一下为什么Spira没人用，而都选择Dijkstra(SPFA)。
我认为主要还是预处理的条件，其实Spira这个算法还是挺容易写的。而且Spira很适合解决APSP问题（要比floyd好）。