1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#pragma once


#include <bit><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.


#include "snippet/iterations.hpp"

#include "internal/dev_env.hpp"
#include "internal/types.hpp"

#include "numeric/bit.hpp"
#include "adaptor/valarray.hpp"

#include "structure/graph.hpp"


namespace uni {


template<class Graph>
struct lowest_common_ancestor {
    using size_type = internal::size_t;
    using graph_type = Graph;
    using edge_type = typename graph_type::edge_type;
    using cost_type = typename graph_type::cost_type;

    valarray<valarray<size_type>> parent;
    valarray<size_type> depth;
    valarray<cost_type> cost;

  private:
    void dfs(const graph_type &G, const edge_type& e) noexcept(NO_EXCEPT) {
        this->parent[0][e.to] = e.from;
        if(e.from >= 0) {
            this->depth[e.to] = this->depth[e.from] + 1;
            this->cost[e.to] = this->cost[e.from] + e.cost;
        }
        ITR(f, G[e.to]) {
            if(f.to != e.from) dfs(G, f);
        }
    }

  public:
    lowest_common_ancestor(const graph_type &G, const size_type root = 0) noexcept(NO_EXCEPT) { this->init(G, root); }

    void init(const graph_type &G, const size_type root = 0) noexcept(NO_EXCEPT) {
        const size_type n = static_cast<size_type>(G.size());
        const size_type d = std::bit_width<std::make_unsigned_t<size_type>>(n);

        this->parent.assign(d, valarray<size_type>(n, -1));
        this->depth.assign(n, 0), this->cost.assign(n, 0);

        this->dfs(G, edge_type(-1, root, 0));

        REP(k, d-1) REP(v, n) {
            if(this->parent[k][v] < 0) this->parent[k+1][v] = -1;
            else this->parent[k+1][v] = this->parent[k][this->parent[k][v]];
        }
    }

    size_type operator()(const size_type u, const size_type v) const noexcept(NO_EXCEPT) {
        return this->find(u, v);
    }

    size_type find(size_type u, size_type v) const noexcept(NO_EXCEPT) {
        if(this->depth[u] < this->depth[v]) std::swap(u, v);
        size_type d = static_cast<size_type>(this->parent.size());

        REP(k, d) {
            if((this->depth[u] - this->depth[v]) >> k & 1) u = this->parent[k][u];
        }

        if(u == v) return u;

        REPD(k, d) {<--- Unsigned positive
            if(this->parent[k][u] != this->parent[k][v]) {
                u = this->parent[k][u];
                v = this->parent[k][v];
            }
        }

        return this->parent[0][u];
    }

    size_type path_length(const size_type u, const size_type v) const noexcept(NO_EXCEPT) {
        return this->depth[u] + this->depth[v] - 2 * this->depth[find(u, v)];
    }

    size_type distance(const size_type u, const size_type v) const noexcept(NO_EXCEPT) {
        return this->cost[u] + this->cost[v] - 2 * this->cost[find(u, v)];
    }
};


} // namespace uni