johnson.hpp

SECTIONGraph INCLUDEnoya/johnson.hpp

Johnson all-pairs shortest paths for a sparse directed graph in O(nm + nm log n); nullopt means the graph contains a negative cycle.

在含负边但无负环的稀疏有向图上求全源最短路，并保留 Dijkstra 在稀疏图上的优势。

Implementation

View on GitHub

SourceExpanded

#ifndef NOYA_JOHNSON_HPP
#define NOYA_JOHNSON_HPP 1

/// @complexity Time: O(VE + V(E + V) log V).
/// Space: O(V^2 + E).

#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>

namespace noya {

/// @brief Johnson all-pairs shortest paths for a sparse directed graph in
/// O(nm + nm log n); nullopt means the graph contains a negative cycle.
template <class Weight>
std::optional<std::vector<std::vector<std::int64_t>>>
johnson(int n, const std::vector<std::tuple<int, int, Weight>> &edges) {
  assert(n >= 0);
  for (auto [from, to, weight] : edges) {
    (void)weight;
    assert(0 <= from && from < n);
    assert(0 <= to && to < n);
  }
  std::vector<__int128> potential(n);
  for (int iteration = 0; iteration < n; iteration++) {
    bool changed = false;
    for (auto [from, to, weight] : edges) {
      __int128 candidate = potential[from] + __int128(weight);
      if (candidate < potential[to]) {
        potential[to] = candidate;
        changed = true;
        if (iteration + 1 == n) {
          return std::nullopt;
        }
      }
    }
    if (!changed) {
      break;
    }
  }

  std::vector<std::vector<std::pair<int, std::int64_t>>> graph(n);
  for (auto [from, to, weight] : edges) {
    __int128 reweighted = __int128(weight) + potential[from] - potential[to];
    assert(reweighted >= 0 &&
           reweighted <= std::numeric_limits<std::int64_t>::max());
    graph[from].emplace_back(to, std::int64_t(reweighted));
  }

  const std::int64_t infinity = std::numeric_limits<std::int64_t>::max();
  std::vector<std::vector<std::int64_t>> result(
      n, std::vector<std::int64_t>(n, infinity));
  for (int source = 0; source < n; source++) {
    using state = std::pair<std::int64_t, int>;
    std::priority_queue<state, std::vector<state>, std::greater<>> queue;
    std::vector<std::int64_t> distance(n, infinity);
    distance[source] = 0;
    queue.emplace(0, source);
    while (!queue.empty()) {
      auto [current, vertex] = queue.top();
      queue.pop();
      if (current != distance[vertex]) {
        continue;
      }
      for (auto [next, weight] : graph[vertex]) {
        if (current <= infinity - weight && current + weight < distance[next]) {
          distance[next] = current + weight;
          queue.emplace(distance[next], next);
        }
      }
    }
    for (int target = 0; target < n; target++) {
      if (distance[target] == infinity) {
        continue;
      }
      __int128 value =
          __int128(distance[target]) - potential[source] + potential[target];
      assert(std::numeric_limits<std::int64_t>::lowest() <= value &&
             value <= std::numeric_limits<std::int64_t>::max());
      result[source][target] = std::int64_t(value);
    }
  }
  return result;
}

} // namespace noya

#endif // NOYA_JOHNSON_HPP

#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <optional>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>

/// @complexity Time: O(VE + V(E + V) log V).
/// Space: O(V^2 + E).

namespace noya {

/// @brief Johnson all-pairs shortest paths for a sparse directed graph in
/// O(nm + nm log n); nullopt means the graph contains a negative cycle.
template <class Weight>
std::optional<std::vector<std::vector<std::int64_t>>>
johnson(int n, const std::vector<std::tuple<int, int, Weight>> &edges) {
  assert(n >= 0);
  for (auto [from, to, weight] : edges) {
    (void)weight;
    assert(0 <= from && from < n);
    assert(0 <= to && to < n);
  }
  std::vector<__int128> potential(n);
  for (int iteration = 0; iteration < n; iteration++) {
    bool changed = false;
    for (auto [from, to, weight] : edges) {
      __int128 candidate = potential[from] + __int128(weight);
      if (candidate < potential[to]) {
        potential[to] = candidate;
        changed = true;
        if (iteration + 1 == n) {
          return std::nullopt;
        }
      }
    }
    if (!changed) {
      break;
    }
  }

  std::vector<std::vector<std::pair<int, std::int64_t>>> graph(n);
  for (auto [from, to, weight] : edges) {
    __int128 reweighted = __int128(weight) + potential[from] - potential[to];
    assert(reweighted >= 0 &&
           reweighted <= std::numeric_limits<std::int64_t>::max());
    graph[from].emplace_back(to, std::int64_t(reweighted));
  }

  const std::int64_t infinity = std::numeric_limits<std::int64_t>::max();
  std::vector<std::vector<std::int64_t>> result(
      n, std::vector<std::int64_t>(n, infinity));
  for (int source = 0; source < n; source++) {
    using state = std::pair<std::int64_t, int>;
    std::priority_queue<state, std::vector<state>, std::greater<>> queue;
    std::vector<std::int64_t> distance(n, infinity);
    distance[source] = 0;
    queue.emplace(0, source);
    while (!queue.empty()) {
      auto [current, vertex] = queue.top();
      queue.pop();
      if (current != distance[vertex]) {
        continue;
      }
      for (auto [next, weight] : graph[vertex]) {
        if (current <= infinity - weight && current + weight < distance[next]) {
          distance[next] = current + weight;
          queue.emplace(distance[next], next);
        }
      }
    }
    for (int target = 0; target < n; target++) {
      if (distance[target] == infinity) {
        continue;
      }
      __int128 value =
          __int128(distance[target]) - potential[source] + potential[target];
      assert(std::numeric_limits<std::int64_t>::lowest() <= value &&
             value <= std::numeric_limits<std::int64_t>::max());
      result[source][target] = std::int64_t(value);
    }
  }
  return result;
}

} // namespace noya