johnson.hpp¶
在含负边但无负环的稀疏有向图上求全源最短路，并保留 Dijkstra 在稀疏图上的优势。
Complexity: Time: \(O(VE + V^{2} + VE \log (E + 2))\). Space: \(O(V^{2} + V + E)\).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(VE + V^2 + VE log(E + 2)). Space: O(V^2 + V + E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Johnson all-pairs shortest paths; nullopt denotes a negative cycle
/// anywhere in the graph. Vertices are in [0,n), and E <= INT_MAX.
/// Integer weights and intermediate arithmetic must fit __int128; every
/// finite shortest distance must fit int64_t. int64_t::max() also encodes
/// unreachable pairs.
/// Intermediate reweighted distances may exceed int64_t. Input is unchanged.
template <class Weight>
std::optional<std::vector<std::vector<std::int64_t>>>
johnson(int n, const std::vector<std::tuple<int, int, Weight>> &es) {
assert(n >= 0);
assert(es.size() <= std::size_t(std::numeric_limits<int>::max()));
for (auto [u, v, w] : es) {
(void)w;
assert(0 <= u && u < n && 0 <= v && v < n);
}
std::vector<__int128> pot(n);
for (int round = 0; round < n; ++round) {
bool changed = false;
for (auto [u, v, w] : es) {
__int128 candidate = pot[u] + __int128(w);
if (candidate < pot[v]) {
pot[v] = candidate;
changed = true;
if (round + 1 == n) return std::nullopt;
}
}
if (!changed) break;
}
std::vector<int> head(std::size_t(n) + 1);
__int128 max_weight = 0;
for (auto [u, v, w] : es) {
++head[u + 1];
__int128 reduced = __int128(w) + pot[u] - pot[v];
assert(reduced >= 0);
max_weight = std::max(max_weight, reduced);
}
for (int u = 0; u < n; ++u) head[u + 1] += head[u];
using matrix = std::vector<std::vector<std::int64_t>>;
constexpr auto inf = std::numeric_limits<std::int64_t>::max();
auto solve = [&]<class Distance>() -> matrix {
std::vector<int> to(es.size()), at = head;
std::vector<Distance> weight(es.size());
for (auto [u, v, w] : es) {
int i = at[u]++;
to[i] = v;
weight[i] = Distance(__int128(w) + pot[u] - pot[v]);
}
const Distance infinity = std::numeric_limits<Distance>::max();
matrix result(n, std::vector<std::int64_t>(n, inf));
std::vector<Distance> wide;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) wide.resize(n);
using state = std::pair<Distance, int>;
std::vector<state> heap;
auto compare = std::greater<state>{};
for (int source = 0; source < n; ++source) {
Distance *dist;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) {
std::fill(wide.begin(), wide.end(), infinity);
dist = wide.data();
} else {
// The result row doubles as the narrow Dijkstra distance buffer.
dist = result[source].data();
}
dist[source] = 0;
heap.emplace_back(0, source);
while (!heap.empty()) {
std::pop_heap(heap.begin(), heap.end(), compare);
auto [cur, u] = heap.back();
heap.pop_back();
if (cur != dist[u]) continue;
for (int i = head[u]; i < head[u + 1]; ++i) {
int v = to[i];
Distance w = weight[i];
if (cur <= infinity - w && cur + w < dist[v]) {
dist[v] = cur + w;
heap.emplace_back(dist[v], v);
std::push_heap(heap.begin(), heap.end(), compare);
}
}
}
for (int target = 0; target < n; ++target) if (dist[target] != infinity) {
__int128 value = __int128(dist[target]) - pot[source] + pot[target];
assert(value >= std::numeric_limits<std::int64_t>::min() && value <= inf);
result[source][target] = std::int64_t(value);
}
}
return result;
};
// A shortest reduced path is simple (at most n-1 edges). The max(1,...)
// also bounds individual stored edges when n <= 1; equality needs wide
// distances because the narrow maximum is reserved for unreachable pairs.
// Division avoids overflowing the bound for very large unused edges.
if (max_weight <= (inf - 1) / std::max(1, n - 1))
return solve.template operator()<std::int64_t>();
return solve.template operator()<__int128>();
}
} // namespace noya
// NOYA_JOHNSON_HPP
#ifndef NOYA_JOHNSON_HPP
#define NOYA_JOHNSON_HPP 1
/// @complexity Time: O(VE + V^2 + VE log(E + 2)). Space: O(V^2 + V + E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Johnson all-pairs shortest paths; nullopt denotes a negative cycle
/// anywhere in the graph. Vertices are in [0,n), and E <= INT_MAX.
/// Integer weights and intermediate arithmetic must fit __int128; every
/// finite shortest distance must fit int64_t. int64_t::max() also encodes
/// unreachable pairs.
/// Intermediate reweighted distances may exceed int64_t. Input is unchanged.
template <class Weight>
std::optional<std::vector<std::vector<std::int64_t>>>
johnson(int n, const std::vector<std::tuple<int, int, Weight>> &es) {
assert(n >= 0);
assert(es.size() <= std::size_t(std::numeric_limits<int>::max()));
for (auto [u, v, w] : es) {
(void)w;
assert(0 <= u && u < n && 0 <= v && v < n);
}
std::vector<__int128> pot(n);
for (int round = 0; round < n; ++round) {
bool changed = false;
for (auto [u, v, w] : es) {
__int128 candidate = pot[u] + __int128(w);
if (candidate < pot[v]) {
pot[v] = candidate;
changed = true;
if (round + 1 == n) return std::nullopt;
}
}
if (!changed) break;
}
std::vector<int> head(std::size_t(n) + 1);
__int128 max_weight = 0;
for (auto [u, v, w] : es) {
++head[u + 1];
__int128 reduced = __int128(w) + pot[u] - pot[v];
assert(reduced >= 0);
max_weight = std::max(max_weight, reduced);
}
for (int u = 0; u < n; ++u) head[u + 1] += head[u];
using matrix = std::vector<std::vector<std::int64_t>>;
constexpr auto inf = std::numeric_limits<std::int64_t>::max();
auto solve = [&]<class Distance>() -> matrix {
std::vector<int> to(es.size()), at = head;
std::vector<Distance> weight(es.size());
for (auto [u, v, w] : es) {
int i = at[u]++;
to[i] = v;
weight[i] = Distance(__int128(w) + pot[u] - pot[v]);
}
const Distance infinity = std::numeric_limits<Distance>::max();
matrix result(n, std::vector<std::int64_t>(n, inf));
std::vector<Distance> wide;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) wide.resize(n);
using state = std::pair<Distance, int>;
std::vector<state> heap;
auto compare = std::greater<state>{};
for (int source = 0; source < n; ++source) {
Distance *dist;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) {
std::fill(wide.begin(), wide.end(), infinity);
dist = wide.data();
} else {
// The result row doubles as the narrow Dijkstra distance buffer.
dist = result[source].data();
}
dist[source] = 0;
heap.emplace_back(0, source);
while (!heap.empty()) {
std::pop_heap(heap.begin(), heap.end(), compare);
auto [cur, u] = heap.back();
heap.pop_back();
if (cur != dist[u]) continue;
for (int i = head[u]; i < head[u + 1]; ++i) {
int v = to[i];
Distance w = weight[i];
if (cur <= infinity - w && cur + w < dist[v]) {
dist[v] = cur + w;
heap.emplace_back(dist[v], v);
std::push_heap(heap.begin(), heap.end(), compare);
}
}
}
for (int target = 0; target < n; ++target) if (dist[target] != infinity) {
__int128 value = __int128(dist[target]) - pot[source] + pot[target];
assert(value >= std::numeric_limits<std::int64_t>::min() && value <= inf);
result[source][target] = std::int64_t(value);
}
}
return result;
};
// A shortest reduced path is simple (at most n-1 edges). The max(1,...)
// also bounds individual stored edges when n <= 1; equality needs wide
// distances because the narrow maximum is reserved for unreachable pairs.
// Division avoids overflowing the bound for very large unused edges.
if (max_weight <= (inf - 1) / std::max(1, n - 1))
return solve.template operator()<std::int64_t>();
return solve.template operator()<__int128>();
}
} // namespace noya
#endif // NOYA_JOHNSON_HPP
/// @complexity Time: O(VE + V^2 + VE log(E + 2)). Space: O(V^2 + V + E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Johnson all-pairs shortest paths; nullopt denotes a negative cycle
/// anywhere in the graph. Vertices are in [0,n), and E <= INT_MAX.
/// Integer weights and intermediate arithmetic must fit __int128; every
/// finite shortest distance must fit int64_t. int64_t::max() also encodes
/// unreachable pairs.
/// Intermediate reweighted distances may exceed int64_t. Input is unchanged.
template <class Weight>
std::optional<std::vector<std::vector<std::int64_t>>>
johnson(int n, const std::vector<std::tuple<int, int, Weight>> &es) {
assert(n >= 0);
assert(es.size() <= std::size_t(std::numeric_limits<int>::max()));
for (auto [u, v, w] : es) {
(void)w;
assert(0 <= u && u < n && 0 <= v && v < n);
}
std::vector<__int128> pot(n);
for (int round = 0; round < n; ++round) {
bool changed = false;
for (auto [u, v, w] : es) {
__int128 candidate = pot[u] + __int128(w);
if (candidate < pot[v]) {
pot[v] = candidate;
changed = true;
if (round + 1 == n) return std::nullopt;
}
}
if (!changed) break;
}
std::vector<int> head(std::size_t(n) + 1);
__int128 max_weight = 0;
for (auto [u, v, w] : es) {
++head[u + 1];
__int128 reduced = __int128(w) + pot[u] - pot[v];
assert(reduced >= 0);
max_weight = std::max(max_weight, reduced);
}
for (int u = 0; u < n; ++u) head[u + 1] += head[u];
using matrix = std::vector<std::vector<std::int64_t>>;
constexpr auto inf = std::numeric_limits<std::int64_t>::max();
auto solve = [&]<class Distance>() -> matrix {
std::vector<int> to(es.size()), at = head;
std::vector<Distance> weight(es.size());
for (auto [u, v, w] : es) {
int i = at[u]++;
to[i] = v;
weight[i] = Distance(__int128(w) + pot[u] - pot[v]);
}
const Distance infinity = std::numeric_limits<Distance>::max();
matrix result(n, std::vector<std::int64_t>(n, inf));
std::vector<Distance> wide;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) wide.resize(n);
using state = std::pair<Distance, int>;
std::vector<state> heap;
auto compare = std::greater<state>{};
for (int source = 0; source < n; ++source) {
Distance *dist;
if constexpr (sizeof(Distance) > sizeof(std::int64_t)) {
std::fill(wide.begin(), wide.end(), infinity);
dist = wide.data();
} else {
// The result row doubles as the narrow Dijkstra distance buffer.
dist = result[source].data();
}
dist[source] = 0;
heap.emplace_back(0, source);
while (!heap.empty()) {
std::pop_heap(heap.begin(), heap.end(), compare);
auto [cur, u] = heap.back();
heap.pop_back();
if (cur != dist[u]) continue;
for (int i = head[u]; i < head[u + 1]; ++i) {
int v = to[i];
Distance w = weight[i];
if (cur <= infinity - w && cur + w < dist[v]) {
dist[v] = cur + w;
heap.emplace_back(dist[v], v);
std::push_heap(heap.begin(), heap.end(), compare);
}
}
}
for (int target = 0; target < n; ++target) if (dist[target] != infinity) {
__int128 value = __int128(dist[target]) - pot[source] + pot[target];
assert(value >= std::numeric_limits<std::int64_t>::min() && value <= inf);
result[source][target] = std::int64_t(value);
}
}
return result;
};
// A shortest reduced path is simple (at most n-1 edges). The max(1,...)
// also bounds individual stored edges when n <= 1; equality needs wide
// distances because the narrow maximum is reserved for unreachable pairs.
// Division avoids overflowing the bound for very large unused edges.
if (max_weight <= (inf - 1) / std::max(1, n - 1))
return solve.template operator()<std::int64_t>();
return solve.template operator()<__int128>();
}
} // namespace noya
// NOYA_JOHNSON_HPP