出来ること

頂点数$V$,辺数$E$の有向グラフが与えられる。
頂点$s$から頂点$t$へ向かう$walk$の中で小さい順に$K$個の経路長を$O((V+E)logV+KlogK)$で列挙する。

walkとは

sから有向辺をたどってtにたどり着く様な経路のこと。
pathは同一頂点を2回以上通ってはならないが、walkは何回通ってもいい。

前提知識

heapとダイクストラ法についての基礎的な事

1. 方針

(1)終点から逆辺にダイクストラを流し、ポテンシャル(終点からの距離)を決める

(2)終点からの最短経路木$T$を求める

(3)ある点$v$から$T$に含まれる辺を通って行ける頂点集合を始点にする、$T$に含まれない辺の集合からなるheap$H_g(v)$を作る

(4)$H_g(v)$の重みを辺に移す(陽に移す必要はない)

(5)$H_g(v)$同士を繋げる辺を作る(陽に繋げる必要はない)

(6)heapから最良優先探索でtop-kを取得

2. Persistent Meldable Heap(Leftist Heap)

$Eppstein's Algorithm$において$H_g(v)$を作るにはpriority_queueでは機能が足りないです。
具体的には以下の機能が求められます

(1)meld:heapAの要素をheapBに挿入する、heapAはそのまま$O(log(要素数))$
(2)insert:heapAに要素xを追加する$O(log(要素数))$

これが出来るのがPersistent Meldable Heapです！！！
今回は作りやすいLeftist Heapを使っていきます！！

それではまず普通のLeftist Heapから作っていきましょう。
hosさんのbunny loves optimized gamesも詳しいのでそちらも良ければ参照してください。

meld

左の子の自身の部分木内の葉からの最短距離が右の子の自身の部分木内の葉からの最短距離以上で有る状態を保ちます。
自身の部分木内の葉からの最短距離を以下sと呼びます。
すると、$s=k$という事は$2^{k-1}$個以上の子(自身も含む)を持つということなのでkは$O(log(要素数))$で抑えられ、右へ右へ降りていった時の葉までの距離は$O(log(要素数))$で抑えられます。

よって次のアルゴリズムが成り立ちます。

heap meld(heap A,heap B){
    if( (Aの最小値) > (Bの最小値) )swap(A,B)
    (Aの右の子) = meld( (Aの右の子) , B )
    if( (Aの左の子)のs < (Aの右の子)のs )swap( (Aの左の子) , (Aの右の子) )
    (Aのs) <- (Aの右の子のs) + 1
    return A
}

これは普通のheapに(Aの左の子)のs < (Aの右の子)のsが逆なら左右を入れ替える操作を付け加えただけです。

insert

すでにheapを作った事がある人はおわかりかと思いますが。
要素数1のheapを作ってmeldするだけです。

永続化

実は木構造の場合簡単にできます。
つくってなぐろ (永続配列/永続Union-Find木/永続セグメント木の作り方と意義、具体例)
meldの部分で、違う木が同じノードを参照しない事を保証すればいいので
meldの部分を

heap meld(heap A,heap B){
    A' <- Aのコピー
    B' <- Bのコピー
    if( (A'の最小値) > (B'の最小値) )swap(A',B')
    (A'の右の子) = meld( (A'の右の子) , B' )
    if( (A'の左の子)のs < (A'の右の子)のs )swap( (A'の左の子) , (A'の右の子) )
    (A'のs) <- (A'の右の子のs) + 1
    return A'
}

に変えるだけでOKです

C++実装例

persistent_leftist_heap.cpp

template<typename T>
struct heap{
    struct node{
        node* ch[2]={0,0};
        int s;
        T val;
        node(T val):s(1),val(val){}
    };
    using np=node*;
    np root=0;
    heap(np t=0):root(t){}
    np meld(np a,np b){
        if(!b)return a?new node(*a):0;
        if(!a)return b?new node(*b):0;
        a=new node(*a);b=new node(*b);
        if(a->val>b->val)swap(a,b);
        a->ch[1]=meld(a->ch[1],b);
        if(!a->ch[0]||a->ch[0]->s<a->ch[1]->s)swap(a->ch[0],a->ch[1]);
        a->s=(a->ch[1]?a->ch[1]->s:0)+1;
        return a;
    }
    heap meld(heap b){
        return heap(meld(root,b.root));
    }
    heap insert(T x){
        return heap(meld(root,new node(x)));
    }
    heap pop(){
        return heap(meld(root->ch[0],root->ch[1]));
    }
    T top(){
        return root?root->val:T();
    }
};

ポテンシャルの概念説明

ダイクストラとポテンシャルのはなしも詳しいです

頂点集合を定義域に持つポテンシャル関数を$p(x)$として(自由に決められる)、辺$(u \rightarrow v)$のポテンシャルは$d'(u\rightarrow v)=d(u\rightarrow v)-p(u)+p(v)$と定義する。
ここで、$d(u\rightarrow v)$は辺の長さである

$s\rightarrow a\rightarrow b\rightarrow c\rightarrow g$と表されるwalkの経路長と辺のポテンシャルの和を比較してみよう

経路長
$D=d(s\rightarrow a)+d(a\rightarrow b)+d(b\rightarrow c)+d(c\rightarrow g)$
辺のポテンシャルの和:
$D'=d'(s\rightarrow a)+d'(a\rightarrow b)+d'(b\rightarrow c)+d'(c\rightarrow g)$
$=d(s\rightarrow a)-p(s)+p(a)+d(a\rightarrow b)-p(a)+p(b)+d(b\rightarrow c)-p(b)+p(c)+d(c\rightarrow g)-p(c)+p(g)$
$=D-p(s)+p(t)$
つまり、和の中抜けによって綺麗に表せるのだ！！！！
また、始点$s$と終点$t$が固定ならば$-p(s)+p(t)$は定数となるので、経路長を求める代わりに辺のポテンシャルの和の小さい順にk個求めれば良い事がわかった

アルゴリズム

(https://qiita.com/hotman78/items/9c643feae1de087e6fc5)
Finding the k Shortest Paths (SICOMP'98)-(iwi) 備忘録
に赤コーダー特有の完結で要点を抑えたお気持ちが書いてある

Eppstein's Algorithm (Find the K shortest paths) 解説と実装 (Python)
に図を踏まえたわかりやすい解説が書いてある。
表記はここに合わせているので両方を見ながら考えて見て欲しい

ポテンシャル関数の設定

方針の欄で書いたが、今回の使用するポテンシャル関数$p(x)$は終点を$t$とした時、$p(x)=d(x,t)$。即ち、$x$から$t$の最短経路の距離である。
これは、$t$から逆辺に対し、ダイクストラすることで求まる。

最短経路木の作成

終点$t$から逆辺に対し、ダイクストラして出来た最短経路木を$T$とする。
ポテンシャルの定義より、$T$に含まれる辺のポテンシャルは0だが、辺のポテンシャルが0でも$T$に含まれるとは限らない(最短経路は複数ある可能性があるので)

heap H_g(v)を作成

ある点$v$から$T$に含まれる辺を通って行ける頂点集合を始点にする、$T$に含まれない辺の集合を$S(v)$とする。
$S(v)$からなるheapを$H_g(v)$と定義する。ただし、要素の大きさは辺のポテンシャルであり、辺の終点をノードに保持しておくと、後の実装が楽である。

この時、$S(v)=S(Tにおけるvの親p)+(Tに含まれない、vを始点とする辺)$
となるので、
最短経路木の根(終点)からDFS/BFSをしながら、
$H_g(v)=meld(H_g(v),H_g(p))$としてから、
(Tに含まれない、vを始点とする辺)を$H_g(v)$に追加すれば良い。
ここで、根からのDFS/BFSなので、$v$について$H_g(v)$を作っている時には、$H_g(p)$は出来上がっている。

heapの重みを頂点から辺に移す

陽にやる必要はない。
heap上で、辺の長さを$(子のポテンシャル)-(親ののポテンシャル)$とすると、$(vのポテンシャル)$は$(rootのポテンシャル)+(rootからvまでの距離)$で表される

heap同士をくっつける

陽にやる必要はない。
全ての経路はこう表せる
$(始点s)$
$\rightarrow (Tに含まれる0個以上の辺) \rightarrow (Tに含まれない1個の辺)$ $\rightarrow (Tに含まれる0個以上の辺)\rightarrow (Tに含まれない1個の辺)$
$\rightarrow ...$
$\rightarrow (Tに含まれる0個以上の辺)$
$\rightarrow (終点t)$

この内、$(Tに含まれない1個の辺)$のみ重みを持つので、heapの各ノードについて、$H_g((辺の終点))$に3本目の辺 ($H_g(辺の終点))$の$root$のポテンシャルを辺の長さとする)を伸ばせば、(最短経路以外の)全経路のheapが出来たことになる。

始点の追加

これも陽にやる必要は無いが、陽にやったほうが実装が楽。
最短経路を経路に付け加えるため、超頂点$p$から始点$s$に長さ最短経路長の辺を貼る。
実際には超頂点$p$に対して、$H_g(p)$を作って、そこに辺の終点$s$,ポテンシャル$(最短経路長$の頂点を挿入する

heapの説明

ヒープの1,2個目の辺を使うという事は、使う辺を入れ替えるという事である。
よって、使う辺のポテンシャルの差分を更新すればいい。
ヒープの3個目の辺を使うということは、現在候補となっている辺の使用を決定し、辺の移動先の最も短い辺$H_g(辺の終点))$の$root$を候補にする事である。
よって$H_g(辺の終点))$の$root$のポテンシャルを差分として更新すればいい。

最良優先探索

$1\sim k-1$番目に短いwalkのどれかと$k$番目に短いwalkはheap上で隣接している。(辺の重みは全て正なので直感的にわかる)
heapは木(DAG)である為、ダイクストラの様に一度通った点は二度と更新される事はない。
よって、priority_queueを使いダイクストラもどきをして、$i$個目に取り出された頂点の始点との距離が$i$番目に小さい経路長である。

勝利！！！！！！

参考文献

hosさんのheap説明
 niuezさんのポテンシャル説明
 つくってなぐろ (永続配列/永続Union-Find木/永続セグメント木の作り方と意義、具体例)(永続化のお気持ちを昔書いたやつ)
Finding the k Shortest Paths (SICOMP'98)-(iwi) 備忘録
 Eppstein's Algorithm (Find the K shortest paths) 解説と実装 (Python)
元論文
他

謝辞

きっかけをくださったnullさん、他色々と助けてくださる皆様ありがとうございます。

verify

K-Shortest Walk -Library Checker(yosupo judge)
314. Shortest Paths -codeforces
回答例(wandbox)
No.1069 電柱 / Pole (Hard) -yukicoder
walkではなくpathなので、O((V+E)logV+KlogK)では恐らく解けない...

Eppstein を path にする論文もありましたが、実用的には高速 (≓ 一般ケースは難しい) みたいなことが書いてありました
— 熨斗袋 (@noshi91) June 5, 2020

C++実装例

eppstein.cpp

#include<bits/stdc++.h>
using namespace::std;
constexpr int INF=1LL<<30;
template<typename T>
struct heap{
    struct node{
        node* ch[2]={0,0};
        int s;
        T val;
        int from,to;
        node(T val,int from,int to):s(1),val(val),from(from),to(to){}
    };
    using np=node*;
    np root=0;
    heap(np t=0):root(t){}
    np meld(np a,np b){
        if(!b)return a?new node(*a):0;
        if(!a)return b?new node(*b):0;
        a=new node(*a);b=new node(*b);
        if(a->val>b->val)swap(a,b);
        a->ch[1]=meld(a->ch[1],b);
        if(!a->ch[0]||a->ch[0]->s<a->ch[1]->s)swap(a->ch[0],a->ch[1]);
        a->s=(a->ch[1]?a->ch[1]->s:0)+1;
        return a;
    }
    heap meld(heap b){
        return heap(meld(root,b.root));
    }
    heap insert(T x,int from,int to){
        return heap(meld(root,new node(x,from,to)));
    }
    heap pop(){
        return heap(meld(root->ch[0],root->ch[1]));
    }
    T top(){
        return root?root->val:T(-1);
    }
    bool empty(){
        return !root;
    }
};

int main(){
    int n,m,k;
    cin>>n>>m>>k;
    int x,y;
    cin>>x>>y;
    x--;y--;
    vector<vector<tuple<int,int,int>>>g(n);
    vector<vector<tuple<int,int,int>>>g2(n);
    for(int i=0;i<m;++i){
        int s,t,c;
        cin>>s>>t>>c;
        s--;t--;
        g[t].emplace_back(s,c,i);
        g2[s].emplace_back(t,c,i);
    }
    vector<int>p(n,INF);
    vector<int>par(n,-1);
    vector<int>idx(n,-1);
    vector<vector<int>>ch(n);
    priority_queue<pair<int,int>,vector<pair<int,int>>,greater<pair<int,int>>>que;
    que.emplace(0,y);
    p[y]=0;
    while(!que.empty()){
        auto [t,n]=que.top();
        que.pop();
        for(auto [e,c,id]:g[n]){
            if(p[e]>t+c){
                p[e]=t+c;
                par[e]=n;
                idx[e]=id;
                que.emplace(p[e],e);
            }
        }
    }
    for(int i=0;i<n;++i){
        if(par[i]!=-1)ch[par[i]].push_back(i);
    }
    vector<heap<int>>h(n);
    for(int i=0;i<n;++i)h[i]=heap<int>();
    queue<int>qq;
    qq.emplace(y);
    while(!qq.empty()){
        auto n=qq.front();
        qq.pop();
        if(par[n]!=-1)h[n]=h[n].meld(h[par[n]]);
        for(auto [e,c,id]:g2[n]){
            if(id!=idx[n])h[n]=h[n].insert(c-p[n]+p[e],n,e);
        }
        for(auto e:ch[n]){
            qq.emplace(e);
        }
    }
    auto comp=[](auto s,auto t){return get<0>(s)>get<0>(t);};
    priority_queue<tuple<int,heap<int>::np>,vector<tuple<int,heap<int>::np>>,decltype(comp)>qqq(comp);
    heap<int>s;
    s=s.insert(p[x],-1,x);
    qqq.emplace(p[x],s.root);
    vector<int>ans;
    while(!qqq.empty()){
        auto [t,now]=qqq.top();
        qqq.pop();
        if(t>=INF)break;
        ans.push_back(t);
        if(int(ans.size())==k)break;
        if(now->ch[0]){
            qqq.emplace(t+now->ch[0]->val-now->val,now->ch[0]);
        }
        if(now->ch[1]){
            qqq.emplace(t+now->ch[1]->val-now->val,now->ch[1]);
        }
        if(h[now->to].root){
            qqq.emplace(t+h[now->to].root->val,h[now->to].root);
        }
    }
    for(int i=0;i<k;++i){
        if(i<(int)ans.size())cout<<ans[i]<<endl;
        else cout<<"NO"<<endl;
    }
}