淺談k短路算法

An Old but Classic Problem

　　給定一個$n$個點，$m$條邊的帶正權有向圖。給定$s$和$t$，詢問$s$到$t$的所有權和為正路徑中，第$k$短的長度。

Notice

　　定義兩條路徑不同，當且僅當它們的邊集中存在一條邊，使得它只在其中的一條路徑上。

Solution#1 Shortest Path & A*

　　對於Dijstra算法，有一個結論就是，當一個點第$k$次出隊的時候，此時路徑長度就是$s$到它的第$k$短路。

　　那為什么還要A*呢？我試了試，寫了個Dijstra，然后交了一發poj 2449，於是MLE了。。。如果您覺得我的Dijstra太丑了，您可以手寫一個交一發。

　　所以用A*來優化優化狀態數（其實就是玄學優化，需要你的人品還有出題人的素質）。

　　首先建反圖，跑一次最短路算法，得到每個點到$t$的最短路的距離。

　　然后用當前走的距離加上到終點的最短路的長度作為優先級進行A*。

　　那如何得到答案？

1. 當一個點第k次出隊時，答案是它的優先級
2. 當終點第k次出隊時，答案是它已經走的路程

　　據說A*可以被卡成$O\left(nk\log n \right )$，只是我不會QAQ。

Code

/**
 * poj
 * Problem#2449
 * Accepted
 * Time: 250ms
 * Memory: 9252k
 */
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <cctype>
#include <cmath>
#include <vector>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <bitset>
using namespace std;
typedef bool boolean;

typedef class Edge {
    public:
        int end;
        int next;
        int w;
        
        Edge(int end = 0, int next = -1, int w = 0):end(end), next(next), w(w) {        }
}Edge;

const int N = 1e3, M = 1e5;

typedef class MapManager {
    public:
        int cnt;
        int h[N + 5];
        Edge edge[M + 5];
        
        MapManager() {        }
        MapManager(int n):cnt(-1) {
//            h = new int[(n + 1)];
//            edge = new Edge[(m + 1)];
            memset(h, -1, sizeof(int) * (n + 1));
        }
        
        inline void addEdge(int u, int v, int w) {
            edge[++cnt] = (Edge(v, h[u], w));
//            h[u] = (signed)edge.size() - 1;
            h[u] = cnt; 
        }
        
        inline int start(int node) {    return h[node];        }
        
        Edge& operator [] (int pos) {
            return edge[pos];
        }
}MapManager;
#define m_endpos -1

int n, m;
MapManager g;
MapManager rg;
int s, t, k;
int ds[N + 5];

inline void init() {
    scanf("%d%d", &n, &m);
    memset(g.h, -1, sizeof(int) * (n + 1));
    memset(rg.h, -1, sizeof(int) * (n + 1));
    for(int i = 1, u, v, w; i <= m; i++) {
        scanf("%d%d%d", &u, &v, &w);
        g.addEdge(u, v, w);
        rg.addEdge(v, u, w);
    }
    scanf("%d%d%d", &s, &t, &k);
//    ds = new int[(n + 1)];
}

#define g rg
#define f ds
#define que que1
boolean vis[N + 5];
queue<int> que;
boolean spfa(int s, int t) {
    memset(f, 0x7f, sizeof(int) * (n + 1));
    memset(vis, false, sizeof(boolean) * (n + 1));
    que.push(s);
    f[s] = 0;
    while(!que.empty()) {
        int e = que.front();
        que.pop();
        vis[e] = false;
        for(int i = g.start(e); i != m_endpos; i = g[i].next) {
            int& eu = g[i].end;
//            cout << e << " " << eu << " " << i <<endl;
            if(f[e] + g[i].w < f[eu]) {
                f[eu] = f[e] + g[i].w;
                if(!vis[eu]) {
                    que.push(eu);
                    vis[eu] = true;
                }
            }
        }
    }
    return (f[t] != 0x7f7f7f7f);
}
#undef g
#undef f
#undef que

typedef class Status {
    public:
        int node;
        int dis;
        int priority;
        
        Status(int node = 0, int dis = 0):node(node), dis(dis), priority(h()) {        }
        
        int h() {
            return dis + ds[node];
        }
        
        boolean operator < (Status b) const {
            return priority > b.priority;
        }
}Status;

int label[N + 5];
priority_queue<Status> que;
int bfs(int s, int t) {
    if(s == t)    k++;
//    label = new int[(n + 1)];
    memset(label, 0, sizeof(int) * (n + 1));
    que.push(Status(s, 0));
    while(!que.empty()) {
        Status e = que.top();
        que.pop();
        label[e.node]++;
        if(e.node == t && label[e.node] == k)
            return e.dis;
        for(int i = g.start(e.node); i != m_endpos; i = g[i].next) {
            if(label[g[i].end] < k) 
                que.push(Status(g[i].end, e.dis + g[i].w));
        } 
    }
    return -1;
}

inline void solve() {
    if(!spfa(t, s)) {
        puts("-1");
        return;
    }
    printf("%d", bfs(s, t));
}

int main() {
    init();
    solve();
    return 0;
}

Solution#2  Shortest Path & Protractable Heap

　　考慮建立反圖，然后跑最短路算法得到以$t$為根的最短路徑生成樹。

　　當走了一條非樹邊$(u, v, w)$意味着什么？

 

　　最終的路徑長度就會因此增加$f[v] - f[u] + w$。

　　對於一條路徑，我們依次將它經過的非樹邊記下來，約定得到的序列是這條路徑的非樹邊序列。

　　考慮對於一個合法的非樹邊序列，我們可以找到唯一的一條$s$到$t$的路徑與之對應。

　　因此，$k$短路的長度就等於第$k$小的代價和加上$s$到$t$的最短路的長度。

　　考慮如何來得到一個合法的非樹邊序列。

1. 找到一條起點在當前點$p$到根$t$的路徑上的非樹邊
2. 令p等於這條邊的終點。

　　我們可以通過這樣的方法來得到所有的非樹邊序列。但是我們並不需要所有的非樹邊序列，因此當找到第$x$短路后再進行拓展狀態，然后用優先隊列來維護。

　　但是這樣每次拓展時間復雜度可達$O(m)$，總時間復雜度可以達到$O\left(mk\log \left (mk  \right ) \right )$。

　　令人無法接受。但是其中真正會被用到的狀態十分少。

　　因此可以像UVa 11997那樣進行優化。

　　當一個非樹邊序列出隊時，代價和比它大的才可能有用。

　　因此，考慮一個非樹邊序列出隊時通過下面的方法來進行得到新的序列:

1. 追加操作：假如最后一條非樹邊的終點為$v$，找到一條起點在$v$到$t$的路徑上代價最小的非樹邊追加在當前非樹邊序列后
2. 替換操作：將最后一條非樹邊更換為代價比它大的1條非樹邊。
   例如圖中橙色虛線是被替換掉的非樹邊，紫色是新加入的非樹邊
   

　　考慮用一些可持久化數據結構（如可持久化斜可並堆，可持久化線段樹，可持久化Treap）來維護起點在點$u$到根的路徑上的非樹邊的代價。

　　對於替換操作，

• 如果用的可持久化堆，那么把最后一條非樹邊替換為它所在的堆（你從哪個堆把它拿出來的）中它的左右子節點代表的邊。
• 如果用的可持久化平衡樹，那么把最后一條非樹邊直接替換為它的后繼
• ......

　　這是一個很穩定的算法，時間復雜度$O\left ( n + m\log m + k\log k \right )$。就是常數有點大，sad.....

　　注意一些細節

• 計算代價時需要考慮終點是否可以到達$t$
• 考慮$s = t$時，要求的$k$短路包不包含0
• $k$短路不存在，隊首為空

（注意！不能用斜堆可持久化，斜堆的時間復雜度是均攤的。這里能過應該只是沒有被卡）

Code

/**
 * poj
 * Problem#2449
 * Accepted
 * Time: 438ms
 * Memory: 15196k 
 */
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <vector>
#include <queue>
using namespace std;
typedef bool boolean;

#define pii pair<int, int>
#define fi first
#define sc second

typedef class Node {
    public:
        int val, ed;
        Node *l, *r;
        
        Node()    {        }
        Node(int val, int ed, Node *l, Node *r):val(val), ed(ed), l(l), r(r) {        }
}Node;

#define Limit 1000000

Node pool[Limit];
Node* top = pool;

Node* newnode(int val, int ed) {
    if(top >= pool + Limit)
        return new Node(val, ed, NULL, NULL);
    top->val = val, top->ed = ed, top->l = top->r = NULL;
    return top++;
}

Node* merge(Node* a, Node* b) {
    if (!a)    return b;
    if (!b)    return a;
    if (a->val > b->val)    swap(a, b);
    Node* p = newnode(a->val, a->ed);
    p->l = a->l, p->r = a->r;
    p->r = merge(p->r, b);
    swap(p->l, p->r);
    return p;
}

typedef class Status {
    public:
        int dist;
        Node* p;
        
        Status(int dist = 0, Node* p = NULL):dist(dist), p(p) {        }

        boolean operator < (Status b) const {
            return dist > b.dist;
        }
}Status;

typedef class Edge {
    public:
        int end, next, w;
        
        Edge(int end = 0, int next = 0, int w = 0):end(end), next(next), w(w) {        }
}Edge;

typedef class MapManager {
    public:
        int ce;
        int* h;
        Edge* es;
        
        MapManager() {            }
        MapManager(int n, int m):ce(0) {
            h = new int[(n + 1)];
            es = new Edge[(m + 5)];
            memset(h, 0, sizeof(int) * (n + 1));
        }
        
        void addEdge(int u, int v, int w) {
            es[++ce] = Edge(v, h[u], w);
            h[u] = ce;
        }
        
        Edge& operator [] (int pos) {
            return es[pos];
        }
}MapManager;

int n, m;
int s, t, k;
MapManager g;
MapManager rg;
boolean *vis;
int* f, *lase;

inline void init() {
    scanf("%d%d", &n, &m);
    g = MapManager(n, m);
    rg = MapManager(n, m);
    for (int i = 1, u, v, w; i <= m; i++) {
        scanf("%d%d%d", &u, &v, &w);
        g.addEdge(u, v, w);
        rg.addEdge(v, u, w);
    }
    scanf("%d%d%d", &s, &t, &k);
}

queue<int> que;
void spfa(MapManager& g, int s) {
    vis = new boolean[(n + 1)];
    f = new int[(n + 1)];
    lase = new int[(n + 1)];
    memset(f, 0x7f, sizeof(int) * (n + 1));
    memset(vis, false, sizeof(boolean) * (n + 1));
    que.push(s);
    f[s] = 0, lase[s] = 0;
    while (!que.empty()) {
        int e = que.front();
        que.pop();
        vis[e] = false;
        for (int i = g.h[e]; i; i = g[i].next) {
            int eu = g[i].end, w = g[i].w;
            if (f[e] + w < f[eu]) {
                f[eu] = f[e] + w, lase[eu] = i;
                if (!vis[eu]) {
                    vis[eu] = true;
                    que.push(eu); 
                } 
            }
        }
    }
 }

Node** hs;
inline void rebuild() {
    for (int i = 1; i <= n; i++)
        for (int j = g.h[i]; j; j = g[j].next) {
            int e = g[j].end;
            if (lase[i] != j)
                g[j].w += f[e] - f[i];
        }
    
    hs = new Node*[(n + 1)];
    que.push(t);
    hs[t] = NULL;
    while (!que.empty()) {
        int e = que.front();
        que.pop();
        if (lase[e])
            hs[e] = hs[g[lase[e]].end];
        for (int i = g.h[e]; i; i = g[i].next)
            if (lase[e] != i && f[g[i].end] != 0x7f7f7f7f)
                hs[e] = merge(hs[e], new Node(g[i].w, g[i].end, NULL, NULL));
        for (int i = rg.h[e]; i; i = rg[i].next) {
            int eu = rg[i].end;
            if (lase[eu] == i)
                que.push(eu);
        }
    }
}

inline int kthpath(int k) {
    if (s == t)
        k++;
    if (f[s] == 0x7f7f7f7f)
        return -1;
    if (k == 1)
        return f[s];
    
    priority_queue<Status> q;
    if (!hs[s])
        return -1;
        
    q.push(Status(hs[s]->val, hs[s]));
    while (--k && !q.empty()) {
        Status e = q.top();
        q.pop();
        
        if(k == 1)
            return e.dist + f[s];
        
        int eu = e.p->ed;
        if (hs[eu])
            q.push(Status(e.dist + hs[eu]->val, hs[eu]));
        if (e.p->l)
            q.push(Status(e.dist - e.p->val + e.p->l->val, e.p->l));
        if (e.p->r)
            q.push(Status(e.dist - e.p->val + e.p->r->val, e.p->r));
    }
    return -1;
}

inline void solve() {
    printf("%d\n", kthpath(k));
}

int main() {
    init();
    spfa(rg, t);
    rebuild();
    solve();
    return 0;
}

 

特別鳴謝

　　YJQ

　　ZJC