Bellman-Ford

graph/bellman_ford.hpp · GitHub で見る

共通ページ

負の重みを持つ辺を含む有向グラフの単一始点最短路を求めます。経路と、始点から到達できる負閉路を復元できます。自己辺・多重辺に対応します。

#include "lib/graph/bellman_ford.hpp"

int main() {
	std::vector<DirectedCSR<int>::Edge> edges{ {0, 1, 4}, {0, 2, 5}, {1, 2, -7}, {2, 3, 2} };
	DirectedCSR<int> g(5, edges);
	auto result = bellman_ford(g, 0);
	if (result.has_negative_cycle()) return 0;
	assert(result.distance[2] == -3);
	assert(result.distance[3] == -1);
	assert(!result.reachable(4));
	auto vertices = result.path(3);       // 0, 1, 2, 3
	auto edges_in_path = result.path_edges(3);
	assert(vertices.size() == 4 && edges_in_path.size() == 3);
}

グラフは DirectedCSR<W> または vector<vector<pair<int,W>>> を渡します。各辺は (行き先,重み) です。頂点番号は 0 始まりで、始点は存在する頂点を指定します。重みは符号付き整数型です。

距離型の既定は std::int64_t。bellman_ford<__int128_t>(g,s) のように指定もできます。距離型は重み型以上の桁数を持つ符号付き整数型を使ってください。加算途中も距離型に収まり、有限値がその型の最大値と一致しないことが必要です。範囲外の計算を検出した場合は std::overflow_error を送出します。値を丸めて計算を継続することはありません。

distance[v] が距離型の最大値なら到達不能です。負の距離、特に -1 も正しい答えになり得るため、到達判定には reachable(v) を使います。始点自身の距離は 0、経路は [s] です。到達不能頂点への経路は空配列です。

path_edges(t)[i] は、path(t)[i] の 隣接リスト内の添字 です。入力辺の通し番号ではありません。g[path[i]][edge_indices[i]] が使用した辺になります。多重辺でも実際に選んだ辺を区別できます。parent_edge[v] も、g[parent[v]] 内の添字です。

到達可能な負閉路が見つかると、has_negative_cycle() が真になります。negative_cycle は先頭と末尾が同じ閉路の頂点列、negative_cycle_edges はその各辺の隣接リスト内の添字です。両者から復元した閉路の重みは負です。始点から到達不能な負閉路は報告しません。

負閉路が見つかった場合、distance・parent は計算途中の値です。最短路の取得には使わず、reachable・path・path_edges も呼ばないでください。負閉路から到達する全頂点を負無限大として列挙する機能ではありません。

bellman_ford は通常の緩和走査を最大 2 回行い、更新が止まれば結果を返します。続く場合は到達可能部分を強連結成分に分けます。成分間はトポロジカル順に処理し、成分内に負辺がなければ Dijkstra、負辺があれば Goldberg–Radzik 法を使います。DAG では負辺があっても線形時間です。重みが全て非負と分かっている場合は、通常の dijkstra も使えます。

goldberg_radzik は前処理を省き、Goldberg–Radzik 法だけで解きます。入出力の型と復元機能は同じです。既定の距離型ではどちらも一般の場合の最悪計算量は \(O(nm)\) で、全ての疎グラフをほぼ線形時間で解く保証はありません。

頂点数を \(n\)、辺数を \(m\)、復元する頂点・辺の数を \(L\) とします。

操作・メンバー	内容・返り値	計算量
bellman_ford<Distance = std::int64_t>(g,s)	最短路と到達可能な負閉路の検出・復元	64 bit 以下は最悪 \(O(n+nm)\)、128 bit は \(O(n+nm+m\log m)\)。DAG は \(O(n+m)\)
goldberg_radzik<Distance = std::int64_t>(g,s)	成分分解と前段の走査を省いて同じ結果を求める	最悪 \(O(n+nm)\)
result.has_negative_cycle()	到達可能な負閉路が見つかったか	\(O(1)\)
result.reachable(v)	負閉路がない場合、v に到達できるか	\(O(1)\)
result.distance[v]	始点からの距離。到達不能なら infinity	\(O(1)\)
result.parent[v]	最短路上の直前の頂点。始点は自身、到達不能は -1	\(O(1)\)
result.parent_edge[v]	直前の頂点の隣接リスト内の辺番号。始点・到達不能は size_t(-1)	\(O(1)\)
result.path(t) / result.path_edges(t)	最短路の頂点列 / 各辺の隣接リスト内の添字。到達不能なら空配列	\(O(L)\)
result.negative_cycle / result.negative_cycle_edges	負閉路の頂点列 / 各辺の隣接リスト内の添字。負閉路がなければ空配列	参照 \(O(1)\)、走査 \(O(L)\)
result.source / result.infinity	指定した始点 / 到達不能を表す値	\(O(1)\)

グラフを除く追加空間は、bellman_ford が \(O(n+m)\)、goldberg_radzik が \(O(n)\) です。64 bit 以下の整数の非負辺成分では既存の radix heap を使い、128 bit の距離型では比較ヒープを使います。ヒープには更新前の古い候補も入るため、辺数に比例する領域を使う場合があります。再帰呼び出しを使わず、長い一本道にも対応します。

依存ライブラリ

• Radix Heap
• 有向グラフの CSR
• Shortest Path SCC

検証コード

• test/bellman_ford.test.cpp
• test/bellman_ford_path.test.cpp
• test/difference_constraints.test.cpp

コード

#ifndef BELLMAN_FORD_HPP
#define BELLMAN_FORD_HPP

#include "../components/shortest_path_scc.hpp"
#include "../components/directed_csr.hpp"
#include "../components/RadixHeap.hpp"
#include <cstdint>
#include <functional>
#include <stdexcept>
#include <type_traits>
#include <utility>

template<class Distance>
struct BellmanFordResult {
	std::vector<Distance> distance;
	std::vector<int> parent, negative_cycle;
	std::vector<std::size_t> parent_edge, negative_cycle_edges;
	Distance infinity;
	int source;
	bool has_negative_cycle() const { return !negative_cycle.empty(); }
	bool reachable(int v) const {
		assert(!has_negative_cycle() && 0 <= v && std::size_t(v) < distance.size());
		return distance[v] != infinity;
	}
	std::vector<int> path(int v) const {
		if (!reachable(v)) return {};
		std::vector<int> result;
		for (; v != source; v = parent[v]) result.push_back(v);
		result.push_back(source);
		std::reverse(result.begin(), result.end());
		return result;
	}
	std::vector<std::size_t> path_edges(int v) const {
		if (!reachable(v)) return {};
		std::vector<std::size_t> result;
		for (; v != source; v = parent[v]) result.push_back(parent_edge[v]);
		std::reverse(result.begin(), result.end());
		return result;
	}
};

namespace negative_path_detail {
template<class Distance, class Weight>
Distance add(Distance a, Weight b) {
	static_assert(std::numeric_limits<Distance>::is_integer && std::numeric_limits<Distance>::is_signed);
	static_assert(std::numeric_limits<Weight>::is_integer && std::numeric_limits<Weight>::is_signed);
	static_assert(std::numeric_limits<Distance>::digits >= std::numeric_limits<Weight>::digits);
	Distance value;
#if defined(__GNUC__) || defined(__clang__)
	if (__builtin_add_overflow(a, Distance(b), &value)) throw std::overflow_error("Bellman-Ford distance");
#else
	Distance w = b;
	if ((w > 0 && a > std::numeric_limits<Distance>::max() - w) ||
		(w < 0 && a < std::numeric_limits<Distance>::lowest() - w)) throw std::overflow_error("Bellman-Ford distance");
	value = a + w;
#endif
	if (value == std::numeric_limits<Distance>::max()) throw std::overflow_error("Bellman-Ford infinity");
	return value;
}

template<bool Decompose, int Passes, class Distance, class Graph>
BellmanFordResult<Distance> solve(const Graph &g, int source) {
	assert(std::size_t(g.size()) <= std::size_t(std::numeric_limits<int>::max()));
	int n = static_cast<int>(g.size());
	assert(0 <= source && source < n);
	Distance inf = std::numeric_limits<Distance>::max();
	BellmanFordResult<Distance> result{std::vector<Distance>(n, inf), std::vector<int>(n, -1), {},
		std::vector<std::size_t>(n, std::size_t(-1)), {}, inf, source};
	result.distance[source] = 0;
	result.parent[source] = source;
	auto &d = result.distance;
	std::vector<int> active, order;
	auto relax = [&](int v, std::size_t i) {
		const auto &[to, w] = g[v][i];
		Distance next = add(d[v], w);
		if (next >= d[to]) return false;
		d[to] = next;
		result.parent[to] = v;
		result.parent_edge[to] = i;
		return true;
	};
	for (int pass = 0; pass < Passes; ++pass) {
		bool changed = false;
		for (int v = 0; v < n; ++v) if (d[v] != inf)
			for (std::size_t i = 0; i < g[v].size(); ++i) changed |= relax(v, i);
		if (!changed) return result;
	}
	std::vector<int> position(n, -1);
	std::vector<unsigned char> labeled(n);
	Components components;
	if constexpr (Decompose) components = reachable_components(g, source);
	else {
		components.id.assign(n, 0);
		components.vertices.resize(n);
		for (int v = 0; v < n; ++v) components.vertices[v] = v;
		components.offset.push_back(n);
	}
	struct Frame { int vertex; std::size_t next, incoming; };
	std::vector<Frame> stack;
	using State = std::pair<Distance, int>;
	constexpr bool radix = std::is_integral_v<Distance> && sizeof(Distance) <= 8;
	std::conditional_t<radix, RadixHeap<Distance>, std::vector<State>> heap;
	for (int c = static_cast<int>(components.offset.size()) - 2; c >= 0; --c) {
		auto begin = components.offset[c], end = components.offset[c + 1];
		auto internal = [&](int to) { return components.id[to] == c; };
		bool negative = !Decompose;
		Distance minimum = inf;
		for (auto j = begin; j < end; ++j) {
			int v = components.vertices[j];
			minimum = std::min(minimum, d[v]);
			if constexpr (Decompose)
				for (const auto &[to, w] : g[v]) if (internal(to) && w < 0) negative = true;
		}
		if (minimum == inf) continue;
		if (negative) {
			active.clear();
			for (auto j = begin; j < end; ++j) {
				int v = components.vertices[j];
				if (d[v] != inf) active.push_back(v);
			}
			while (!active.empty()) {
				order.clear();
				for (int v : active) labeled[v] = 0;
				for (int root : active) {
					if (position[root] != -1) continue;
					bool useful = false;
					for (const auto &[to, w] : g[root])
						if (internal(to) && add(d[root], w) < d[to]) { useful = true; break; }
					if (!useful) continue;
					position[root] = 0;
					stack.push_back({root, 0, std::size_t(-1)});
					while (!stack.empty()) {
						int v = stack.back().vertex;
						if (stack.back().next == g[v].size()) {
							position[v] = -2;
							order.push_back(v);
							stack.pop_back();
							continue;
						}
						std::size_t i = stack.back().next++;
						int to = g[v][i].first;
						if (!internal(to) || !relax(v, i)) continue;
						if (position[to] >= 0) {
							for (auto j = std::size_t(position[to]); j < stack.size(); ++j) {
								result.negative_cycle.push_back(stack[j].vertex);
								if (j + 1 < stack.size()) result.negative_cycle_edges.push_back(stack[j + 1].incoming);
							}
							result.negative_cycle.push_back(to);
							result.negative_cycle_edges.push_back(i);
							return result;
						}
						if (position[to] == -1) {
							position[to] = static_cast<int>(stack.size());
							stack.push_back({to, 0, i});
						}
					}
				}
				active.clear();
				for (auto it = order.rbegin(); it != order.rend(); ++it) {
					int v = *it;
					for (std::size_t i = 0; i < g[v].size(); ++i) {
						int to = g[v][i].first;
						if (internal(to) && relax(v, i) && !labeled[to]) active.push_back(to), labeled[to] = 1;
					}
				}
				for (int v : order) position[v] = -1;
			}
		} else if (end - begin > 1) {
			auto push = [&](Distance cost, int v) {
				if constexpr (radix) heap.push({cost, v});
				else { heap.emplace_back(cost, v); std::push_heap(heap.begin(), heap.end(), std::greater<State>()); }
			};
			if constexpr (radix) heap.clear(minimum);
			for (auto j = begin; j < end; ++j) {
				int v = components.vertices[j];
				if (d[v] != inf) push(d[v], v);
			}
			while (!heap.empty()) {
				auto [cost, v] = [&] {
					if constexpr (radix) return heap.pop_value();
					else { std::pop_heap(heap.begin(), heap.end(), std::greater<State>()); auto state = heap.back(); heap.pop_back(); return state; }
				}();
				if (cost != d[v]) continue;
				for (std::size_t i = 0; i < g[v].size(); ++i) {
					int to = g[v][i].first;
					if (internal(to) && relax(v, i)) push(d[to], to);
				}
			}
		}
		if constexpr (Decompose)
			for (auto j = begin; j < end; ++j) {
				int v = components.vertices[j];
				for (std::size_t i = 0; i < g[v].size(); ++i) if (!internal(g[v][i].first)) relax(v, i);
			}
	}
	return result;
}
}

template<class Distance = std::int64_t, class Graph>
BellmanFordResult<Distance> bellman_ford(const Graph &g, int source) {
	return negative_path_detail::solve<true, 2, Distance>(g, source);
}

template<class Distance = std::int64_t, class Graph>
BellmanFordResult<Distance> goldberg_radzik(const Graph &g, int source) {
	return negative_path_detail::solve<false, 0, Distance>(g, source);
}

#endif

Back to top page