ワーシャルフロイド法

graph/warshall_floyd.hpp · GitHub で見る

共通ページ

重み付き有向グラフの全点対最短距離を求めます。負辺・多重辺・自己辺に対応し、グラフ全体の負閉路を検出します。頂点番号は 0 始まりです。

#include "lib/graph/warshall_floyd.hpp"

int main() {
	WarshallFloyd<> graph(4);
	graph.add_edge(0, 1, 5);
	graph.add_edge(1, 2, -7);
	graph.add_edge(0, 2, 3);
	if (!graph.solve()) return 0;
	assert(graph(0, 2) == -2);
	assert(graph(3, 0) == graph.infinity());
}

距離型 T の既定は std::int64_t です。符号付きの 32 bit・64 bit 整数型を指定できます。到達不能を表す inf の既定は std::numeric_limits<T>::max()。有限の辺重み・最短距離は T の範囲内かつ inf 未満である必要があります。有限の負距離には -1 も使えるため、到達不能の判定は infinity() との比較で行います。

solve() は負閉路がなければ true、あれば false を返します。false の後の距離は計算途中の値で、最短距離として使えません。負無限大になる全ての頂点対を列挙する機能や経路復元はありません。

行列から構築したり、既存の行列を更新したりすることもできます。対角成分は 0 と入力値の小さい方、到達不能な成分は inf としてください。

#include "lib/graph/warshall_floyd.hpp"

std::int64_t inf = std::numeric_limits<std::int64_t>::max();
std::vector<std::vector<std::int64_t>> d{
	{0, 5, inf},
	{inf, 0, -7},
	{inf, inf, 0}
};
bool ok = warshall_floyd(d);
assert(ok && d[0][2] == -2);

add_edge は既存の辺との最小値を保持します。成功した solve() の後に辺を追加して再計算することもできます。辺の削除や重みの増加には、新しいオブジェクトを構築してください。

負閉路がなく有限の最短距離が T の範囲外または inf 以上なら std::overflow_error を送出します。

n < 0、inf <= 0、正方行列でない入力、weight >= inf は std::invalid_argument です。頂点番号は範囲内である必要があります。

WARSHALL_FLOYD_DISABLE_AVX512、WARSHALL_FLOYD_DISABLE_SIMD、WARSHALL_FLOYD_DISABLE_NARROWING を include 前に定義すると、それぞれ AVX-512・明示的な SIMD・64 bit から 32 bit への変換を無効にできます。

実行時間は CPU・重みの範囲・負辺の有無・入出力に依存します。

頂点数を \(n\)、32 の倍数に切り上げた頂点数を \(p=32\lceil n/32\rceil\) とします。

操作	内容・返り値	計算量
WarshallFloyd<T>(n,inf)	辺なし、対角 0 の行列を構築。inf は省略可	時間・空間 \(O(p^2)\)
WarshallFloyd<T>(matrix,inf)	正方行列から構築	時間・空間 \(O(p^2)\)
size() / infinity()	頂点数 / 到達不能値	\(O(1)\)
add_edge(from,to,weight)	有向辺を追加	\(O(1)\)
solve()	全点対最短距離を計算し、負閉路がなければ true	時間 \(O(p^3)\)、型変換時の追加空間 \(O(p^2)\)
operator()(i,j)	現在の距離を値で取得	\(O(1)\)
distances()	現在の行列を vector<vector<T>> で取得	時間・追加空間 \(O(n^2)\)
warshall_floyd(matrix,inf)	行列を更新し、負閉路がなければ true	時間 \(O(p^3)\)、追加空間 \(O(p^2)\)

本体は sizeof(T) * p * p バイトを保持します。型変換時は内部型の同じ要素数の作業領域が加わります。例えば 64 bit から 32 bit へ変換する場合、本体と作業領域で合計約 12 * p * p バイトです。行列版は呼び出し元の行列も別に保持するため、大きい入力では辺から直接 WarshallFloyd を構築するとメモリを抑えられます。

依存ライブラリ

• 整列した配列のアロケータ
• ワーシャルフロイド法の内部演算

検証コード

• test/warshall_floyd.test.cpp

コード

#ifndef WARSHALL_FLOYD_HPP
#define WARSHALL_FLOYD_HPP

#include "../components/warshall_floyd_kernel.hpp"
#include <stdexcept>

template<class T = std::int64_t>
class WarshallFloyd {
	static_assert(std::is_integral_v<T> && std::is_signed_v<T> && (sizeof(T) == 4 || sizeof(T) == 8));
	static constexpr int block = 32;
	int n_ = 0, blocks_ = 0;
	T infinity_ = std::numeric_limits<T>::max();
	warshall_floyd_detail::Storage<T> data_;
	static int checked_size(std::size_t n) {
		if (n > std::size_t(std::numeric_limits<int>::max())) throw std::length_error("Warshall-Floyd size");
		return int(n);
	}
	std::size_t index(int i, int j) const {
		return (std::size_t(i / block) * blocks_ + j / block) * block * block + i % block * block + j % block;
	}
	template<class U> bool run_as(bool negative) {
		const U inf = std::numeric_limits<U>::max();
		warshall_floyd_detail::Storage<U> work(data_.size());
		for (std::size_t i = 0; i < data_.size(); ++i) work[i] = data_[i] == infinity_ ? inf : U(data_[i]);
		bool success = negative ? warshall_floyd_detail::run<U, block, true>(work.data(), blocks_, inf)
			: warshall_floyd_detail::run<U, block, false>(work.data(), blocks_, inf);
		if (success)
			for (U value : work) if (value != inf && (__int128_t(value) < std::numeric_limits<T>::lowest() || __int128_t(value) >= infinity_))
				throw std::overflow_error("Warshall-Floyd distance");
		for (std::size_t i = 0; i < data_.size(); ++i) {
			if (work[i] == inf) data_[i] = infinity_;
			else data_[i] = T(std::clamp(__int128_t(work[i]), __int128_t(std::numeric_limits<T>::lowest()), __int128_t(infinity_) - 1));
		}
		return success;
	}

  public:
	using value_type = T;
	explicit WarshallFloyd(int n = 0, T inf = std::numeric_limits<T>::max()) : n_(n), infinity_(inf) {
		if (n < 0 || inf <= 0) throw std::invalid_argument("Warshall-Floyd shape or infinity");
		blocks_ = n / block + (n % block != 0);
		std::size_t padded = std::size_t(blocks_) * block;
		data_.assign(padded * padded, infinity_);
		for (std::size_t i = 0; i < padded; ++i) data_[index(int(i), int(i))] = 0;
	}
	explicit WarshallFloyd(const std::vector<std::vector<T>>& distance, T inf = std::numeric_limits<T>::max())
		: WarshallFloyd(checked_size(distance.size()), inf) {
		for (int i = 0; i < n_; ++i) {
			if (distance[i].size() != std::size_t(n_)) throw std::invalid_argument("Warshall-Floyd matrix");
			for (int j = 0; j < n_; ++j) {
				if (distance[i][j] == infinity_) continue;
				add_edge(i, j, distance[i][j]);
			}
		}
	}
	int size() const { return n_; }
	T infinity() const { return infinity_; }
	T operator()(int i, int j) const {
		assert(0 <= i && i < n_ && 0 <= j && j < n_);
		return data_[index(i, j)];
	}
	void add_edge(int from, int to, T weight) {
		assert(0 <= from && from < n_ && 0 <= to && to < n_);
		if (weight >= infinity_) throw std::invalid_argument("Warshall-Floyd edge");
		auto& value = data_[index(from, to)];
		value = std::min(value, weight);
	}
	bool solve() {
		__int128_t maximum = 0;
		std::size_t finite = 0;
		bool negative = false;
		for (T value : data_) if (value != infinity_) {
			++finite;
			__int128_t x = value;
			negative |= x < 0;
			maximum = std::max(maximum, x < 0 ? -x : x);
		}
		__int128_t bound = 2 * (__int128_t(n_) + 1) * maximum;
		bool complete = finite == std::size_t(n_) * n_ + std::size_t(blocks_) * block - n_;
		if constexpr (sizeof(T) == 8) {
#ifndef WARSHALL_FLOYD_DISABLE_NARROWING
			__int128_t limit = std::min(__int128_t(infinity_), __int128_t(std::numeric_limits<std::int32_t>::max()));
			if (negative ? bound < limit : maximum < limit && (complete || (__int128_t(n_) - 1) * maximum < limit))
				return run_as<std::int32_t>(negative);
#endif
		}
		if (negative ? bound >= infinity_ : !complete && (__int128_t(n_) - 1) * maximum >= infinity_) {
			if constexpr (sizeof(T) == 4) {
				if (bound < std::numeric_limits<std::int64_t>::max()) return run_as<std::int64_t>(negative);
			}
			return run_as<__int128_t>(negative);
		}
		return negative ? warshall_floyd_detail::run<T, block, true>(data_.data(), blocks_, infinity_)
			: warshall_floyd_detail::run<T, block, false>(data_.data(), blocks_, infinity_);
	}
	std::vector<std::vector<T>> distances() const {
		std::vector<std::vector<T>> result(n_, std::vector<T>(n_));
		for (int i = 0; i < n_; ++i) for (int j = 0; j < n_; ++j) result[i][j] = data_[index(i, j)];
		return result;
	}
};

template<class T>
bool warshall_floyd(std::vector<std::vector<T>>& distance, T inf = std::numeric_limits<T>::max()) {
	WarshallFloyd<T> work(distance, inf);
	bool success = work.solve();
	distance = work.distances();
	return success;
}

#endif

Back to top page