Functional Graph
(library/graph/functional_graph.hpp)

View this file on GitHub
Last update: 2022-06-05 19:58:08+09:00
Include: #include "library/graph/functional_graph.hpp"
Link: https://noshi91.hatenablog.com/entry/2019/09/22/114149

Functional Graph

全ての頂点の出次数がちょうど $1$ のグラフ。$[N] := \{0, \ldots, N - 1\}$ としたとき、$N$ 頂点の Functional Graph の集合と関数 $f : [N] \to [N]$ の集合は 1 対 1 対応する。

`Functional Graph` がサポートする関数

std::vector<int> kth_iterate(long long k) const

$f$ の $k$ 回合成 $\underbrace{f\circ f \circ \cdots \circ f}_{k}$ を $f^k$ として、各 $i\in[N]$ に対する $f ^ k(i)$ を計算する。詳細は https://noshi91.hatenablog.com/entry/2019/09/22/114149 を参照するとよい。

計算量: $O(N)$
std::vector<InfinitePath> infinite_paths() const

Functional Graph の特徴として、任意の $i\in [N]$ に対してある非負整数 $k _ i$ と正整数 $t _ i$ が存在して、$p \geq k _ i$ ならば $f ^ p(i) = f ^ {p + t _ i}(i)$ を満たす (即ち、ある時点からループに入る)。そのような $k _ i, t _ i$ のうち最小のものを $k _ i’, t _ i’$ とする。本関数では、全ての $i$ に対する $k _ i’, t _ i’$ および $f ^ {k _ i’}(i)$ を計算する。

返り値は std::vector<InfinitePath> であるが、InfinitePath は以下のメンバを持つ。
- head_v : $i$ を表す
- heav_len : $k _ i’$ を表す
- loop_v : $f ^ {k _ i’}(i)$ を表す
- loop_len : $t _ i’$ を表す
計算量: $O(N)$
Doubling doubling(long long max_step) const

ダブリングテーブルを構築する。引数の max_step はダブリングにおける最大ステップ数。返り値の Doubling については後述。

計算量 : $O(N \log \max \mathrm{step})$
template <typename T, T(*op)(T, T), T(*e)()> DoublingSum<T, op, e> doubling(long long max_step, const std::vector<T>& dat) const

パス上の辺の値を順番通りに fold した値を計算できるダブリングテーブルを構築する。

テンプレート引数
- T : 辺に書かれた値の型
- op : 二項演算
- e : 単位元 (を与える関数)
引数
- max_step: ダブリングにおける最大ステップ数
- dat: 各辺 $(i,f(i))$ に書かれた値を表す std::vector
ただし、$(T,op)$ がモノイドであることを要求する。返り値の型 DoublingSum については後述。

計算量 : $O(N \log \max \mathrm{step})$

`Doubling` がサポートする関数

int query(int u, long long d) const

$f ^ d(u)$ を計算して返す。

計算量: $O(\log d)$
template <typename Pred> BinarySearchResult max_step(int u, Pred &&pred) const

引数に与えた $\mathrm{pred} : [N] \to \{\mathrm{true},\mathrm{false}\}$ に対して、$x := \max \{d \mid \mathrm{pred}(f ^ d (u)) = \mathrm{true}\}$ を計算する。$x$ が doubling 関数に与えた max_step 以上であるような場合は、max_step 以上の不定値を返す。

注意点として、$\mathrm{pred}(u) = \mathrm{true}$ および $\mathrm{pred}(f ^ d(u)) = \mathrm{false} \Rightarrow \mathrm{pred}(f ^ {d + 1}(u)) = \mathrm{false}$ が成り立つこと (単調性) を要求する。

返り値の BinarySearchResult は以下のメンバを持つ。
- v : $f ^ x(u)$ を表す。
- step : $x$ を表す。$x \geq \max \mathrm{step}$ の場合は、max_step 以上の不定値となる。
計算量: $O(T \cdot \log \max \mathrm{step})$ ($T$ は $\mathrm{pred}$ の計算時間)
template <typename Pred> std::optional<BinarySearchResult> step_until(int u, Pred &&f) const

引数に与えた $\mathrm{pred} : [N] \to \{\mathrm{true},\mathrm{false}\}$ に対して、$x := \min \{d \mid \mathrm{pred}(f ^ d (u)) = \mathrm{true}\}$ を計算する。$x$ が doubling 関数に与えた max_step 以上であるような場合は、std::nullopt を返す。

注意点として、$\mathrm{pred}(f ^ d(u)) = \mathrm{true} \Rightarrow \mathrm{pred}(f ^ {d + 1}(u)) = \mathrm{true}$ が成り立つこと (単調性) を要求する。特に、$\mathrm{pred}(u) = \mathrm{false}$ でなくてもよい。

計算量: $O(T \cdot \log \max \mathrm{step})$ ($T$ は $\mathrm{pred}$ の計算時間)

`DoublingSum` がサポートする関数

辺 $(i, f(i))\in [N]$ に書かれた値を $A _ i$ と書く。また、$op(A _ {i _ 1}, A _ {i _ 2})$ を $A _ {i _ 1} \odot A _ {i _ 2}$ と書く。さらに、$\displaystyle \bigodot _ {j = 1} ^ k A _ {i _ j}$ は $A _ {i _ 1} \odot A _ {i _ 2} \odot \cdots \odot A _ {i _ k}$ を表す。また、単位元を $e$ と書く。

Result query(int u, long long d) const

$f ^ d(u)$ および $\displaystyle \bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}$ を計算する。返り値の Result は以下のメンバを持つ。
- v : $f ^ d(u)$ を表す。
- sum : $\displaystyle \bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}$ を表す。
計算量: $O(T\cdot \log d)$ ($T$ は二項演算 $\odot$ の計算量)
template <typename Pred> BinarySearchResult max_step(int u, Pred &&pred) const

引数に与えた $\mathrm{pred} : T \to \{\mathrm{true},\mathrm{false}\}$ ($T$ は $A _ i$ の型) に対して、$x := \max \{d \mid \mathrm{pred}(\bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}) = \mathrm{true}\}$ を計算する。$x$ が doubling 関数に与えた max_step 以上であるような場合は、max_step 以上の不定値を返す。

注意点として、$\mathrm{pred}(e) = \mathrm{true}$ および $\mathrm{pred}(\bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}) = \mathrm{false} \Rightarrow \mathrm{pred}(\bigodot _ {i = 0} ^ {d} A _ {f ^ i(u)}) = \mathrm{false}$ が成り立つこと (単調性) を要求する。

返り値の BinarySearchResult は以下のメンバを持つ。
- v : $f ^ x(u)$ を表す。
- sum : $\displaystyle \bigodot _ {i = 0} ^ {x - 1} A _ {f ^ i(u)}$ を表す。
- step : $x$ を表す。$x \geq \max \mathrm{step}$ の場合は、max_step 以上の不定値となる。
template <typename Pred> std::optional<BinarySearchResult> step_until(int u, Pred &&f) const

引数に与えた $\mathrm{pred} : T \to \{\mathrm{true},\mathrm{false}\}$ に対して、$x := \min \{d \mid \mathrm{pred}(\bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}) = \mathrm{true}\}$ を計算する。$x$ が doubling 関数に与えた max_step 以上であるような場合は、std::nullopt を返す。

注意点として、$\mathrm{pred}(\bigodot _ {i = 0} ^ {d - 1} A _ {f ^ i(u)}) = \mathrm{true} \Rightarrow \mathrm{pred}(\bigodot _ {i = 0} ^ {d} A _ {f ^ i(u)}) = \mathrm{true}$ が成り立つこと (単調性) を要求する。特に、$\mathrm{pred}(e) = \mathrm{false}$ でなくてもよい。

Verified with

Code

#ifndef SUISEN_DOUBLING
#define SUISEN_DOUBLING

#include <cassert>
#include <cstdint>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>

namespace suisen {
    struct FunctionalGraph {
        struct Doubling;
        template <typename T, T(*)(T, T), T(*)()>
        struct DoublingSum;
        friend struct Doubling;
        template <typename T, T(*op)(T, T), T(*e)()>
        friend struct DoublingSum;

        FunctionalGraph() : FunctionalGraph(0) {}
        FunctionalGraph(int n) : _n(n), _nxt(n) {}
        FunctionalGraph(const std::vector<int>& nxt) : _n(nxt.size()), _nxt(nxt) {}

        const int& operator[](int u) const {
            return _nxt[u];
        }
        int& operator[](int u) {
            return _nxt[u];
        }

        struct Doubling {
            friend struct FunctionalGraph;

            int query(int u, long long d) const {
                for (int l = _log; l >= 0; --l) if ((d >> l) & 1) u = _nxt[l][u];
                return u;
            }

            struct BinarySearchResult {
                int v;
                long long step;
                operator std::pair<int, long long>() const { return std::pair<int, long long>{ v, step }; }
            };

            template <typename Pred>
            auto max_step(int u, Pred &&f) const {
                assert(f(u));
                long long step = 0;
                for (int l = _log; l >= 0; --l) if (int nxt_u = _nxt[l][u]; f(nxt_u)) {
                    u = nxt_u, step |= 1LL << l;
                }
                return BinarySearchResult{ u, step };
            }

            template <typename Pred>
            std::optional<BinarySearchResult> step_until(int u, Pred &&f) const {
                if (f(u)) return BinarySearchResult { u, 0 };
                auto [v, step] = max_step(u, [&](int v) { return not f(v); });
                v = _nxt[0][v], ++step;
                if (not f(v)) return std::nullopt;
                return BinarySearchResult{ v, step };
            }

        private:
            int _n, _log;
            std::vector<std::vector<int>> _nxt;

            Doubling(const std::vector<int>& nxt, long long max_step) : _n(nxt.size()), _log(floor_log2(max_step)), _nxt(_log + 1, std::vector<int>(_n)) {
                _nxt[0] = nxt;
                for (int i = 1; i <= _log; ++i) for (int j = 0; j < _n; ++j) {
                    _nxt[i][j] = _nxt[i - 1][_nxt[i - 1][j]];
                }
            }
        };

        template <typename T, T(*op)(T, T), T(*e)()>
        struct DoublingSum : private Doubling {
            friend struct FunctionalGraph;

            struct Result {
                int v;
                T sum;
                operator std::pair<int, T>() const { return std::pair<int, T>{ v, sum }; }
            };

            auto query(int u, long long d) const {
                T sum = e();
                for (int l = _log; l >= 0; --l) if ((d >> l) & 1) sum = op(sum, _dat[l][std::exchange(u, _nxt[l][u])]);
                return Result{ u, sum };
            }

            struct BinarySearchResult {
                int v;
                T sum;
                long long step;
                operator std::tuple<int, T, long long>() const { return std::tuple<int, T, long long>{ v, sum, step }; }
            };

            template <typename Pred>
            auto max_step(int u, Pred &&f) const {
                assert(f(e()));
                long long step = 0;
                T sum = e();
                for (int l = _log; l >= 0; --l) {
                    if (T nxt_sum = op(sum, _dat[l][u]); f(nxt_sum)) {
                        sum = std::move(nxt_sum), u = _nxt[l][u], step |= 1LL << l;
                    }
                }
                return BinarySearchResult{ u, sum, step };
            }
            template <typename Pred>
            std::optional<BinarySearchResult> step_until(int u, Pred &&f) const {
                if (f(e())) return BinarySearchResult { u, e(), 0 };
                auto [v, sum, step] = max_step(u, [&](const T& v) { return not f(v); });
                sum = op(sum, _dat[0][v]), v = _nxt[0][v], ++step;
                if (not f(sum)) return std::nullopt;
                return BinarySearchResult{ v, sum, step };
            }

        private:
            std::vector<std::vector<T>> _dat;

            DoublingSum(const std::vector<int>& nxt, long long max_step, const std::vector<T>& dat) : Doubling(nxt, max_step), _dat(_log + 1, std::vector<T>(_n, e())) {
                _dat[0] = dat;
                for (int i = 1; i <= _log; ++i) for (int j = 0; j < _n; ++j) {
                    _dat[i][j] = op(_dat[i - 1][j], _dat[i - 1][_nxt[i - 1][j]]);
                }
            }
        };

        Doubling doubling(long long max_step) const {
            return Doubling(_nxt, max_step);
        }

        template <typename T, T(*op)(T, T), T(*e)()>
        DoublingSum<T, op, e> doubling(long long max_step, const std::vector<T>& dat) const {
            return DoublingSum<T, op, e>(_nxt, max_step, dat);
        }

        struct InfinitePath {
            int head_v;
            int head_len;
            int loop_v;
            int loop_len;
            InfinitePath() = default;
            InfinitePath(int head_v, int head_len, int loop_v, int loop_len) : head_v(head_v), head_len(head_len), loop_v(loop_v), loop_len(loop_len) {}
        };

        std::vector<InfinitePath> infinite_paths() const {
            std::vector<InfinitePath> res(_n);

            std::vector<int> vis(_n, _n);
            std::vector<int> dep(_n, 0);

            int time = 0;
            auto dfs = [&](auto dfs, int u) -> int {
                vis[u] = time;
                int v = _nxt[u];
                if (vis[v] == vis[u]) { // found cycle
                    int loop_len = dep[u] - dep[v] + 1;
                    res[u] = { u, 0, u, loop_len };
                    return loop_len - 1;
                } else if (vis[v] < vis[u]) {
                    res[u] = { u, res[v].head_len + 1, res[v].loop_v, res[v].loop_len };
                    return 0;
                } else {
                    dep[v] = dep[u] + 1;
                    int c = dfs(dfs, v);
                    if (c > 0) { // in cycle
                        res[u] = { u, 0, u, res[v].loop_len };
                        return c - 1;
                    } else { // out of cycle
                        res[u] = { u, res[v].head_len + 1, res[v].loop_v, res[v].loop_len };
                        return 0;
                    }
                }
            };
            for (int i = 0; i < _n; ++i, ++time) if (vis[i] == _n) dfs(dfs, i);
            return res;
        }

        /**
         * Calculates k'th iterate: f(f(f(...f(i)))) for all 0 <= i < N in O(N) time.
         * Reference: https://noshi91.hatenablog.com/entry/2019/09/22/114149
         */
        std::vector<int> kth_iterate(const long long k) const {
            assert(k >= 0);
            std::vector<int> res(_n);
            std::vector<int> forest_roots;
            std::vector<std::vector<int>> forest(_n);
            std::vector<std::vector<std::pair<long long, int>>> qs(_n);
            for (const auto& path : infinite_paths()) {
                const int v = path.head_v;
                (path.head_len == 0 ? forest_roots : forest[_nxt[v]]).push_back(v);
                if (path.head_len >= k) continue;
                qs[path.loop_v].emplace_back(k - path.head_len, v);
            }
            std::vector<int> dfs_path(_n);
            auto dfs = [&](auto dfs, int u, int d) -> void {
                dfs_path[d] = u;
                if (d >= k) res[u] = dfs_path[d - k];
                for (int v : forest[u]) dfs(dfs, v, d + 1);
            };
            for (int root : forest_roots) dfs(dfs, root, 0);
            std::vector<int8_t> seen(_n, false);
            for (int root : forest_roots) {
                if (seen[root]) continue;
                std::vector<int> cycle{ root };
                for (int v = _nxt[root]; v != root; v = _nxt[v]) cycle.push_back(v);
                const int len = cycle.size();
                for (int i = 0; i < len; ++i) {
                    const int s = cycle[i];
                    seen[s] = true;
                    for (const auto& [rem, res_index] : qs[s]) {
                        res[res_index] = cycle[(i + rem) % len];
                    }
                }
            }
            return res;
        }

    private:
        int _n;
        std::vector<int> _nxt;

        static int floor_log2(long long v) {
            int l = 0;
            while (1LL << (l + 1) <= v) ++l;
            return l;
        }
    };
} // namespace suisen


#endif // SUISEN_DOUBLING

#line 1 "library/graph/functional_graph.hpp"



#include <cassert>
#include <cstdint>
#include <optional>
#include <tuple>
#include <utility>
#include <vector>

namespace suisen {
    struct FunctionalGraph {
        struct Doubling;
        template <typename T, T(*)(T, T), T(*)()>
        struct DoublingSum;
        friend struct Doubling;
        template <typename T, T(*op)(T, T), T(*e)()>
        friend struct DoublingSum;

        FunctionalGraph() : FunctionalGraph(0) {}
        FunctionalGraph(int n) : _n(n), _nxt(n) {}
        FunctionalGraph(const std::vector<int>& nxt) : _n(nxt.size()), _nxt(nxt) {}

        const int& operator[](int u) const {
            return _nxt[u];
        }
        int& operator[](int u) {
            return _nxt[u];
        }

        struct Doubling {
            friend struct FunctionalGraph;

            int query(int u, long long d) const {
                for (int l = _log; l >= 0; --l) if ((d >> l) & 1) u = _nxt[l][u];
                return u;
            }

            struct BinarySearchResult {
                int v;
                long long step;
                operator std::pair<int, long long>() const { return std::pair<int, long long>{ v, step }; }
            };

            template <typename Pred>
            auto max_step(int u, Pred &&f) const {
                assert(f(u));
                long long step = 0;
                for (int l = _log; l >= 0; --l) if (int nxt_u = _nxt[l][u]; f(nxt_u)) {
                    u = nxt_u, step |= 1LL << l;
                }
                return BinarySearchResult{ u, step };
            }

            template <typename Pred>
            std::optional<BinarySearchResult> step_until(int u, Pred &&f) const {
                if (f(u)) return BinarySearchResult { u, 0 };
                auto [v, step] = max_step(u, [&](int v) { return not f(v); });
                v = _nxt[0][v], ++step;
                if (not f(v)) return std::nullopt;
                return BinarySearchResult{ v, step };
            }

        private:
            int _n, _log;
            std::vector<std::vector<int>> _nxt;

            Doubling(const std::vector<int>& nxt, long long max_step) : _n(nxt.size()), _log(floor_log2(max_step)), _nxt(_log + 1, std::vector<int>(_n)) {
                _nxt[0] = nxt;
                for (int i = 1; i <= _log; ++i) for (int j = 0; j < _n; ++j) {
                    _nxt[i][j] = _nxt[i - 1][_nxt[i - 1][j]];
                }
            }
        };

        template <typename T, T(*op)(T, T), T(*e)()>
        struct DoublingSum : private Doubling {
            friend struct FunctionalGraph;

            struct Result {
                int v;
                T sum;
                operator std::pair<int, T>() const { return std::pair<int, T>{ v, sum }; }
            };

            auto query(int u, long long d) const {
                T sum = e();
                for (int l = _log; l >= 0; --l) if ((d >> l) & 1) sum = op(sum, _dat[l][std::exchange(u, _nxt[l][u])]);
                return Result{ u, sum };
            }

            struct BinarySearchResult {
                int v;
                T sum;
                long long step;
                operator std::tuple<int, T, long long>() const { return std::tuple<int, T, long long>{ v, sum, step }; }
            };

            template <typename Pred>
            auto max_step(int u, Pred &&f) const {
                assert(f(e()));
                long long step = 0;
                T sum = e();
                for (int l = _log; l >= 0; --l) {
                    if (T nxt_sum = op(sum, _dat[l][u]); f(nxt_sum)) {
                        sum = std::move(nxt_sum), u = _nxt[l][u], step |= 1LL << l;
                    }
                }
                return BinarySearchResult{ u, sum, step };
            }
            template <typename Pred>
            std::optional<BinarySearchResult> step_until(int u, Pred &&f) const {
                if (f(e())) return BinarySearchResult { u, e(), 0 };
                auto [v, sum, step] = max_step(u, [&](const T& v) { return not f(v); });
                sum = op(sum, _dat[0][v]), v = _nxt[0][v], ++step;
                if (not f(sum)) return std::nullopt;
                return BinarySearchResult{ v, sum, step };
            }

        private:
            std::vector<std::vector<T>> _dat;

            DoublingSum(const std::vector<int>& nxt, long long max_step, const std::vector<T>& dat) : Doubling(nxt, max_step), _dat(_log + 1, std::vector<T>(_n, e())) {
                _dat[0] = dat;
                for (int i = 1; i <= _log; ++i) for (int j = 0; j < _n; ++j) {
                    _dat[i][j] = op(_dat[i - 1][j], _dat[i - 1][_nxt[i - 1][j]]);
                }
            }
        };

        Doubling doubling(long long max_step) const {
            return Doubling(_nxt, max_step);
        }

        template <typename T, T(*op)(T, T), T(*e)()>
        DoublingSum<T, op, e> doubling(long long max_step, const std::vector<T>& dat) const {
            return DoublingSum<T, op, e>(_nxt, max_step, dat);
        }

        struct InfinitePath {
            int head_v;
            int head_len;
            int loop_v;
            int loop_len;
            InfinitePath() = default;
            InfinitePath(int head_v, int head_len, int loop_v, int loop_len) : head_v(head_v), head_len(head_len), loop_v(loop_v), loop_len(loop_len) {}
        };

        std::vector<InfinitePath> infinite_paths() const {
            std::vector<InfinitePath> res(_n);

            std::vector<int> vis(_n, _n);
            std::vector<int> dep(_n, 0);

            int time = 0;
            auto dfs = [&](auto dfs, int u) -> int {
                vis[u] = time;
                int v = _nxt[u];
                if (vis[v] == vis[u]) { // found cycle
                    int loop_len = dep[u] - dep[v] + 1;
                    res[u] = { u, 0, u, loop_len };
                    return loop_len - 1;
                } else if (vis[v] < vis[u]) {
                    res[u] = { u, res[v].head_len + 1, res[v].loop_v, res[v].loop_len };
                    return 0;
                } else {
                    dep[v] = dep[u] + 1;
                    int c = dfs(dfs, v);
                    if (c > 0) { // in cycle
                        res[u] = { u, 0, u, res[v].loop_len };
                        return c - 1;
                    } else { // out of cycle
                        res[u] = { u, res[v].head_len + 1, res[v].loop_v, res[v].loop_len };
                        return 0;
                    }
                }
            };
            for (int i = 0; i < _n; ++i, ++time) if (vis[i] == _n) dfs(dfs, i);
            return res;
        }

        /**
         * Calculates k'th iterate: f(f(f(...f(i)))) for all 0 <= i < N in O(N) time.
         * Reference: https://noshi91.hatenablog.com/entry/2019/09/22/114149
         */
        std::vector<int> kth_iterate(const long long k) const {
            assert(k >= 0);
            std::vector<int> res(_n);
            std::vector<int> forest_roots;
            std::vector<std::vector<int>> forest(_n);
            std::vector<std::vector<std::pair<long long, int>>> qs(_n);
            for (const auto& path : infinite_paths()) {
                const int v = path.head_v;
                (path.head_len == 0 ? forest_roots : forest[_nxt[v]]).push_back(v);
                if (path.head_len >= k) continue;
                qs[path.loop_v].emplace_back(k - path.head_len, v);
            }
            std::vector<int> dfs_path(_n);
            auto dfs = [&](auto dfs, int u, int d) -> void {
                dfs_path[d] = u;
                if (d >= k) res[u] = dfs_path[d - k];
                for (int v : forest[u]) dfs(dfs, v, d + 1);
            };
            for (int root : forest_roots) dfs(dfs, root, 0);
            std::vector<int8_t> seen(_n, false);
            for (int root : forest_roots) {
                if (seen[root]) continue;
                std::vector<int> cycle{ root };
                for (int v = _nxt[root]; v != root; v = _nxt[v]) cycle.push_back(v);
                const int len = cycle.size();
                for (int i = 0; i < len; ++i) {
                    const int s = cycle[i];
                    seen[s] = true;
                    for (const auto& [rem, res_index] : qs[s]) {
                        res[res_index] = cycle[(i + rem) % len];
                    }
                }
            }
            return res;
        }

    private:
        int _n;
        std::vector<int> _nxt;

        static int floor_log2(long long v) {
            int l = 0;
            while (1LL << (l + 1) <= v) ++l;
            return l;
        }
    };
} // namespace suisen

Functional Graph (library/graph/functional_graph.hpp)