m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Shortest Path
(graph/shortest_path.hpp)

Overview

graph/shortest_path.hpp includes the shortest-path algorithms whose behavior respects the adjacency stored in Graph<T>.

Most of these algorithms are direction-respecting and can be used on directed graphs as written or on undirected graphs built with add_edge. The exception is dag_shortest_path, which is specifically for directed acyclic graphs.

Included Headers

Header Graph orientation Contents
graph/bfs.hpp Direction-respecting Shortest paths by number of edges.
graph/zero_one_bfs.hpp Direction-respecting Shortest paths with edge costs 0 or 1.
graph/dag_shortest_path.hpp Directed DAG only Shortest paths in a DAG, including negative edge costs.
graph/dijkstra.hpp Direction-respecting Non-negative weighted shortest paths.
graph/k_shortest_walk.hpp Direction-respecting The first k walk lengths with non-negative edge costs.
graph/bellman_ford.hpp Direction-respecting Shortest paths with negative edges.
graph/warshall_floyd.hpp Direction-respecting All-pairs shortest paths.

Complexity

This header is an include bundle and provides no runtime operation by itself. See the included algorithm pages for public interfaces and complexities.

Depends on

Required by

Verified with

Code

#ifndef M1UNE_GRAPH_SHORTEST_PATH_HPP
#define M1UNE_GRAPH_SHORTEST_PATH_HPP 1

#include "bellman_ford.hpp"
#include "bfs.hpp"
#include "cow_game.hpp"
#include "dag_shortest_path.hpp"
#include "dijkstra.hpp"
#include "k_shortest_walk.hpp"
#include "warshall_floyd.hpp"
#include "zero_one_bfs.hpp"

#endif  // M1UNE_GRAPH_SHORTEST_PATH_HPP
#line 1 "graph/shortest_path.hpp"



#line 1 "graph/bellman_ford.hpp"



#include <algorithm>
#include <cassert>
#include <limits>
#include <queue>
#include <vector>

#line 1 "graph/graph.hpp"



#include <array>
#line 6 "graph/graph.hpp"
#include <utility>
#line 8 "graph/graph.hpp"

namespace m1une {
namespace graph {

template <class T = int>
struct Edge {
    using cost_type = T;

    int from;
    int to;
    T cost;
    int id;
    bool alive;

    Edge() : from(-1), to(-1), cost(T()), id(-1), alive(true) {}
    Edge(int from_, int to_, T cost_ = T(1), int id_ = -1, bool alive_ = true)
        : from(from_), to(to_), cost(cost_), id(id_), alive(alive_) {}

    int other(int v) const {
        assert(v == from || v == to);
        return from ^ to ^ v;
    }
};

template <class T = int>
struct Graph {
    using edge_type = Edge<T>;
    using cost_type = T;

   private:
    struct EdgePositions {
        std::array<std::pair<int, int>, 2> value{};
        int size = 0;

        void push_back(std::pair<int, int> position) {
            assert(size < 2);
            value[size++] = position;
        }
    };

    int _n;
    int _edge_count;
    std::vector<std::vector<edge_type>> _g;
    std::vector<EdgePositions> _edge_positions;

   public:
    Graph() : _n(0), _edge_count(0) {}
    explicit Graph(int n) : _n(n), _edge_count(0), _g(n) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    bool empty() const {
        return _n == 0;
    }

    int edge_count() const {
        return _edge_count;
    }

    int add_vertex() {
        _g.emplace_back();
        return _n++;
    }

    int add_directed_edge(int from, int to, T cost = T(1)) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        int id = _edge_count++;
        int idx = int(_g[from].size());
        _g[from].push_back(edge_type(from, to, cost, id));
        _edge_positions.emplace_back();
        _edge_positions.back().push_back({from, idx});
        return id;
    }

    int add_edge(int u, int v, T cost = T(1)) {
        assert(0 <= u && u < _n);
        assert(0 <= v && v < _n);
        int id = _edge_count++;
        int u_idx = int(_g[u].size());
        _g[u].push_back(edge_type(u, v, cost, id));
        int v_idx = int(_g[v].size());
        _g[v].push_back(edge_type(v, u, cost, id));
        _edge_positions.emplace_back();
        _edge_positions.back().push_back({u, u_idx});
        _edge_positions.back().push_back({v, v_idx});
        return id;
    }

    void set_edge_alive(int id, bool alive) {
        assert(0 <= id && id < _edge_count);
        for (int i = 0; i < _edge_positions[id].size; ++i) {
            auto [v, idx] = _edge_positions[id].value[i];
            _g[v][idx].alive = alive;
        }
    }

    void erase_edge(int id) {
        set_edge_alive(id, false);
    }

    void revive_edge(int id) {
        set_edge_alive(id, true);
    }

    bool is_edge_alive(int id) const {
        assert(0 <= id && id < _edge_count);
        assert(_edge_positions[id].size != 0);
        auto [v, idx] = _edge_positions[id].value[0];
        return _g[v][idx].alive;
    }

    const std::vector<edge_type>& operator[](int v) const {
        assert(0 <= v && v < _n);
        return _g[v];
    }

    std::vector<edge_type>& operator[](int v) {
        assert(0 <= v && v < _n);
        return _g[v];
    }

    const std::vector<std::vector<edge_type>>& adjacency() const {
        return _g;
    }

    std::vector<std::vector<edge_type>>& adjacency() {
        return _g;
    }

    std::vector<edge_type> edges(bool include_inactive = false) const {
        std::vector<edge_type> result;
        result.reserve(_edge_count);
        std::vector<char> used(_edge_count, false);
        for (int v = 0; v < _n; v++) {
            for (const auto& e : _g[v]) {
                if (!include_inactive && !e.alive) continue;
                if (0 <= e.id && e.id < _edge_count) {
                    if (used[e.id]) continue;
                    used[e.id] = true;
                }
                result.push_back(e);
            }
        }
        return result;
    }

    Graph reversed() const {
        Graph result(_n);
        result._edge_count = _edge_count;
        result._edge_positions.assign(_edge_count, {});
        for (int v = 0; v < _n; v++) {
            for (const auto& e : _g[v]) {
                int idx = int(result._g[e.to].size());
                result._g[e.to].push_back(edge_type(e.to, e.from, e.cost, e.id, e.alive));
                if (0 <= e.id && e.id < _edge_count) result._edge_positions[e.id].push_back({e.to, idx});
            }
        }
        return result;
    }
};

}  // namespace graph
}  // namespace m1une


#line 11 "graph/bellman_ford.hpp"

namespace m1une {
namespace graph {

template <class T>
struct BellmanFordResult {
    std::vector<T> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<bool> negative;
    T inf;
    bool has_negative_cycle;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    bool affected_by_negative_cycle(int v) const {
        assert(0 <= v && v < int(negative.size()));
        return negative[v];
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        assert(!affected_by_negative_cycle(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
BellmanFordResult<T> bellman_ford(const Graph<T>& g, const std::vector<int>& sources,
                                  T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    BellmanFordResult<T> result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.negative.assign(n, false);
    result.inf = inf;
    result.has_negative_cycle = false;

    for (int s : sources) {
        assert(0 <= s && s < n);
        result.dist[s] = T(0);
    }

    std::vector<int> relaxed_vertices;
    for (int iter = 0; iter < n; iter++) {
        bool updated = false;
        for (int v = 0; v < n; v++) {
            if (result.dist[v] == inf) continue;
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                T nd = result.dist[v] + e.cost;
                if (result.dist[e.to] <= nd) continue;
                result.dist[e.to] = nd;
                result.parent[e.to] = v;
                result.parent_edge[e.to] = e.id;
                updated = true;
                if (iter == n - 1) relaxed_vertices.push_back(e.to);
            }
        }
        if (!updated) break;
    }

    std::queue<int> que;
    for (int v : relaxed_vertices) {
        if (result.negative[v]) continue;
        result.negative[v] = true;
        que.push(v);
    }
    while (!que.empty()) {
        int v = que.front();
        que.pop();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (result.negative[e.to]) continue;
            result.negative[e.to] = true;
            que.push(e.to);
        }
    }

    for (bool x : result.negative) result.has_negative_cycle = result.has_negative_cycle || x;
    return result;
}

template <class T>
BellmanFordResult<T> bellman_ford(const Graph<T>& g, int s, T inf = std::numeric_limits<T>::max() / T(4)) {
    return bellman_ford(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/bfs.hpp"



#line 6 "graph/bfs.hpp"
#include <concepts>
#include <functional>
#line 11 "graph/bfs.hpp"

#line 13 "graph/bfs.hpp"

namespace m1une {
namespace graph {

struct BfsResult {
    std::vector<int> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != -1;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

namespace bfs_detail {

template <class Callback>
concept BfsCallback =
    std::invocable<Callback&, int, int> ||
    std::invocable<Callback&, int>;

template <BfsCallback Callback>
void invoke_callback(Callback& callback, int vertex, int parent) {
    if constexpr (std::invocable<Callback&, int, int>) {
        std::invoke(callback, vertex, parent);
    } else {
        std::invoke(callback, vertex);
    }
}

template <class T, class Callback>
BfsResult run_bfs(
    const Graph<T>& g,
    const std::vector<int>& sources,
    Callback& callback
) {
    int n = g.size();
    BfsResult result;
    result.dist.assign(n, -1);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);

    std::queue<int> que;
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] != -1) continue;
        result.dist[s] = 0;
        invoke_callback(callback, s, -1);
        que.push(s);
    }

    while (!que.empty()) {
        int v = que.front();
        que.pop();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (result.dist[e.to] != -1) continue;
            result.dist[e.to] = result.dist[v] + 1;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            invoke_callback(callback, e.to, v);
            que.push(e.to);
        }
    }

    return result;
}

}  // namespace bfs_detail

template <class T>
BfsResult bfs(const Graph<T>& g, const std::vector<int>& sources) {
    auto callback = [](int) {};
    return bfs_detail::run_bfs(g, sources, callback);
}

template <class T>
BfsResult bfs(const Graph<T>& g, int s) {
    return bfs(g, std::vector<int>{s});
}

template <class T, class Callback>
requires bfs_detail::BfsCallback<Callback>
BfsResult bfs(
    const Graph<T>& g,
    const std::vector<int>& sources,
    Callback&& callback
) {
    return bfs_detail::run_bfs(g, sources, callback);
}

template <class T, class Callback>
requires bfs_detail::BfsCallback<Callback>
BfsResult bfs(const Graph<T>& g, int source, Callback&& callback) {
    return bfs(
        g,
        std::vector<int>{source},
        std::forward<Callback>(callback)
    );
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/cow_game.hpp"



#line 6 "graph/cow_game.hpp"
#include <optional>
#include <type_traits>
#line 10 "graph/cow_game.hpp"

namespace m1une {
namespace graph {

template <class T>
struct CowGameConstraint {
    int a;
    int b;
    T upper_bound;
};

template <class T>
struct CowGameSolution {
    bool feasible = false;
    std::vector<T> value;

    bool is_feasible() const {
        return feasible;
    }
};

template <class T>
struct CowGameUpperBounds {
    bool feasible;
    std::vector<T> upper_bound;
    T inf;

    bool is_feasible() const {
        return feasible;
    }

    bool bounded(int variable) const {
        assert(0 <= variable && variable < int(upper_bound.size()));
        return feasible && upper_bound[variable] != inf;
    }
};

template <class T>
struct CowGameDifferenceBounds {
    bool feasible;
    std::optional<T> lower_bound;
    std::optional<T> upper_bound;

    bool is_feasible() const {
        return feasible;
    }

    bool bounded_below() const {
        return feasible && lower_bound.has_value();
    }

    bool bounded_above() const {
        return feasible && upper_bound.has_value();
    }
};

template <class T>
class CowGame {
    static_assert(std::is_arithmetic_v<T> && std::is_signed_v<T>);

    struct RelaxationResult {
        bool has_negative_cycle;
        std::vector<T> dist;
    };

    int _n;
    std::vector<CowGameConstraint<T>> _constraints;
    std::vector<std::vector<int>> _outgoing_constraints;
    bool _has_negative_upper_bound = false;
    mutable bool _solution_cached = false;
    mutable CowGameSolution<T> _cached_solution;

    void assert_variable(int variable) const {
        (void)variable;
        assert(0 <= variable && variable < _n);
    }

    T negate(T value) const {
        assert(value != std::numeric_limits<T>::lowest());
        return -value;
    }

    RelaxationResult check_feasibility() const {
        std::vector<T> dist(_n, T());
        for (int iteration = 0; iteration < _n; iteration++) {
            bool updated = false;
            for (const auto& constraint : _constraints) {
                T candidate = dist[constraint.b] + constraint.upper_bound;
                if (dist[constraint.a] <= candidate) continue;
                dist[constraint.a] = candidate;
                updated = true;
                if (iteration == _n - 1) return RelaxationResult{true, std::move(dist)};
            }
            if (!updated) break;
        }
        return RelaxationResult{false, std::move(dist)};
    }

    std::vector<T> shortest_paths(int source, T inf) const {
        const auto& potential = _cached_solution.value;
        std::vector<T> dist(_n, inf);
        std::vector<int> heap;
        // -1 is unseen, -2 is fixed, and every other value is a heap index.
        std::vector<int> position(_n, -1);
        heap.reserve(_n);

        auto swap_heap = [&](int i, int j) {
            std::swap(heap[i], heap[j]);
            position[heap[i]] = i;
            position[heap[j]] = j;
        };
        auto sift_up = [&](int i) {
            while (i > 0) {
                int parent = (i - 1) / 2;
                if (dist[heap[parent]] <= dist[heap[i]]) break;
                swap_heap(parent, i);
                i = parent;
            }
        };
        auto sift_down = [&](int i) {
            while (2 * i + 1 < int(heap.size())) {
                int child = 2 * i + 1;
                if (child + 1 < int(heap.size()) &&
                    dist[heap[child + 1]] < dist[heap[child]]) {
                    child++;
                }
                if (dist[heap[i]] <= dist[heap[child]]) break;
                swap_heap(i, child);
                i = child;
            }
        };

        dist[source] = T();
        position[source] = 0;
        heap.push_back(source);

        while (!heap.empty()) {
            int b = heap[0];
            position[b] = -2;
            int last = heap.back();
            heap.pop_back();
            if (!heap.empty()) {
                heap[0] = last;
                position[last] = 0;
                sift_down(0);
            }

            for (int id : _outgoing_constraints[b]) {
                const auto& constraint = _constraints[id];
                T cost = constraint.upper_bound + potential[b] -
                         potential[constraint.a];
                assert(cost >= T());
                T candidate = dist[b] + cost;
                if (dist[constraint.a] <= candidate) continue;
                dist[constraint.a] = candidate;
                assert(position[constraint.a] != -2);
                if (position[constraint.a] == -1) {
                    position[constraint.a] = int(heap.size());
                    heap.push_back(constraint.a);
                }
                sift_up(position[constraint.a]);
            }
        }

        for (int v = 0; v < _n; v++) {
            if (dist[v] == inf) continue;
            dist[v] = dist[v] - potential[source] + potential[v];
        }
        return dist;
    }

   public:
    CowGame() : CowGame(0) {}

    explicit CowGame(int variable_count)
        : _n(variable_count),
          _outgoing_constraints(variable_count < 0 ? 0 : variable_count) {
        assert(variable_count >= 0);
    }

    int size() const {
        return _n;
    }

    int constraint_count() const {
        return int(_constraints.size());
    }

    const CowGameConstraint<T>& get_constraint(int id) const {
        assert(0 <= id && id < int(_constraints.size()));
        return _constraints[id];
    }

    const std::vector<CowGameConstraint<T>>& constraints() const {
        return _constraints;
    }

    bool can_use_dijkstra() const {
        return !_has_negative_upper_bound ||
               (_solution_cached && _cached_solution.feasible);
    }

    int add_upper_bound(int a, int b, T upper_bound) {
        assert_variable(a);
        assert_variable(b);
        int id = int(_constraints.size());
        _constraints.push_back(CowGameConstraint<T>{a, b, upper_bound});
        _outgoing_constraints[b].push_back(id);
        _has_negative_upper_bound = _has_negative_upper_bound || upper_bound < T();
        _solution_cached = false;
        return id;
    }

    int add_constraint(int a, int b, T upper_bound) {
        return add_upper_bound(a, b, upper_bound);
    }

    int add_lower_bound(int a, int b, T lower_bound) {
        return add_upper_bound(b, a, negate(lower_bound));
    }

    void add_bounds(int a, int b, T lower_bound, T upper_bound) {
        assert(lower_bound <= upper_bound);
        add_lower_bound(a, b, lower_bound);
        add_upper_bound(a, b, upper_bound);
    }

    void add_equality(int a, int b, T difference) {
        add_bounds(a, b, difference, difference);
    }

    CowGameSolution<T> solve() const {
        if (_solution_cached) return _cached_solution;

        _cached_solution.feasible = true;
        _cached_solution.value.assign(_n, T());
        if (_has_negative_upper_bound) {
            auto result = check_feasibility();
            _cached_solution.feasible = !result.has_negative_cycle;
            _cached_solution.value.clear();
            if (_cached_solution.feasible) {
                _cached_solution.value = std::move(result.dist);
            }
        }
        _solution_cached = true;
        return _cached_solution;
    }

    bool is_feasible() const {
        if (!_solution_cached) (void)solve();
        return _cached_solution.feasible;
    }

    CowGameUpperBounds<T> tightest_upper_bounds(int source) const {
        assert_variable(source);
        T inf = std::numeric_limits<T>::max() / T(4);
        CowGameUpperBounds<T> result;
        result.feasible = is_feasible();
        result.inf = inf;
        result.upper_bound.assign(_n, inf);
        if (!result.feasible) return result;

        result.upper_bound = shortest_paths(source, inf);
        return result;
    }

    CowGameDifferenceBounds<T> difference_bounds(int a, int b) const {
        assert_variable(a);
        assert_variable(b);
        T inf = std::numeric_limits<T>::max() / T(4);
        CowGameDifferenceBounds<T> result;
        result.feasible = is_feasible();
        if (!result.feasible) return result;

        auto upper = shortest_paths(b, inf);
        if (upper[a] != inf) result.upper_bound = upper[a];

        auto lower = shortest_paths(a, inf);
        if (lower[b] != inf) result.lower_bound = negate(lower[b]);
        return result;
    }
};

template <class T>
using DifferenceConstraints = CowGame<T>;

}  // namespace graph
}  // namespace m1une


#line 1 "graph/dag_shortest_path.hpp"



#line 9 "graph/dag_shortest_path.hpp"

#line 1 "graph/topological_sort.hpp"



#line 7 "graph/topological_sort.hpp"

#line 9 "graph/topological_sort.hpp"

namespace m1une {
namespace graph {

template <class T>
std::optional<std::vector<int>> topological_sort(const Graph<T>& g) {
    int n = g.size();
    std::vector<int> indeg(n, 0);
    for (int v = 0; v < n; v++) {
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            indeg[e.to]++;
        }
    }

    std::queue<int> que;
    for (int v = 0; v < n; v++) {
        if (indeg[v] == 0) que.push(v);
    }

    std::vector<int> order;
    order.reserve(n);
    while (!que.empty()) {
        int v = que.front();
        que.pop();
        order.push_back(v);
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            indeg[e.to]--;
            if (indeg[e.to] == 0) que.push(e.to);
        }
    }

    if (int(order.size()) != n) return std::nullopt;
    return order;
}

template <class T>
bool is_dag(const Graph<T>& g) {
    return topological_sort(g).has_value();
}

}  // namespace graph
}  // namespace m1une


#line 12 "graph/dag_shortest_path.hpp"

namespace m1une {
namespace graph {

template <class T>
struct DagShortestPathResult {
    std::vector<T> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> topological_order;
    T inf;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
std::optional<DagShortestPathResult<T>> dag_shortest_path(
    const Graph<T>& g, const std::vector<int>& sources, T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    auto order = topological_sort(g);
    if (!order) return std::nullopt;

    DagShortestPathResult<T> result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.topological_order = *order;
    result.inf = inf;

    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] == T(0)) continue;
        result.dist[s] = T(0);
    }

    for (int v : *order) {
        if (result.dist[v] == inf) continue;
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            T nd = result.dist[v] + e.cost;
            if (result.dist[e.to] <= nd) continue;
            result.dist[e.to] = nd;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
        }
    }

    return result;
}

template <class T>
std::optional<DagShortestPathResult<T>> dag_shortest_path(
    const Graph<T>& g, int s, T inf = std::numeric_limits<T>::max() / T(4)) {
    return dag_shortest_path(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/dijkstra.hpp"



#line 8 "graph/dijkstra.hpp"

#line 10 "graph/dijkstra.hpp"

namespace m1une {
namespace graph {

template <class T>
struct DijkstraResult {
    std::vector<T> dist;
    std::vector<char> reached;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    T inf = T();

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return reached[v];
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

namespace internal {

template <class T>
class DijkstraHeap {
   private:
    const std::vector<T>& dist_;
    std::vector<int> heap_;
    std::vector<int> position_;

    bool less(int first, int second) const {
        return dist_[heap_[first]] < dist_[heap_[second]];
    }

    void swap_nodes(int first, int second) {
        std::swap(heap_[first], heap_[second]);
        position_[heap_[first]] = first;
        position_[heap_[second]] = second;
    }

    void sift_up(int index) {
        while (index != 0) {
            const int parent = (index - 1) / 2;
            if (!less(index, parent)) break;
            swap_nodes(index, parent);
            index = parent;
        }
    }

    void sift_down(int index) {
        while (2 * index + 1 < int(heap_.size())) {
            int child = 2 * index + 1;
            if (child + 1 < int(heap_.size()) && less(child + 1, child)) {
                ++child;
            }
            if (!less(child, index)) break;
            swap_nodes(index, child);
            index = child;
        }
    }

   public:
    DijkstraHeap(const std::vector<T>& dist, int size)
        : dist_(dist), position_(size, -1) {
        heap_.reserve(size);
    }

    bool empty() const {
        return heap_.empty();
    }

    void push_or_decrease(int vertex) {
        int& position = position_[vertex];
        if (position == -1) {
            position = int(heap_.size());
            heap_.push_back(vertex);
        }
        sift_up(position);
    }

    int pop_min() {
        const int result = heap_.front();
        position_[result] = -1;
        if (heap_.size() == 1) {
            heap_.pop_back();
            return result;
        }
        heap_.front() = heap_.back();
        position_[heap_.front()] = 0;
        heap_.pop_back();
        sift_down(0);
        return result;
    }
};

}  // namespace internal

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g,
                           const std::vector<int>& sources) {
    int n = g.size();
    DijkstraResult<T> result;
    result.dist.resize(n);
    result.reached.assign(n, false);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);

    internal::DijkstraHeap<T> que(result.dist, n);
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.reached[s]) continue;
        result.reached[s] = true;
        result.dist[s] = T();
        que.push_or_decrease(s);
    }

    while (!que.empty()) {
        const int current = que.pop_min();
        for (const auto& e : g[current]) {
            if (!e.alive) continue;
            T nd = result.dist[current] + e.cost;
            if (result.reached[e.to] && !(nd < result.dist[e.to])) continue;
            result.reached[e.to] = true;
            result.dist[e.to] = std::move(nd);
            result.parent[e.to] = current;
            result.parent_edge[e.to] = e.id;
            que.push_or_decrease(e.to);
        }
    }

    return result;
}

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g, int s) {
    return dijkstra(g, std::vector<int>{s});
}

// Compatibility overload: unreachable distances are replaced by inf after the
// search. Reachability itself never depends on this sentinel.
template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g,
                           const std::vector<int>& sources, const T& inf) {
    DijkstraResult<T> result = dijkstra(g, sources);
    result.inf = inf;
    for (int v = 0; v < int(result.dist.size()); v++) {
        if (!result.reachable(v)) result.dist[v] = inf;
    }
    return result;
}

template <class T>
DijkstraResult<T> dijkstra(const Graph<T>& g, int s, const T& inf) {
    return dijkstra(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/k_shortest_walk.hpp"



#line 10 "graph/k_shortest_walk.hpp"

#line 12 "graph/k_shortest_walk.hpp"

namespace m1une {
namespace graph {

namespace internal {

template <class T>
class KShortestWalkHeap {
    struct Node {
        T key;
        int to;
        int left;
        int right;
        int rank;
    };

    std::vector<Node> _nodes;

    int rank(int root) const {
        return root == -1 ? 0 : _nodes[root].rank;
    }

   public:
    int make_node(T key, int to) {
        int result = int(_nodes.size());
        _nodes.push_back(Node{key, to, -1, -1, 1});
        return result;
    }

    int meld_mutable(int first, int second) {
        if (first == -1) return second;
        if (second == -1) return first;
        if (_nodes[second].key < _nodes[first].key) std::swap(first, second);
        _nodes[first].right = meld_mutable(_nodes[first].right, second);
        if (rank(_nodes[first].left) < rank(_nodes[first].right)) {
            std::swap(_nodes[first].left, _nodes[first].right);
        }
        _nodes[first].rank = rank(_nodes[first].right) + 1;
        return first;
    }

    int meld_persistent(int first, int second) {
        if (first == -1) return second;
        if (second == -1) return first;
        if (_nodes[second].key < _nodes[first].key) std::swap(first, second);
        int result = int(_nodes.size());
        _nodes.push_back(_nodes[first]);
        _nodes[result].right = meld_persistent(_nodes[result].right, second);
        if (rank(_nodes[result].left) < rank(_nodes[result].right)) {
            std::swap(_nodes[result].left, _nodes[result].right);
        }
        _nodes[result].rank = rank(_nodes[result].right) + 1;
        return result;
    }

    const Node& operator[](int index) const {
        return _nodes[index];
    }
};

}  // namespace internal

template <class T>
std::vector<T> k_shortest_walk(
    const Graph<T>& g,
    int s,
    int t,
    int k,
    T inf = std::numeric_limits<T>::max() / T(4)
) {
    int n = g.size();
    assert(0 <= s && s < n);
    assert(0 <= t && t < n);
    assert(0 <= k);
    if (k == 0) return {};

    struct ReverseEdge {
        int from;
        int index;
        T cost;
    };
    std::vector<std::vector<ReverseEdge>> reverse_graph(n);
    for (int from = 0; from < n; from++) {
        for (int index = 0; index < int(g[from].size()); index++) {
            const auto& edge = g[from][index];
            if (!edge.alive) continue;
            assert(T(0) <= edge.cost);
            reverse_graph[edge.to].push_back(ReverseEdge{from, index, edge.cost});
        }
    }

    std::vector<T> dist(n, inf);
    std::vector<int> tree_edge(n, -1);
    std::vector<int> order;
    order.reserve(n);
    using QueueEntry = std::pair<T, int>;
    std::priority_queue<QueueEntry, std::vector<QueueEntry>, std::greater<QueueEntry>> queue;
    dist[t] = T(0);
    queue.emplace(T(0), t);
    while (!queue.empty()) {
        auto [current_dist, vertex] = queue.top();
        queue.pop();
        if (dist[vertex] != current_dist) continue;
        order.push_back(vertex);
        for (const auto& edge : reverse_graph[vertex]) {
            T next_dist = current_dist + edge.cost;
            if (dist[edge.from] <= next_dist) continue;
            dist[edge.from] = next_dist;
            tree_edge[edge.from] = edge.index;
            queue.emplace(next_dist, edge.from);
        }
    }
    if (dist[s] == inf) return {};

    internal::KShortestWalkHeap<T> heap_pool;
    std::vector<int> local_heap(n, -1);
    for (int vertex : order) {
        for (int index = 0; index < int(g[vertex].size()); index++) {
            const auto& edge = g[vertex][index];
            if (!edge.alive || dist[edge.to] == inf || index == tree_edge[vertex]) continue;
            T extra = edge.cost + dist[edge.to] - dist[vertex];
            assert(T(0) <= extra);
            int node = heap_pool.make_node(extra, edge.to);
            local_heap[vertex] = heap_pool.meld_mutable(local_heap[vertex], node);
        }
    }

    std::vector<int> path_heap(n, -1);
    for (int vertex : order) {
        int inherited = -1;
        if (tree_edge[vertex] != -1) inherited = path_heap[g[vertex][tree_edge[vertex]].to];
        path_heap[vertex] = heap_pool.meld_persistent(inherited, local_heap[vertex]);
    }

    std::vector<T> result;
    result.reserve(k);
    result.push_back(dist[s]);
    std::priority_queue<QueueEntry, std::vector<QueueEntry>, std::greater<QueueEntry>> candidates;
    if (path_heap[s] != -1) {
        candidates.emplace(dist[s] + heap_pool[path_heap[s]].key, path_heap[s]);
    }
    while (int(result.size()) < k && !candidates.empty()) {
        auto [cost, node_index] = candidates.top();
        candidates.pop();
        result.push_back(cost);
        const auto& node = heap_pool[node_index];
        if (node.left != -1) {
            candidates.emplace(cost - node.key + heap_pool[node.left].key, node.left);
        }
        if (node.right != -1) {
            candidates.emplace(cost - node.key + heap_pool[node.right].key, node.right);
        }
        int next_heap = path_heap[node.to];
        if (next_heap != -1) {
            candidates.emplace(cost + heap_pool[next_heap].key, next_heap);
        }
    }
    return result;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/warshall_floyd.hpp"



#line 8 "graph/warshall_floyd.hpp"

#line 10 "graph/warshall_floyd.hpp"

namespace m1une {
namespace graph {

template <class T>
std::vector<std::vector<T>> warshall_floyd(std::vector<std::vector<T>> dist,
                                           T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    for (int k = 0; k < n; k++) {
        for (int i = 0; i < n; i++) {
            if (dist[i][k] == inf) continue;
            for (int j = 0; j < n; j++) {
                if (dist[k][j] == inf) continue;
                T nd = dist[i][k] + dist[k][j];
                if (nd < dist[i][j]) dist[i][j] = nd;
            }
        }
    }
    return dist;
}

template <class T>
std::vector<std::vector<T>> warshall_floyd(const Graph<T>& g, T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = g.size();
    std::vector<std::vector<T>> dist(n, std::vector<T>(n, inf));
    for (int i = 0; i < n; i++) dist[i][i] = T(0);
    for (int v = 0; v < n; v++) {
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (e.cost < dist[e.from][e.to]) dist[e.from][e.to] = e.cost;
        }
    }
    return warshall_floyd(std::move(dist), inf);
}

template <class T>
bool warshall_floyd_add_directed_edge(std::vector<std::vector<T>>& dist, int from, int to, T cost,
                                      T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    assert(0 <= from && from < n);
    assert(0 <= to && to < n);

    std::vector<T> to_from(n), from_to(n);
    for (int i = 0; i < n; i++) {
        to_from[i] = dist[i][from];
        from_to[i] = dist[to][i];
    }

    bool updated = false;
    for (int i = 0; i < n; i++) {
        if (to_from[i] == inf) continue;
        for (int j = 0; j < n; j++) {
            if (from_to[j] == inf) continue;
            T nd = to_from[i] + cost + from_to[j];
            if (nd < dist[i][j]) {
                dist[i][j] = nd;
                updated = true;
            }
        }
    }
    return updated;
}

template <class T>
bool warshall_floyd_add_undirected_edge(std::vector<std::vector<T>>& dist, int u, int v, T cost,
                                        T inf = std::numeric_limits<T>::max() / T(4)) {
    int n = int(dist.size());
    assert(0 <= u && u < n);
    assert(0 <= v && v < n);

    std::vector<T> to_u(n), from_u(n), to_v(n), from_v(n);
    for (int i = 0; i < n; i++) {
        to_u[i] = dist[i][u];
        from_u[i] = dist[u][i];
        to_v[i] = dist[i][v];
        from_v[i] = dist[v][i];
    }

    bool updated = false;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (to_u[i] != inf && from_v[j] != inf) {
                T nd = to_u[i] + cost + from_v[j];
                if (nd < dist[i][j]) {
                    dist[i][j] = nd;
                    updated = true;
                }
            }
            if (to_v[i] != inf && from_u[j] != inf) {
                T nd = to_v[i] + cost + from_u[j];
                if (nd < dist[i][j]) {
                    dist[i][j] = nd;
                    updated = true;
                }
            }
        }
    }
    return updated;
}

template <class T>
bool has_negative_cycle(const std::vector<std::vector<T>>& dist) {
    int n = int(dist.size());
    for (int i = 0; i < n; i++) {
        if (dist[i][i] < T(0)) return true;
    }
    return false;
}

}  // namespace graph
}  // namespace m1une


#line 1 "graph/zero_one_bfs.hpp"



#line 6 "graph/zero_one_bfs.hpp"
#include <deque>
#line 9 "graph/zero_one_bfs.hpp"

#line 11 "graph/zero_one_bfs.hpp"

namespace m1une {
namespace graph {

struct ZeroOneBfsResult {
    std::vector<int> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    int inf;

    bool reachable(int v) const {
        assert(0 <= v && v < int(dist.size()));
        return dist[v] != inf;
    }

    std::vector<int> path(int t) const {
        assert(reachable(t));
        std::vector<int> result;
        for (int v = t; v != -1; v = parent[v]) result.push_back(v);
        std::reverse(result.begin(), result.end());
        return result;
    }
};

template <class T>
ZeroOneBfsResult zero_one_bfs(const Graph<T>& g, const std::vector<int>& sources,
                              int inf = std::numeric_limits<int>::max() / 2) {
    int n = g.size();
    ZeroOneBfsResult result;
    result.dist.assign(n, inf);
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);
    result.inf = inf;

    std::deque<int> deq;
    for (int s : sources) {
        assert(0 <= s && s < n);
        if (result.dist[s] == 0) continue;
        result.dist[s] = 0;
        deq.push_back(s);
    }

    while (!deq.empty()) {
        int v = deq.front();
        deq.pop_front();
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            int w;
            if (e.cost == T(0)) {
                w = 0;
            } else {
                assert(e.cost == T(1));
                w = 1;
            }
            int nd = result.dist[v] + w;
            if (result.dist[e.to] <= nd) continue;
            result.dist[e.to] = nd;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            if (w == 0) {
                deq.push_front(e.to);
            } else {
                deq.push_back(e.to);
            }
        }
    }

    return result;
}

template <class T>
ZeroOneBfsResult zero_one_bfs(const Graph<T>& g, int s, int inf = std::numeric_limits<int>::max() / 2) {
    return zero_one_bfs(g, std::vector<int>{s}, inf);
}

}  // namespace graph
}  // namespace m1une


#line 12 "graph/shortest_path.hpp"
Back to top page