Day 124-Q4 — Eppstein's K-Shortest Paths Algorithm（サイドトラック辺とpersistent leftist heapによるk最短ウォーク列挙）

問題

$N$ 頂点 $M$ 辺の有向グラフ（辺重みは正の整数）と、始点 $s$・終点 $t$・整数 $K$ が与えられる。$s$ から $t$ へのウォーク（頂点や辺の重複を許す経路）のうち、コストが小さい方から $K$ 番目のコストを出力せよ。存在しない場合（$t$ にそもそも到達できない、または到達可能なウォークの総数が $K$ 未満）は -1 を出力せよ。

入力形式

N M K
s t
u_1 v_1 w_1
...
u_M v_M w_M

制約

入出力例

4 5 3
1 4
1 2 1
2 4 1
1 3 2
3 4 3
1 4 10

10

概念図

ヒント（段階的開示）

# leftist heap のマージ（persistent = 経路コピー）
def merge(a, b):
    if a is None: return b
    if b is None: return a
    if a.cost > b.cost: a, b = b, a
    node = LNode(a.edge, a.cost)
    node.l = a.l
    node.r = merge(a.r, b)
    if node.l is None or (node.r and node.l.npl < node.r.npl):
        node.l, node.r = node.r, node.l
    node.npl = (node.r.npl + 1) if node.r else 0
    return node

模範解答 (Python)

📋 コピーimport sys, heapq

def solve():
    data = sys.stdin.read().split()
    idx = 0
    N = int(data[idx]); idx += 1
    M = int(data[idx]); idx += 1
    K = int(data[idx]); idx += 1
    s = int(data[idx]) - 1; idx += 1
    t = int(data[idx]) - 1; idx += 1
    edges = []
    rev = [[] for _ in range(N)]
    for _ in range(M):
        u = int(data[idx]) - 1; idx += 1
        v = int(data[idx]) - 1; idx += 1
        w = int(data[idx]); idx += 1
        edges.append((u, v, w))
        rev[v].append((u, w))

    INF = float('inf')
    dist = [INF] * N
    dist[t] = 0
    succ = [-1] * N
    pq = [(0, t)]
    while pq:
        d, u = heapq.heappop(pq)
        if d > dist[u]:
            continue
        for v, w in rev[u]:  # 元グラフの辺 v -> u
            nd = d + w
            if nd < dist[v]:
                dist[v] = nd
                succ[v] = u
                heapq.heappush(pq, (nd, v))

    if dist[s] == INF:
        print(-1)
        return

    class LNode:
        __slots__ = ['edge', 'cost', 'l', 'r', 'npl']
        def __init__(self, edge, cost):
            self.edge = edge
            self.cost = cost
            self.l = None
            self.r = None
            self.npl = 0

    def merge(a, b):
        if a is None:
            return b
        if b is None:
            return a
        if a.cost > b.cost:
            a, b = b, a
        node = LNode(a.edge, a.cost)
        node.l = a.l
        node.r = merge(a.r, b)
        if node.l is None or (node.r is not None and node.l.npl < node.r.npl):
            node.l, node.r = node.r, node.l
        node.npl = (node.r.npl + 1) if node.r is not None else 0
        return node

    # 各頂点の出辺をグループ化し、ツリー辺(1本だけ)を除いてサイドトラック辺集合を作る
    out_edges = [[] for _ in range(N)]
    for (u, v, w) in edges:
        if dist[u] == INF or dist[v] == INF:
            continue
        out_edges[u].append((u, v, w))

    H_out = [None] * N
    for v in range(N):
        skip_tree_edge = True  # ツリー辺を1本だけ除外する
        h = None
        for (u, vv, w) in out_edges[v]:
            if skip_tree_edge and vv == succ[v] and w == dist[v] - dist[vv]:
                skip_tree_edge = False
                continue
            delta = w + dist[vv] - dist[v]
            h = merge(h, LNode((u, vv, w), delta))
        H_out[v] = h

    order = sorted((v for v in range(N) if dist[v] < INF), key=lambda v: dist[v])
    D = [None] * N
    for v in order:
        if v == t:
            D[v] = H_out[v]
        else:
            D[v] = merge(H_out[v], D[succ[v]])

    # K番目のウォークコストを生成
    results = [dist[s]]
    counter = 0
    pq2 = []
    if D[s] is not None:
        heapq.heappush(pq2, (dist[s] + D[s].cost, counter, D[s]))
        counter += 1

    while pq2 and len(results) < K:
        cost, _, node = heapq.heappop(pq2)
        results.append(cost)
        if node.l is not None:
            heapq.heappush(pq2, (cost - node.cost + node.l.cost, counter, node.l))
            counter += 1
        if node.r is not None:
            heapq.heappush(pq2, (cost - node.cost + node.r.cost, counter, node.r))
            counter += 1
        v_land = node.edge[1]
        if D[v_land] is not None:
            heapq.heappush(pq2, (cost + D[v_land].cost, counter, D[v_land]))
            counter += 1

    if len(results) < K:
        print(-1)
    else:
        print(results[K - 1])

sys.setrecursionlimit(50000)
solve()

Step-by-Step 解説

よくあるミス

次のステップ

• 発展: 単純パスに制限したK最短路が必要な場合はYenのアルゴリズム（$O(K \cdot N \cdot (\text{Dijkstra}))$）を使う。負辺がある場合はBellman-Ford系の前処理でポテンシャルを作ってから同様の手法を適用できる。
• 次回予告: Berlekamp's Q-Matrix Algorithm（有限体上の多項式因数分解）

自己評価

自分の回答

気づき・メモ