m1une's library

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

View on GitHub

:heavy_check_mark: Tree All
(graph/tree/all.hpp)

Overview

graph/tree/all.hpp includes every tree header in this directory. Most algorithms use m1une::graph::Graph<T> as their input container, and sequence-derived helpers such as Cartesian tree can convert their result to the same graph container.

Included Headers

Header Contents
graph/tree/cartesian_tree.hpp Cartesian tree construction from an array.
graph/tree/tree.hpp Core rooted tree and diameter bundle.
graph/tree/cumulative_sum.hpp Static commutative-group products and additive sums on vertex- or edge-weighted tree paths.
graph/tree/euler_tour.hpp Lightweight rooted-tree preorder, subtree ranges, and parent/depth metadata.
graph/tree/rooted_tree.hpp Rooted metadata, Euler intervals, LCA, jumps, paths, and distances.
graph/tree/sparse_table_lca.hpp Euler-tour sparse-table LCA with $O(1)$ queries.
graph/tree/heavy_light_decomposition.hpp HLD order, path segments, subtree ranges, LCA, and jumps.
graph/tree/mo_on_tree.hpp Offline Mo ordering for static vertex-path or edge-path queries.
graph/tree/range_contour_query.hpp Commutative-group contour products and applications by unweighted tree distance, with additive wrappers.
graph/tree/diameter.hpp Weighted tree/forest diameter path.
graph/tree/distance_frequency.hpp Counts unordered vertex pairs at every unweighted tree distance.
graph/tree/tree_hash.hpp Probabilistic rooted-subtree and unrooted-tree isomorphism hashes.
graph/tree/dsu_on_tree.hpp Iterative small-to-large subtree processing with user callbacks.
graph/tree/rerooting_dp.hpp Generic rerooting DP helper.
graph/tree/static_top_tree.hpp Dynamic one-root tree DP on a fixed tree using static top-tree clusters.
graph/tree/rerooting_static_top_tree.hpp Bidirectional static top-tree clusters for dynamic rerooting-DP solutions.
graph/tree/centroid_decomposition.hpp Centroid decomposition for trees and forests.
graph/tree/virtual_tree.hpp Virtual-tree compression for a selected vertex set.
graph/tree/zero_one_on_tree.hpp Minimum inversions in an ancestor-respecting order of binary-labeled vertices.

Example

#include "graph/graph.hpp"
#include "graph/tree/all.hpp"
#include <iostream>

int main() {
    m1une::graph::Graph<int> g(3);
    g.add_edge(0, 1);
    g.add_edge(1, 2);

    auto diameter = m1une::tree::tree_diameter(g);
    std::cout << diameter.edge_count << "\n";
}

Depends on

Required by

Verified with

Code

#ifndef M1UNE_TREE_ALL_HPP
#define M1UNE_TREE_ALL_HPP 1

#include "cartesian_tree.hpp"
#include "centroid_decomposition.hpp"
#include "cumulative_sum.hpp"
#include "diameter.hpp"
#include "distance_frequency.hpp"
#include "dsu_on_tree.hpp"
#include "euler_tour.hpp"
#include "heavy_light_decomposition.hpp"
#include "mo_on_tree.hpp"
#include "range_contour_query.hpp"
#include "rerooting_dp.hpp"
#include "rerooting_static_top_tree.hpp"
#include "rooted_tree.hpp"
#include "sparse_table_lca.hpp"
#include "static_top_tree.hpp"
#include "tree.hpp"
#include "tree_hash.hpp"
#include "virtual_tree.hpp"
#include "zero_one_on_tree.hpp"

#endif  // M1UNE_TREE_ALL_HPP
#line 1 "graph/tree/all.hpp"



#line 1 "graph/tree/cartesian_tree.hpp"



#include <cassert>
#include <cstddef>
#include <functional>
#include <limits>
#include <utility>
#include <vector>

#line 1 "graph/graph.hpp"



#include <array>
#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 12 "graph/tree/cartesian_tree.hpp"

namespace m1une {
namespace tree {

struct CartesianTree {
    int root;
    std::vector<int> parent;
    std::vector<int> left;
    std::vector<int> right;

   private:
    int _n;

    void check_vertex(int v) const {
        assert(0 <= v && v < _n);
    }

   public:
    CartesianTree() : root(-1), _n(0) {}

    template <class T, class Compare = std::less<T>>
    explicit CartesianTree(const std::vector<T>& a, Compare comp = Compare()) : root(-1), _n(0) {
        build(a, comp);
    }

    template <class T, class Compare = std::less<T>>
    void build(const std::vector<T>& a, Compare comp = Compare()) {
        assert(a.size() <= static_cast<std::size_t>(std::numeric_limits<int>::max()));
        _n = int(a.size());
        root = -1;
        parent.assign(_n, -1);
        left.assign(_n, -1);
        right.assign(_n, -1);

        std::vector<int> stack;
        stack.reserve(_n);
        for (int i = 0; i < _n; i++) {
            int last = -1;
            while (!stack.empty() && comp(a[i], a[stack.back()])) {
                last = stack.back();
                stack.pop_back();
            }
            if (last != -1) {
                left[i] = last;
                parent[last] = i;
            }
            if (!stack.empty()) {
                right[stack.back()] = i;
                parent[i] = stack.back();
            }
            stack.push_back(i);
        }

        if (!stack.empty()) root = stack.front();
    }

    int size() const {
        return _n;
    }

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

    int parent_or_self(int v) const {
        check_vertex(v);
        return parent[v] == -1 ? v : parent[v];
    }

    std::vector<int> parent_with_root_self() const {
        std::vector<int> result = parent;
        if (root != -1) result[root] = root;
        return result;
    }

    std::vector<std::pair<int, int>> edges() const {
        std::vector<std::pair<int, int>> result;
        if (_n == 0) return result;
        result.reserve(_n - 1);
        for (int v = 0; v < _n; v++) {
            if (parent[v] != -1) result.emplace_back(parent[v], v);
        }
        return result;
    }

    m1une::graph::Graph<int> to_graph() const {
        m1une::graph::Graph<int> g(_n);
        for (int v = 0; v < _n; v++) {
            if (parent[v] != -1) g.add_edge(parent[v], v);
        }
        return g;
    }
};

template <class T, class Compare = std::less<T>>
CartesianTree cartesian_tree(const std::vector<T>& a, Compare comp = Compare()) {
    CartesianTree result;
    result.build(a, comp);
    return result;
}

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/centroid_decomposition.hpp"



#include <algorithm>
#line 6 "graph/tree/centroid_decomposition.hpp"

#line 8 "graph/tree/centroid_decomposition.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct CentroidDecomposition {
    int n;
    std::vector<int> parent;
    std::vector<int> depth;
    std::vector<int> order;
    std::vector<int> roots;
    std::vector<std::vector<int>> children;

   private:
    std::vector<int> _subtree_size;
    std::vector<int> _work_parent;
    std::vector<char> _removed;

    void build_component(const m1une::graph::Graph<T>& g, int start, int p, int d) {
        std::vector<int> nodes;
        std::vector<int> stack = {start};
        _work_parent[start] = -2;
        while (!stack.empty()) {
            int v = stack.back();
            stack.pop_back();
            nodes.push_back(v);
            for (const auto& e : g[v]) {
                if (!e.alive || _removed[e.to]) continue;
                if (_work_parent[e.to] != -1) continue;
                _work_parent[e.to] = v;
                stack.push_back(e.to);
            }
        }

        for (int v : nodes) _subtree_size[v] = 1;
        for (int i = int(nodes.size()) - 1; i >= 0; i--) {
            int v = nodes[i];
            if (_work_parent[v] >= 0) _subtree_size[_work_parent[v]] += _subtree_size[v];
        }

        int total = int(nodes.size());
        int centroid = start;
        int best = total + 1;
        for (int v : nodes) {
            int largest = total - _subtree_size[v];
            for (const auto& e : g[v]) {
                if (!e.alive || _removed[e.to]) continue;
                if (_work_parent[e.to] == v) largest = std::max(largest, _subtree_size[e.to]);
            }
            if (largest < best) {
                best = largest;
                centroid = v;
            }
        }

        for (int v : nodes) _work_parent[v] = -1;

        parent[centroid] = p;
        depth[centroid] = d;
        order.push_back(centroid);
        if (p == -1) {
            roots.push_back(centroid);
        } else {
            children[p].push_back(centroid);
        }
        _removed[centroid] = true;

        for (const auto& e : g[centroid]) {
            if (!e.alive || _removed[e.to]) continue;
            build_component(g, e.to, centroid, d + 1);
        }
    }

   public:
    CentroidDecomposition() : n(0) {}
    explicit CentroidDecomposition(const m1une::graph::Graph<T>& g) {
        build(g);
    }

    void build(const m1une::graph::Graph<T>& g) {
        n = g.size();
        parent.assign(n, -1);
        depth.assign(n, -1);
        order.clear();
        order.reserve(n);
        roots.clear();
        children.assign(n, {});
        _subtree_size.assign(n, 0);
        _work_parent.assign(n, -1);
        _removed.assign(n, false);

        for (int v = 0; v < n; v++) {
            if (depth[v] == -1) build_component(g, v, -1, 0);
        }
    }

    int size() const {
        return n;
    }

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

    int root() const {
        return roots.empty() ? -1 : roots[0];
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/cumulative_sum.hpp"



#line 8 "graph/tree/cumulative_sum.hpp"

#line 1 "monoid/add.hpp"



namespace m1une {
namespace monoid {

// Monoid for addition (Range Sum).
template <typename T>
struct Add {
    using value_type = T;
    static constexpr bool commutative = true;

    // Returns the identity element for addition, which is 0.
    static constexpr T id() {
        return T(0);
    }

    // Returns the sum of a and b.
    static constexpr T op(const T& a, const T& b) {
        return a + b;
    }

    static constexpr T inv(const T& x) {
        return -x;
    }
};

}  // namespace monoid
}  // namespace m1une


#line 1 "monoid/concept.hpp"



#include <concepts>

namespace m1une {
namespace monoid {

// Concept to check if a type satisfies the requirements of a Monoid.
// A Monoid must have a `value_type`, an identity element `id()`, and an associative binary operation `op()`.
template <typename M>
concept IsMonoid = requires(typename M::value_type a, typename M::value_type b) {
    // 1. Must define `value_type`
    typename M::value_type;

    // 2. Must have a static method `id()` returning `value_type`
    { M::id() } -> std::same_as<typename M::value_type>;

    // 3. Must have a static method `op(a, b)` returning `value_type`
    { M::op(a, b) } -> std::same_as<typename M::value_type>;
};

// Concept for groups. A type satisfying this concept must also obey the group
// laws; concepts can check the interface but not the algebraic properties.
template <typename M>
concept IsGroup = IsMonoid<M> && requires(typename M::value_type a) {
    { M::inv(a) } -> std::same_as<typename M::value_type>;
};

// Concept for commutative groups. Commutativity is a semantic requirement and
// cannot be checked by a C++ concept.
template <typename M>
concept IsCommutativeGroup = IsGroup<M>;

}  // namespace monoid
}  // namespace m1une


#line 12 "graph/tree/cumulative_sum.hpp"

namespace m1une {
namespace tree {

// Static cumulative products on root paths. Values are attached to vertices by
// default; set EdgeValues to true to index them by graph edge id instead.
template <m1une::monoid::IsCommutativeGroup Group, bool EdgeValues = false>
class TreeCumulativeProduct {
   public:
    using value_type = typename Group::value_type;

   private:
    int _n = 0;
    int _root = -1;
    std::vector<int> _parent;
    std::vector<int> _depth;
    std::vector<int> _head;
    std::vector<value_type> _prefix;

    void check_vertex(int vertex) const {
        assert(0 <= vertex && vertex < _n);
    }

   public:
    TreeCumulativeProduct() = default;

    template <class EdgeCost>
    explicit TreeCumulativeProduct(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<value_type>& values,
        int root = 0
    ) {
        build(graph, values, root);
    }

    template <class EdgeCost>
    void build(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<value_type>& values,
        int root = 0
    ) {
        _n = graph.size();
        _root = _n == 0 ? -1 : root;
        assert(
            int(values.size())
            == (EdgeValues ? graph.edge_count() : graph.size())
        );

        _parent.assign(_n, -2);
        _depth.assign(_n, 0);
        _head.assign(_n, -1);
        _prefix.assign(_n, Group::id());
        if (_n == 0) return;
        assert(0 <= root && root < _n);

        std::vector<int> parent_edge(_n, -1);
        std::vector<int> order;
        order.reserve(_n);
        std::vector<int> stack = {root};
        _parent[root] = -1;
        while (!stack.empty()) {
            int vertex = stack.back();
            stack.pop_back();
            order.push_back(vertex);
            for (const auto& edge : graph[vertex]) {
                if (!edge.alive || _parent[edge.to] != -2) continue;
                _parent[edge.to] = vertex;
                parent_edge[edge.to] = edge.id;
                _depth[edge.to] = _depth[vertex] + 1;
                stack.push_back(edge.to);
            }
        }
        assert(int(order.size()) == _n);

        std::vector<int> subtree_size(_n, 1);
        std::vector<int> heavy(_n, -1);
        for (int index = _n - 1; index > 0; index--) {
            int vertex = order[index];
            int parent = _parent[vertex];
            subtree_size[parent] += subtree_size[vertex];
            if (
                heavy[parent] == -1
                || subtree_size[heavy[parent]] < subtree_size[vertex]
            ) {
                heavy[parent] = vertex;
            }
        }

        std::vector<std::pair<int, int>> starts;
        starts.emplace_back(root, root);
        while (!starts.empty()) {
            auto [start, head] = starts.back();
            starts.pop_back();
            for (
                int vertex = start;
                vertex != -1;
                vertex = heavy[vertex]
            ) {
                _head[vertex] = head;
                for (const auto& edge : graph[vertex]) {
                    if (
                        edge.alive && _parent[edge.to] == vertex
                        && edge.to != heavy[vertex]
                    ) {
                        starts.emplace_back(edge.to, edge.to);
                    }
                }
            }
        }

        if constexpr (!EdgeValues) _prefix[root] = values[root];
        for (int vertex : order) {
            if (vertex == root) continue;
            if constexpr (EdgeValues) {
                assert(0 <= parent_edge[vertex]);
                _prefix[vertex] = Group::op(
                    _prefix[_parent[vertex]],
                    values[parent_edge[vertex]]
                );
            } else {
                _prefix[vertex] = Group::op(
                    _prefix[_parent[vertex]],
                    values[vertex]
                );
            }
        }
    }

    int size() const {
        return _n;
    }

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

    int root() const {
        return _root;
    }

    int lca(int first, int second) const {
        check_vertex(first);
        check_vertex(second);
        while (_head[first] != _head[second]) {
            if (_depth[_head[first]] < _depth[_head[second]]) {
                std::swap(first, second);
            }
            first = _parent[_head[first]];
        }
        return _depth[first] < _depth[second] ? first : second;
    }

    // Product on the root-to-vertex path. The root vertex is included for
    // vertex values; no edge lies above it in edge-value mode.
    value_type prod(int vertex) const {
        check_vertex(vertex);
        return _prefix[vertex];
    }

    // Product on the simple path from first to second. Both endpoints are
    // included for vertex values.
    value_type prod(int first, int second) const {
        int ancestor = lca(first, second);
        value_type result = Group::op(_prefix[first], _prefix[second]);
        result = Group::op(result, Group::inv(_prefix[ancestor]));
        if constexpr (EdgeValues) {
            result = Group::op(result, Group::inv(_prefix[ancestor]));
        } else if (_parent[ancestor] != -1) {
            result = Group::op(
                result,
                Group::inv(_prefix[_parent[ancestor]])
            );
        }
        return result;
    }
};

template <m1une::monoid::IsCommutativeGroup Group>
using TreeEdgeCumulativeProduct = TreeCumulativeProduct<Group, true>;

template <class T, bool EdgeValues = false>
class TreeCumulativeSum
    : public TreeCumulativeProduct<m1une::monoid::Add<T>, EdgeValues> {
   private:
    using Base =
        TreeCumulativeProduct<m1une::monoid::Add<T>, EdgeValues>;

   public:
    using Base::Base;

    T sum(int vertex) const {
        return Base::prod(vertex);
    }

    T sum(int first, int second) const {
        return Base::prod(first, second);
    }
};

template <class T>
using TreeEdgeCumulativeSum = TreeCumulativeSum<T, true>;

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/diameter.hpp"



#line 6 "graph/tree/diameter.hpp"

#line 8 "graph/tree/diameter.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct TreeDiameter {
    T cost;
    int edge_count;
    int from;
    int to;
    std::vector<int> vertices;
    std::vector<int> edge_ids;

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

namespace internal {

template <class T>
struct FarthestResult {
    int vertex;
    std::vector<char> seen;
    std::vector<T> dist;
    std::vector<int> parent;
    std::vector<int> parent_edge;
};

template <class T>
FarthestResult<T> farthest_from(const m1une::graph::Graph<T>& g, int start) {
    int n = g.size();
    FarthestResult<T> result;
    result.vertex = start;
    result.seen.assign(n, false);
    result.dist.assign(n, T(0));
    result.parent.assign(n, -1);
    result.parent_edge.assign(n, -1);

    std::vector<int> stack = {start};
    result.seen[start] = true;
    while (!stack.empty()) {
        int v = stack.back();
        stack.pop_back();
        if (result.dist[result.vertex] < result.dist[v]) result.vertex = v;
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (result.seen[e.to]) continue;
            result.seen[e.to] = true;
            result.dist[e.to] = result.dist[v] + e.cost;
            result.parent[e.to] = v;
            result.parent_edge[e.to] = e.id;
            stack.push_back(e.to);
        }
    }
    return result;
}

}  // namespace internal

template <class T>
TreeDiameter<T> tree_diameter(const m1une::graph::Graph<T>& g) {
    int n = g.size();
    TreeDiameter<T> best;
    best.cost = T(0);
    best.edge_count = 0;
    best.from = -1;
    best.to = -1;
    if (n == 0) return best;

    std::vector<char> done(n, false);
    for (int start = 0; start < n; start++) {
        if (done[start]) continue;
        auto first = internal::farthest_from(g, start);
        for (int v = 0; v < n; v++) {
            if (first.seen[v]) done[v] = true;
        }
        auto second = internal::farthest_from(g, first.vertex);
        int a = first.vertex;
        int b = second.vertex;
        T cost = second.dist[b];
        if (best.from != -1 && !(best.cost < cost)) continue;

        best.cost = cost;
        best.from = a;
        best.to = b;
        best.vertices.clear();
        best.edge_ids.clear();
        for (int v = b; v != -1; v = second.parent[v]) {
            best.vertices.push_back(v);
            if (v != a) best.edge_ids.push_back(second.parent_edge[v]);
        }
        std::reverse(best.vertices.begin(), best.vertices.end());
        std::reverse(best.edge_ids.begin(), best.edge_ids.end());
        best.edge_count = int(best.edge_ids.size());
    }

    return best;
}

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/distance_frequency.hpp"



#line 7 "graph/tree/distance_frequency.hpp"
#include <cstdint>
#line 10 "graph/tree/distance_frequency.hpp"

#line 1 "math/fps/convolution.hpp"



#line 8 "math/fps/convolution.hpp"
#include <cstring>
#include <new>
#include <type_traits>
#line 13 "math/fps/convolution.hpp"

#if defined(__GNUC__) && !defined(__clang__) && \
    (defined(__x86_64__) || defined(__i386__)) && \
    !defined(M1UNE_FPS_DISABLE_X86_SIMD)
#include <immintrin.h>
#define M1UNE_FPS_HAS_X86_SIMD 1
#pragma GCC push_options
#pragma GCC target("avx2,bmi")
#endif

#line 1 "math/fps/internal/ntt998_faster.hpp"



#ifdef M1UNE_FPS_HAS_X86_SIMD

#line 9 "math/fps/internal/ntt998_faster.hpp"

#include <immintrin.h>

namespace m1une {
namespace fps {
namespace internal {
namespace fast998_v2 {

// Fixed-modulus AVX2 transform with an in-register degree-8 residue product.

using u32=unsigned;
using u64=unsigned long long;
using idt=std::size_t;
using I256=__m256i;
inline void store256(void*p,I256 x){
    _mm256_store_si256((I256*)p,x);
}
inline I256 load256(const void*p){
    return _mm256_load_si256((const I256*)p);
}
constexpr u32 shrk(u32 x,u32 M){
    return std::min(x,x-M);
}
constexpr u32 dilt(u32 x,u32 M){
    return std::min(x,x+M);
}
constexpr u32 reduce(u64 x,u32 niv,u32 M){
    return (x+u64(u32(x)*niv)*M)>>32;
}
constexpr u32 mul(u32 x,u32 y,u32 niv,u32 M){
    return reduce(u64(x)*y,niv,M);
}
constexpr u32 mul_s(u32 x,u32 y,u32 niv,u32 M){
    return shrk(reduce(u64(x)*y,niv,M),M);
}
constexpr u32 qpw(u32 a,u32 b,u32 niv,u32 M,u32 r){
    for(;b;b>>=1,a=mul(a,a,niv,M)){
        if(b&1){
            r=mul(r,a,niv,M);
        }
    }
    return r;
}
constexpr u32 qpw_s(u32 a,u32 b,u32 niv,u32 M,u32 r){
    return shrk(qpw(a,b,niv,M,r),M);
}
inline I256 shrk32(I256 x,I256 M){
    return _mm256_min_epu32(x,_mm256_sub_epi32(x,M));
}
inline I256 dilt32(I256 x,I256 M){
    return _mm256_min_epu32(x,_mm256_add_epi32(x,M));
}
inline I256 Ladd32(I256 x,I256 y,I256){
    return _mm256_add_epi32(x,y);
}
inline I256 Lsub32(I256 x,I256 y,I256 M){
    return _mm256_add_epi32(_mm256_sub_epi32(x,y),M);
}
inline I256 add32(I256 x,I256 y,I256 M){
    return shrk32(_mm256_add_epi32(x,y),M);
}
inline I256 sub32(I256 x,I256 y,I256 M){
    return dilt32(_mm256_sub_epi32(x,y),M);
}
template<int msk>inline I256 neg32_m(I256 x,I256 M){
    return _mm256_blend_epi32(x,_mm256_sub_epi32(M,x),msk);
}
inline I256 reduce(I256 a,I256 b,I256 niv,I256 M){
    I256 c=_mm256_mul_epu32(a,niv),d=_mm256_mul_epu32(b,niv);
    c=_mm256_mul_epu32(c,M),d=_mm256_mul_epu32(d,M);
    return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(a,c),32),_mm256_add_epi64(b,d),0xaa);
}
inline I256 mul(I256 a,I256 b,I256 niv,I256 M){
    return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32)),niv,M);
}
inline I256 mul_s(I256 a,I256 b,I256 niv,I256 M){
    return shrk32(mul(a,b,niv,M),M);
}
inline I256 mul_bsm(I256 a,I256 b,I256 niv,I256 M){
    return reduce(_mm256_mul_epu32(a,b),_mm256_mul_epu32(_mm256_srli_epi64(a,32),b),niv,M);
}
inline I256 mul_bsmfxd(I256 a,I256 b,I256 bniv,I256 M){
    I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),bniv);
    I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),b);
    cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
    return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_bfxd(I256 a,I256 b,I256 bniv,I256 M){
    I256 cc=_mm256_mul_epu32(a,bniv),dd=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(bniv,32));
    I256 c=_mm256_mul_epu32(a,b),d=_mm256_mul_epu32(_mm256_srli_epi64(a,32),_mm256_srli_epi64(b,32));
    cc=_mm256_mul_epu32(cc,M),dd=_mm256_mul_epu32(dd,M);
    return _mm256_blend_epi32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),_mm256_add_epi64(d,dd),0xaa);
}
inline I256 mul_upd_rt(I256 a,I256 bu,I256 M){
    I256 cc=_mm256_mul_epu32(a,bu),c=_mm256_mul_epu32(a,_mm256_srli_epi64(bu,32));
    cc=_mm256_mul_epu32(cc,M);
    return shrk32(_mm256_srli_epi64(_mm256_add_epi64(c,cc),32),M);
}
constexpr auto _mxlg=26,_lg_itth=6;
constexpr auto _itth=idt(1)<<_lg_itth;
static_assert(_lg_itth%2==0);
struct FNTT32_info{
    u32 mod,mod2,niv,one,r2,r3,img,imgniv,RT1[_mxlg];
    alignas(32) std::array<u32,8> rt3[_mxlg-2],rt3i[_mxlg-2],bwbr,bwb,bwbi,rt4[_mxlg-3],rt4niv[_mxlg-3],rt4i[_mxlg-3],rt4iniv[_mxlg-3],pr2,pr4,pr2niv,pr4niv,pr2i,pr2iniv,pr4i,pr4iniv;
    constexpr FNTT32_info(const u32 m):mod(m),mod2(m*2),niv([&]{u32 n=2+m;for(int i=0;i<4;++i){n*=2+m*n;}return n;}()),one((-m)%m),r2((-u64(m))%m),r3(mul_s(r2,r2,niv,m)),img{},imgniv{},RT1{},rt3{},rt3i{},bwbr{},bwb{},bwbi{},rt4{},rt4niv{},rt4i{},rt4iniv{},pr2{},pr4{},pr2niv{},pr4niv{},pr2i{},pr2iniv{},pr4i{},pr4iniv{}{
        const int k=__builtin_ctz(m-1);
		u32 _g=mul(3,r2,niv,mod);
        for(;;++_g){
            if(qpw_s(_g,mod>>1,niv,mod,one)!=one){
                break;
            }
        }
		_g=qpw(_g,mod>>k,niv,mod,one);
        u32 rt1[_mxlg-1],rt1i[_mxlg-1];
        rt1[k-2]=_g,rt1i[k-2]=qpw(_g,mod-2,niv,mod,one);
        for(int i=k-2;i>0;--i){
            rt1[i-1]=mul(rt1[i],rt1[i],niv,mod);
            rt1i[i-1]=mul(rt1i[i],rt1i[i],niv,mod);
        }
        RT1[k-1]=qpw_s(_g,3,niv,mod,one);
        for(int i=k-1;i>0;--i){
			RT1[i-1]=mul_s(RT1[i],RT1[i],niv,mod);
        }
        img=rt1[0],imgniv=img*niv;
        bwbr={one,0,one,0,one};
        bwb={rt1[1],0,rt1[0],0,mod-mul_s(rt1[0],rt1[1],niv,mod)};
        bwbi={rt1i[1],0,rt1i[0],0,mul_s(rt1i[0],rt1i[1],niv,mod)};
        u32 pr=one,pri=one;
        for(int i=0;i<k-2;++i){
            const u32 r=mul_s(pr,rt1[i+1],niv,mod),ri=mul_s(pri,rt1i[i+1],niv,mod);
            const u32 r2=mul_s(r,r,niv,mod),r2i=mul_s(ri,ri,niv,mod);
            const u32 r3=mul_s(r,r2,niv,mod),r3i=mul_s(ri,r2i,niv,mod);
            rt3[i]={r*niv,r,r2*niv,r2,r3*niv,r3};
            rt3i[i]={ri*niv,ri,r2i*niv,r2i,r3i*niv,r3i};
            pr=mul(pr,rt1i[i+1],niv,mod),pri=mul(pri,rt1[i+1],niv,mod);
        }
        pr=one,pri=one;
        for(int i=0;i<k-3;++i){
            const u32 r=mul_s(pr,rt1[i+2],niv,mod),ri=mul_s(pri,rt1i[i+2],niv,mod);
            rt4[i][0]=rt4i[i][0]=one;
            for(int j=1;j<8;++j){
                rt4[i][j]=mul_s(rt4[i][j-1],r,niv,mod);
                rt4i[i][j]=mul_s(rt4i[i][j-1],ri,niv,mod);
            }
            for(int j=0;j<8;++j){
                rt4niv[i][j]=rt4[i][j]*niv;
                rt4iniv[i][j]=rt4i[i][j]*niv;
            }
            pr=mul(pr,rt1i[i+2],niv,mod),pri=mul(pri,rt1[i+2],niv,mod);
        }
        pr2={one,one,one,img,one,one,one,img};
        pr4={one,one,one,one,one,rt1[1],img,mul_s(img,rt1[1],niv,mod)};
        const u32 nr2=mod-r2,imgr2=mul_s(img,r2,niv,mod);
        pr2i={nr2,nr2,nr2,imgr2,nr2,nr2,nr2,imgr2};
        pr4i={one,one,one,one,one,rt1i[1],rt1i[0],mul_s(rt1i[0],rt1i[1],niv,mod)};
        for(int j=0;j<8;++j){
            pr2niv[j]=pr2[j]*niv,pr4niv[j]=pr4[j]*niv;
            pr2iniv[j]=pr2i[j]*niv,pr4iniv[j]=pr4i[j]*niv;
        }
    }
};
inline void vector_dif(I256*const f,const idt n,const FNTT32_info*info){
    alignas(32) std::array<u32,8> st_1[_mxlg>>1];
    const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
    const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
    const int lgn=__builtin_ctzll(n);
    std::fill(st_1,st_1+(lgn>>1),info->bwb);
    const idt nn=n>>(lgn&1),m=std::min(n,_itth),mm=std::min(nn,_itth);
    // I256 rr=_mm256_set1_epi32(info->one);
    if(nn!=n){
        for(idt i=0;i<nn;++i){
            auto const p0=f+i,p1=f+nn+i;
            const auto f0=load256(p0),f1=load256(p1);
            const auto g0=add32(f0,f1,Mod2),g1=Lsub32(f0,f1,Mod2);
            store256(p0,g0),store256(p1,g1);
        }
    }
    for(idt L=nn>>2;L>0;L>>=2){
        for(idt i=0;i<L;++i){
            auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
            const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
            const auto g3=mul_bsmfxd(Lsub32(f1,f3,Mod2),Img,ImgNiv,Mod),g1=add32(f1,f3,Mod2);
            const auto g0=add32(f0,f2,Mod2),g2=sub32(f0,f2,Mod2);
            const auto h0=add32(g0,g1,Mod2),h1=Lsub32(g0,g1,Mod2);
            const auto h2=Ladd32(g2,g3,Mod2),h3=Lsub32(g2,g3,Mod2);
            store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
        }
    }
    for(idt j=0;j<n;j+=m){
        int t=((j==0)?std::min(_lg_itth,lgn):__builtin_ctzll(j))&-2,p=(t-2)>>1;
        for(idt L=(idt(1)<<t)>>2;L>=_itth;L>>=2,t-=2,--p){
            auto rt=load256(st_1+p);
            const auto r1=_mm256_permutevar8x32_epi32(rt,id);
            const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
            rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~j>>t)),Mod);
            const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
            const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),nr3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
            store256(st_1+p,rt);
            for(idt i=0;i<L;++i){
                auto const p0=f+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
                const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
                const auto g1=mul_bsmfxd(f1,r1,r1Niv,Mod),ng3=mul_bsmfxd(f3,nr3,nr3Niv,Mod);
                const auto g2=mul_bsmfxd(f2,r2,r2Niv,Mod),g0=shrk32(f0,Mod2);
                const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
                const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
                const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
                const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
                store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
            }
        }
        I256*const g=f+j;
        for(idt l=mm,L=mm>>2;L;l=L,L>>=2,t-=2,--p){
            auto rt=load256(st_1+p);
            for(idt i=(j==0?l:0),k=(j+i)>>t;i<m;i+=l,++k){
                const auto r1=_mm256_permutevar8x32_epi32(rt,id);
                const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
                const auto nr3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
                for(idt j=0;j<L;++j){
                    auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
                    const auto f1=load256(p1),f3=load256(p3),f2=load256(p2),f0=load256(p0);
                    const auto g1=mul_bsm(f1,r1,Niv,Mod),ng3=mul_bsm(f3,nr3,Niv,Mod);
                    const auto g2=mul_bsm(f2,r2,Niv,Mod),g0=shrk32(f0,Mod2);
                    const auto h3=mul_bsmfxd(Ladd32(g1,ng3,Mod2),Img,ImgNiv,Mod),h1=sub32(g1,ng3,Mod2);
                    const auto h0=add32(g0,g2,Mod2),h2=sub32(g0,g2,Mod2);
                    const auto u0=Ladd32(h0,h1,Mod2),u1=Lsub32(h0,h1,Mod2);
                    const auto u2=Ladd32(h2,h3,Mod2),u3=Lsub32(h2,h3,Mod2);
                    store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
                }
                rt=mul_upd_rt(rt,load256(info->rt3+__builtin_ctzll(~k)),Mod);
            }
            store256(st_1+p,rt);
        }
        // const auto pr2=load256(&info->pr2),pr4=load256(&info->pr4);
        // const auto pr2Niv=load256(&info->pr2niv),pr4Niv=load256(&info->pr4niv);
        // for(idt i=j;i<j+m;++i){
        //     auto fi=load256(f+i);
        //     fi=mul(fi,rr,Niv,Mod);
        //     rr=shrk32(mul_bfxd(rr,load256(info->rt4+__builtin_ctzll(~i)),load256(info->rt4niv+__builtin_ctzll(~i)),Mod),Mod);
        //     fi=mul_bfxd(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),pr4,pr4Niv,Mod);
        //     fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr2,pr2Niv,Mod);
        //     fi=sub32(_mm256_shuffle_epi32(fi,0xb1),neg32_m<0x55>(fi,Mod2),Mod2);
        //     store256(f+i,fi);
        // }
    }
}
template<bool shrk=false>inline void vector_dit(I256*const f,idt n,const FNTT32_info*const info){
    alignas(32) std::array<u32,8> st_1[_mxlg>>1];
    const I256 Mod=_mm256_set1_epi32(info->mod),Mod2=_mm256_set1_epi32(info->mod2),Niv=_mm256_set1_epi32(info->niv);
    const I256 Img=_mm256_set1_epi32(info->img),ImgNiv=_mm256_set1_epi32(info->imgniv),id=_mm256_setr_epi32(0,2,0,4,0,2,0,4);
    const int lgn=__builtin_ctzll(n);
    std::fill(st_1,st_1+(_lg_itth>>1),info->bwbr);
    std::fill(st_1+(_lg_itth>>1),st_1+(_mxlg>>1),info->bwbi);
    const idt nn=n>>(lgn&1),mm=std::min(nn,_itth);
    // I256 rr=_mm256_set1_epi32((info->mod-1)>>(lgn+3));
    for(idt j=0;j<n;j+=mm){
        // const auto pr2=load256(&info->pr2i),pr4=load256(&info->pr4i);
        // const auto pr2Niv=load256(&info->pr2iniv),pr4Niv=load256(&info->pr4iniv);
        // for(idt i=j;i<j+mm;++i){
        //     auto fi=load256(f+i);
        //     const auto rt=rr;
        //     rr=shrk32(mul_bfxd(rr,load256(info->rt4i+__builtin_ctzll(~i)),load256(info->rt4iniv+__builtin_ctzll(~i)),Mod),Mod);
        //     fi=mul_bfxd(Ladd32(neg32_m<0xaa>(fi,Mod2),_mm256_shuffle_epi32(fi,0xb1),Mod2),pr2,pr2Niv,Mod);
        //     fi=mul_bfxd(Ladd32(neg32_m<0xcc>(fi,Mod2),_mm256_shuffle_epi32(fi,0x4e),Mod2),pr4,pr4Niv,Mod);
        //     fi=mul(Ladd32(neg32_m<0xf0>(fi,Mod2),_mm256_permute2x128_si256(fi,fi,1),Mod2),rt,Niv,Mod);
        //     store256(f+i,fi);
        // }
        I256*const g=f+j;
        int t=2,p=0;
        for(idt l=4,L=1;l<=mm;L=l,l<<=2,t+=2,++p){
            auto rt=load256(st_1+p);
            for(idt i=0,k=j>>t;i<mm;i+=l,++k){
                const auto r1=_mm256_permutevar8x32_epi32(rt,id);
                const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB);
                const auto r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
                for(idt j=0;j<L;++j){
                    auto const p0=g+i+j,p1=p0+L,p2=p1+L,p3=p2+L;
                    const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
                    const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
                    const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
                    const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
                    const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
                    const auto u0=shrk32(h0,Mod2),u1=mul_bsm(h1,r1,Niv,Mod);
                    const auto u2=mul_bsm(h2,r2,Niv,Mod),u3=mul_bsm(h3,r3,Niv,Mod);
                    store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
                }
                rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~k)),Mod);
            }
            store256(st_1+p,rt);
        }
        int tt=std::min(__builtin_ctzll(~(j>>_lg_itth))+_lg_itth,lgn);
        for(idt L=_itth,l=L<<2;t<=tt;L=l,l<<=2,t+=2,++p){
            if((j+_itth)==l){
                if(shrk && l==n){
                    for(idt i=0;i<L;++i){
                        auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
                        const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
                        const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
                        const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
                        const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
                        const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
                        const auto u0=shrk32(h0,Mod),u1=shrk32(h1,Mod);
                        const auto u2=shrk32(h2,Mod),u3=shrk32(h3,Mod);
                        store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
                    }
                }
                else{
                    for(idt i=0;i<L;++i){
                        auto const p0=f+i,p1=p0+L,p2=p1+L,p3=p2+L;
                        const auto f2=load256(p2),f3=load256(p3),f0=load256(p0),f1=load256(p1);
                        const auto g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod),g2=add32(f2,f3,Mod2);
                        const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
                        const auto h0=add32(g0,g2,Mod2),h1=add32(g1,g3,Mod2);
                        const auto h2=sub32(g0,g2,Mod2),h3=sub32(g1,g3,Mod2);
                        store256(p0,h0),store256(p1,h1),store256(p2,h2),store256(p3,h3);
                    }
                }
            }
            else{
                auto rt=load256(st_1+p);
                const auto r1=_mm256_permutevar8x32_epi32(rt,id);
                const auto r1Niv=_mm256_permutevar8x32_epi32(_mm256_mul_epu32(rt,Niv),id);
                rt=mul_upd_rt(rt,load256(info->rt3i+__builtin_ctzll(~j>>t)),Mod);
                const auto r2=_mm256_shuffle_epi32(r1,_MM_PERM_BBBB),r3=_mm256_shuffle_epi32(r1,_MM_PERM_DDDD);
                const auto r2Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_BBBB),r3Niv=_mm256_shuffle_epi32(r1Niv,_MM_PERM_DDDD);
                store256(st_1+p,rt);
                for(idt i=0;i<L;++i){
                    auto const p0=f+j+_itth-l+i,p1=p0+L,p2=p1+L,p3=p2+L;
                    const auto f0=load256(p0),f1=load256(p1),f2=load256(p2),f3=load256(p3);
                    const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
                    const auto g2=add32(f2,f3,Mod2),g3=mul_bsmfxd(Lsub32(f3,f2,Mod2),Img,ImgNiv,Mod);
                    const auto h0=Ladd32(g0,g2,Mod2),h1=Ladd32(g1,g3,Mod2);
                    const auto h2=Lsub32(g0,g2,Mod2),h3=Lsub32(g1,g3,Mod2);
                    const auto u0=shrk32(h0,Mod2),u1=mul_bsmfxd(h1,r1,r1Niv,Mod);
                    const auto u2=mul_bsmfxd(h2,r2,r2Niv,Mod),u3=mul_bsmfxd(h3,r3,r3Niv,Mod);
                    store256(p0,u0),store256(p1,u1),store256(p2,u2),store256(p3,u3);
                }
            }
        }
    }
    if(shrk && nn==n && n<=_itth){
        for(idt i=0;i<n;++i){
            const auto f0=load256(f+i);
            store256(f+i,shrk32(f0,Mod));
        }
    }
    if(nn!=n){
        for(idt i=0;i<nn;++i){
            auto const p0=f+i,p1=f+nn+i;
            const auto f0=load256(p0),f1=load256(p1);
            const auto g0=add32(f0,f1,Mod2),g1=sub32(f0,f1,Mod2);
            if constexpr(shrk){
                const auto h0=shrk32(g0,Mod),h1=shrk32(g1,Mod);
                store256(p0,h0),store256(p1,h1);
            }
            else{
                store256(p0,g0),store256(p1,g1);
            }
        }
    }
}
// Returns fx * f[0,8) * g[0,8) (mod x^8 - ww).
[[gnu::always_inline]] inline I256 convolve8(const I256*f,const I256*g,I256 ww,I256 fx,I256 Niv,I256 Mod,I256 Mod2){
    const auto raa=load256(f),rbb=load256(g);
    const auto taa=shrk32(raa,Mod2),bb=shrk32(mul_bsm(rbb,fx,Niv,Mod),Mod);
    const auto aw=shrk32(mul_bsm(taa,ww,Niv,Mod),Mod);
    const auto aa=shrk32(taa,Mod);
    const auto awa=_mm256_permute2x128_si256(aa,aw,3);
    
    const auto b0=_mm256_permute4x64_epi64(bb,0x00),b1=_mm256_shuffle_epi32(b0,_MM_PERM_CDAB);
    const auto a0=aa,a1=_mm256_srli_epi64(a0,32);
    const auto aw7=_mm256_alignr_epi8(aa,awa,12);
    auto res00=_mm256_mul_epu32(a0,b0);
    auto res01=_mm256_mul_epu32(a1,b0);
    auto res10=_mm256_mul_epu32(aw7,b1);
    auto res11=_mm256_mul_epu32(a0,b1);

    const auto b2=_mm256_permute4x64_epi64(bb,0x55),b3=_mm256_shuffle_epi32(b2,_MM_PERM_CDAB);
    const auto aw6=_mm256_alignr_epi8(aa,awa,8);
    const auto aw5=_mm256_alignr_epi8(aa,awa,4);
    res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw6,b2));
    res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw7,b2));
    res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw5,b3));
    res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw6,b3));

    const auto b4=_mm256_permute4x64_epi64(bb,0xaa),b5=_mm256_shuffle_epi32(b4,_MM_PERM_CDAB);
    const auto aw3=_mm256_alignr_epi8(awa,aw,12);
    res00=_mm256_add_epi64(res00,_mm256_mul_epu32(awa,b4));
    res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw5,b4));
    res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw3,b5));
    res11=_mm256_add_epi64(res11,_mm256_mul_epu32(awa,b5));

    const auto b6=_mm256_permute4x64_epi64(bb,0xff),b7=_mm256_shuffle_epi32(b6,_MM_PERM_CDAB);
    const auto aw2=_mm256_alignr_epi8(awa,aw,8);
    const auto aw1=_mm256_alignr_epi8(awa,aw,4);
    res00=_mm256_add_epi64(res00,_mm256_mul_epu32(aw2,b6));
    res01=_mm256_add_epi64(res01,_mm256_mul_epu32(aw3,b6));
    res10=_mm256_add_epi64(res10,_mm256_mul_epu32(aw1,b7));
    res11=_mm256_add_epi64(res11,_mm256_mul_epu32(aw2,b7));

    res00=_mm256_add_epi64(res00,res10);
    res01=_mm256_add_epi64(res01,res11);

    return shrk32(reduce(res00,res01,Niv,Mod),Mod2);
}
inline void vector_convolution_direct(I256*f,const I256*g,idt lm,const FNTT32_info*const info){
    u32 RR=info->one;
    const auto mod=info->mod,niv=info->niv;
    const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
    const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
    for(idt i=0;i<lm;++i){
        store256(f+i,convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2));
        RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
    }
}
inline void vector_convolution_accumulate(I256*const result,const I256*const f,
                                          const I256*const g,idt lm,
                                          const FNTT32_info*const info){
    u32 RR=info->one;
    const auto mod=info->mod,niv=info->niv;
    const auto Fx=_mm256_set1_epi32(mul_s((mod-((mod-1)>>(__builtin_ctzll(lm)))),info->r3,niv,mod));
    const auto Niv=_mm256_set1_epi32(niv),Mod=_mm256_set1_epi32(mod),Mod2=_mm256_set1_epi32(info->mod2);
    for(idt i=0;i<lm;++i){
        const auto product=convolve8(f+i,g+i,_mm256_set1_epi32(RR),Fx,Niv,Mod,Mod2);
        store256(result+i,add32(load256(result+i),product,Mod2));
        RR=mul(RR,info->RT1[__builtin_ctzll(~i)],niv,mod);
    }
}

}  // namespace fast998_v2
}  // namespace internal
}  // namespace fps
}  // namespace m1une

#endif  // M1UNE_FPS_HAS_X86_SIMD


#line 24 "math/fps/convolution.hpp"
#ifdef M1UNE_FPS_HAS_X86_SIMD
#pragma GCC pop_options
#endif

#line 1 "math/modint.hpp"



#line 6 "math/modint.hpp"
#include <iostream>
#line 9 "math/modint.hpp"

namespace m1une {
namespace math {

template <uint32_t Modulus>
struct ModInt {
    static_assert(0 < Modulus, "Modulus must be positive");

   private:
    uint32_t _v;

   public:
    static constexpr uint32_t mod() {
        return Modulus;
    }

    static constexpr ModInt raw(uint32_t v) noexcept {
        ModInt x;
        x._v = v;
        return x;
    }

    constexpr ModInt() noexcept : _v(0) {}

    template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
    constexpr ModInt(Integer v) noexcept {
        if constexpr (std::is_signed_v<Integer>) {
            int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(Modulus);
            if (x < 0) x += Modulus;
            _v = static_cast<uint32_t>(x);
        } else {
            _v = static_cast<uint32_t>(static_cast<uint64_t>(v) % Modulus);
        }
    }

    constexpr uint32_t val() const noexcept {
        return _v;
    }

    constexpr ModInt& operator++() noexcept {
        _v++;
        if (_v == Modulus) _v = 0;
        return *this;
    }

    constexpr ModInt& operator--() noexcept {
        if (_v == 0) _v = Modulus;
        _v--;
        return *this;
    }

    constexpr ModInt operator++(int) noexcept {
        ModInt res = *this;
        ++*this;
        return res;
    }

    constexpr ModInt operator--(int) noexcept {
        ModInt res = *this;
        --*this;
        return res;
    }

    constexpr ModInt& operator+=(const ModInt& rhs) noexcept {
        _v += rhs._v;
        if (_v >= Modulus) _v -= Modulus;
        return *this;
    }

    constexpr ModInt& operator-=(const ModInt& rhs) noexcept {
        _v -= rhs._v;
        if (_v >= Modulus) _v += Modulus;
        return *this;
    }

    constexpr ModInt& operator*=(const ModInt& rhs) noexcept {
        uint64_t z = _v;
        z *= rhs._v;
        _v = static_cast<uint32_t>(z % Modulus);
        return *this;
    }

    constexpr ModInt& operator/=(const ModInt& rhs) noexcept {
        return *this *= rhs.inv();
    }

    constexpr ModInt operator+(const ModInt& rhs) const noexcept {
        return ModInt(*this) += rhs;
    }
    constexpr ModInt operator-(const ModInt& rhs) const noexcept {
        return ModInt(*this) -= rhs;
    }
    constexpr ModInt operator*(const ModInt& rhs) const noexcept {
        return ModInt(*this) *= rhs;
    }
    constexpr ModInt operator/(const ModInt& rhs) const noexcept {
        return ModInt(*this) /= rhs;
    }

    constexpr bool operator==(const ModInt& rhs) const noexcept {
        return _v == rhs._v;
    }
    constexpr bool operator!=(const ModInt& rhs) const noexcept {
        return _v != rhs._v;
    }

    constexpr ModInt pow(long long n) const noexcept {
        ModInt res = raw(1 % Modulus);
        ModInt x = n < 0 ? inv() : *this;
        uint64_t exponent = n < 0 ? uint64_t(-(n + 1)) + 1 : uint64_t(n);
        while (exponent > 0) {
            if (exponent & 1) res *= x;
            x *= x;
            exponent >>= 1;
        }
        return res;
    }

    constexpr ModInt inv() const noexcept {
        int64_t a = _v, b = Modulus, u = 1, v = 0;
        while (b) {
            int64_t t = a / b;
            a -= t * b;
            std::swap(a, b);
            u -= t * v;
            std::swap(u, v);
        }
        assert(a == 1);
        u %= Modulus;
        if (u < 0) u += Modulus;
        return raw(static_cast<uint32_t>(u));
    }

    friend std::ostream& operator<<(std::ostream& os, const ModInt& rhs) {
        return os << rhs._v;
    }

    friend std::istream& operator>>(std::istream& is, ModInt& rhs) {
        long long v;
        is >> v;
        rhs = ModInt(v);
        return is;
    }
};

using modint998244353 = ModInt<998244353>;
using modint1000000007 = ModInt<1000000007>;

template <int Id = 0>
struct DynamicModInt {
   private:
    uint32_t _v;
    inline static uint32_t _mod = 1;

   public:
    static uint32_t mod() noexcept {
        return _mod;
    }

    static void set_mod(uint32_t modulus) noexcept {
        assert(modulus > 0);
        assert(modulus <= uint32_t(1) << 31);
        _mod = modulus;
    }

    static DynamicModInt raw(uint32_t v) noexcept {
        assert(v < _mod);
        DynamicModInt x;
        x._v = v;
        return x;
    }

    DynamicModInt() noexcept : _v(0) {}

    template <class Integer, std::enable_if_t<std::is_integral_v<Integer>, int> = 0>
    DynamicModInt(Integer v) noexcept {
        if constexpr (std::is_signed_v<Integer>) {
            int64_t x = static_cast<int64_t>(v) % static_cast<int64_t>(_mod);
            if (x < 0) x += _mod;
            _v = static_cast<uint32_t>(x);
        } else {
            _v = static_cast<uint32_t>(static_cast<uint64_t>(v) % _mod);
        }
    }

    uint32_t val() const noexcept {
        return _v;
    }

    DynamicModInt& operator++() noexcept {
        _v++;
        if (_v == _mod) _v = 0;
        return *this;
    }

    DynamicModInt& operator--() noexcept {
        if (_v == 0) _v = _mod;
        _v--;
        return *this;
    }

    DynamicModInt operator++(int) noexcept {
        DynamicModInt result = *this;
        ++*this;
        return result;
    }

    DynamicModInt operator--(int) noexcept {
        DynamicModInt result = *this;
        --*this;
        return result;
    }

    DynamicModInt& operator+=(const DynamicModInt& rhs) noexcept {
        _v += rhs._v;
        if (_v >= _mod) _v -= _mod;
        return *this;
    }

    DynamicModInt& operator-=(const DynamicModInt& rhs) noexcept {
        _v -= rhs._v;
        if (_v >= _mod) _v += _mod;
        return *this;
    }

    DynamicModInt& operator*=(const DynamicModInt& rhs) noexcept {
        _v = static_cast<uint32_t>(uint64_t(_v) * rhs._v % _mod);
        return *this;
    }

    DynamicModInt& operator/=(const DynamicModInt& rhs) noexcept {
        return *this *= rhs.inv();
    }

    DynamicModInt operator+(const DynamicModInt& rhs) const noexcept {
        return DynamicModInt(*this) += rhs;
    }

    DynamicModInt operator-(const DynamicModInt& rhs) const noexcept {
        return DynamicModInt(*this) -= rhs;
    }

    DynamicModInt operator*(const DynamicModInt& rhs) const noexcept {
        return DynamicModInt(*this) *= rhs;
    }

    DynamicModInt operator/(const DynamicModInt& rhs) const noexcept {
        return DynamicModInt(*this) /= rhs;
    }

    bool operator==(const DynamicModInt& rhs) const noexcept {
        return _v == rhs._v;
    }

    bool operator!=(const DynamicModInt& rhs) const noexcept {
        return _v != rhs._v;
    }

    DynamicModInt pow(long long exponent) const noexcept {
        DynamicModInt result = raw(1 % _mod);
        DynamicModInt base = exponent < 0 ? inv() : *this;
        uint64_t magnitude =
            exponent < 0 ? uint64_t(-(exponent + 1)) + 1 : uint64_t(exponent);
        while (magnitude > 0) {
            if (magnitude & 1) result *= base;
            base *= base;
            magnitude >>= 1;
        }
        return result;
    }

    DynamicModInt inv() const noexcept {
        int64_t a = _v, b = _mod, u = 1, v = 0;
        while (b) {
            int64_t quotient = a / b;
            a -= quotient * b;
            std::swap(a, b);
            u -= quotient * v;
            std::swap(u, v);
        }
        assert(a == 1);
        u %= _mod;
        if (u < 0) u += _mod;
        return raw(static_cast<uint32_t>(u));
    }

    friend std::ostream& operator<<(std::ostream& os, const DynamicModInt& rhs) {
        return os << rhs._v;
    }

    friend std::istream& operator>>(std::istream& is, DynamicModInt& rhs) {
        long long value;
        is >> value;
        rhs = DynamicModInt(value);
        return is;
    }
};

}  // namespace math
}  // namespace m1une


#line 29 "math/fps/convolution.hpp"

namespace m1une {
namespace fps {

namespace internal {

template <class Mint, class = void>
struct has_static_modulus : std::false_type {};

template <class Mint>
struct has_static_modulus<
    Mint, std::void_t<decltype(std::integral_constant<uint32_t, Mint::mod()>{})>>
    : std::true_type {};

constexpr uint32_t primitive_root_constexpr(uint32_t mod) {
    if (mod == 2) return 1;
    if (mod == 167772161) return 3;
    if (mod == 469762049) return 3;
    if (mod == 754974721) return 11;
    if (mod == 998244353) return 3;
    if (mod == 1224736769) return 3;

    uint32_t divisors[32] = {};
    int count = 0;
    uint32_t x = mod - 1;
    for (uint32_t p = 2; uint64_t(p) * p <= x; p++) {
        if (x % p != 0) continue;
        divisors[count++] = p;
        while (x % p == 0) x /= p;
    }
    if (x > 1) divisors[count++] = x;

    for (uint32_t g = 2;; g++) {
        bool ok = true;
        for (int i = 0; i < count; i++) {
            uint64_t value = 1;
            uint64_t base = g;
            uint32_t exponent = (mod - 1) / divisors[i];
            while (exponent > 0) {
                if (exponent & 1) value = value * base % mod;
                base = base * base % mod;
                exponent >>= 1;
            }
            if (value == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return g;
    }
}

constexpr int two_adic_order(uint32_t x) {
    int result = 0;
    while ((x & 1) == 0) {
        x >>= 1;
        result++;
    }
    return result;
}

template <class Mint>
struct NttRoots {
    static constexpr int max_base = two_adic_order(Mint::mod() - 1);
    std::array<Mint, max_base + 1> root;
    std::array<Mint, max_base + 1> inverse_root;
    std::array<Mint, max_base> rate;
    std::array<Mint, max_base> inverse_rate;
    std::array<Mint, max_base> rate_radix4;
    std::array<Mint, max_base> inverse_rate_radix4;

    NttRoots() {
        constexpr uint32_t primitive_root = primitive_root_constexpr(Mint::mod());
        for (int level = 1; level <= max_base; level++) {
            root[level] = Mint(primitive_root).pow((Mint::mod() - 1) >> level);
            inverse_root[level] = root[level].inv();
        }
        Mint product = 1;
        Mint inverse_product = 1;
        for (int i = 0; i + 1 < max_base; i++) {
            rate[i] = root[i + 2] * product;
            inverse_rate[i] = inverse_root[i + 2] * inverse_product;
            product *= inverse_root[i + 2];
            inverse_product *= root[i + 2];
        }
        product = 1;
        inverse_product = 1;
        for (int i = 0; i + 2 < max_base; i++) {
            rate_radix4[i] = root[i + 3] * product;
            inverse_rate_radix4[i] = inverse_root[i + 3] * inverse_product;
            product *= inverse_root[i + 3];
            inverse_product *= root[i + 3];
        }
    }
};

template <class Mint>
const NttRoots<Mint>& ntt_roots() {
    static const NttRoots<Mint> roots;
    return roots;
}

template <class Mint>
void ntt(std::vector<Mint>& a, bool inverse, bool normalize = true) {
    const int n = int(a.size());
    assert(n > 0 && (n & (n - 1)) == 0);
    assert((Mint::mod() - 1) % uint32_t(n) == 0);

    const auto& roots = ntt_roots<Mint>();
    const int height = two_adic_order(uint32_t(n));
    if (!inverse) {
        int phase = 0;
        while (phase < height) {
            if (height - phase == 1) {
                const int width = 1 << (height - phase - 1);
                Mint twiddle = 1;
                for (int block = 0; block < (1 << phase); block++) {
                    const int offset = block << (height - phase);
                    for (int i = 0; i < width; i++) {
                        const Mint left = a[offset + i];
                        const Mint right = a[offset + i + width] * twiddle;
                        a[offset + i] = left + right;
                        a[offset + i + width] = left - right;
                    }
                    if (block + 1 != (1 << phase))
                        twiddle *= roots.rate[__builtin_ctz(~uint32_t(block))];
                }
                phase++;
                continue;
            }

            const int width = 1 << (height - phase - 2);
            Mint twiddle = 1;
            const Mint imaginary = roots.root[2];
            for (int block = 0; block < (1 << phase); block++) {
                const Mint twiddle2 = twiddle * twiddle;
                const Mint twiddle3 = twiddle2 * twiddle;
                const int offset = block << (height - phase);
                for (int i = 0; i < width; i++) {
                    const uint64_t mod2 = uint64_t(Mint::mod()) * Mint::mod();
                    const uint64_t a0 = a[offset + i].val();
                    const uint64_t a1 = uint64_t(a[offset + i + width].val()) * twiddle.val();
                    const uint64_t a2 =
                        uint64_t(a[offset + i + 2 * width].val()) * twiddle2.val();
                    const uint64_t a3 =
                        uint64_t(a[offset + i + 3 * width].val()) * twiddle3.val();
                    const uint64_t a1na3i =
                        uint64_t(Mint(a1 + mod2 - a3).val()) * imaginary.val();
                    const uint64_t negative_a2 = mod2 - a2;
                    a[offset + i] = Mint(a0 + a2 + a1 + a3);
                    a[offset + i + width] = Mint(a0 + a2 + 2 * mod2 - a1 - a3);
                    a[offset + i + 2 * width] = Mint(a0 + negative_a2 + a1na3i);
                    a[offset + i + 3 * width] = Mint(a0 + negative_a2 + mod2 - a1na3i);
                }
                if (block + 1 != (1 << phase))
                    twiddle *= roots.rate_radix4[__builtin_ctz(~uint32_t(block))];
            }
            phase += 2;
        }
    } else {
        int phase = height;
        while (phase > 0) {
            if (phase == 1) {
                const int width = 1 << (height - phase);
                Mint twiddle = 1;
                for (int block = 0; block < (1 << (phase - 1)); block++) {
                    const int offset = block << (height - phase + 1);
                    for (int i = 0; i < width; i++) {
                        const Mint left = a[offset + i];
                        const Mint right = a[offset + i + width];
                        a[offset + i] = left + right;
                        a[offset + i + width] = (left - right) * twiddle;
                    }
                    if (block + 1 != (1 << (phase - 1)))
                        twiddle *= roots.inverse_rate[__builtin_ctz(~uint32_t(block))];
                }
                phase--;
                continue;
            }

            const int width = 1 << (height - phase);
            Mint twiddle = 1;
            const Mint inverse_imaginary = roots.inverse_root[2];
            for (int block = 0; block < (1 << (phase - 2)); block++) {
                const Mint twiddle2 = twiddle * twiddle;
                const Mint twiddle3 = twiddle2 * twiddle;
                const int offset = block << (height - phase + 2);
                for (int i = 0; i < width; i++) {
                    const uint64_t a0 = a[offset + i].val();
                    const uint64_t a1 = a[offset + i + width].val();
                    const uint64_t a2 = a[offset + i + 2 * width].val();
                    const uint64_t a3 = a[offset + i + 3 * width].val();
                    const uint64_t a2na3i =
                        uint64_t(Mint((Mint::mod() + a2 - a3) * inverse_imaginary.val()).val());
                    a[offset + i] = Mint(a0 + a1 + a2 + a3);
                    a[offset + i + width] =
                        Mint((a0 + Mint::mod() - a1 + a2na3i) * twiddle.val());
                    a[offset + i + 2 * width] = Mint(
                        (a0 + a1 + 2ULL * Mint::mod() - a2 - a3) * twiddle2.val());
                    a[offset + i + 3 * width] = Mint(
                        (a0 + Mint::mod() - a1 + Mint::mod() - a2na3i) * twiddle3.val());
                }
                if (block + 1 != (1 << (phase - 2)))
                    twiddle *= roots.inverse_rate_radix4[__builtin_ctz(~uint32_t(block))];
            }
            phase -= 2;
        }
        if (normalize) {
            const Mint inverse_n = Mint(n).inv();
            for (Mint& value : a) value *= inverse_n;
        }
    }
}

#ifdef M1UNE_FPS_HAS_X86_SIMD

#pragma GCC push_options
#pragma GCC target("avx2,bmi")

template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_simd(const std::vector<Mint>& a,
                                             const std::vector<Mint>& b) {
    const int result_size = int(a.size() + b.size() - 1);
    int n = 1;
    while (n < result_size) n <<= 1;
    const bool squaring = &a == &b;
    auto* transformed_a = static_cast<uint32_t*>(
        ::operator new[](sizeof(uint32_t) * n, std::align_val_t(32)));
    auto* transformed_b = squaring
                              ? transformed_a
                              : static_cast<uint32_t*>(::operator new[](
                                    sizeof(uint32_t) * n, std::align_val_t(32)));
    if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
        static_assert(sizeof(Mint) == sizeof(uint32_t) && std::is_trivially_copyable_v<Mint>);
        std::memcpy(transformed_a, a.data(), sizeof(uint32_t) * a.size());
        if (!squaring)
            std::memcpy(transformed_b, b.data(), sizeof(uint32_t) * b.size());
    } else {
        for (int i = 0; i < int(a.size()); i++) transformed_a[i] = a[i].val();
        if (!squaring)
            for (int i = 0; i < int(b.size()); i++) transformed_b[i] = b[i].val();
    }
    std::memset(transformed_a + a.size(), 0, sizeof(uint32_t) * (n - a.size()));
    if (!squaring)
        std::memset(transformed_b + b.size(), 0, sizeof(uint32_t) * (n - b.size()));

    static constexpr fast998_v2::FNTT32_info transform(998244353);
    const std::size_t vector_size = std::size_t(n) >> 3;
    fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_a), vector_size, &transform);
    if (!squaring)
        fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed_b), vector_size,
                              &transform);
    fast998_v2::vector_convolution_direct(
        reinterpret_cast<__m256i*>(transformed_a),
        reinterpret_cast<const __m256i*>(transformed_b), vector_size, &transform);
    fast998_v2::vector_dit<true>(reinterpret_cast<__m256i*>(transformed_a), vector_size,
                                 &transform);

    std::vector<Mint> result(result_size);
    for (int j = 0; j < result_size; j++) result[j] = Mint::raw(transformed_a[j]);
    ::operator delete[](transformed_a, std::align_val_t(32));
    if (!squaring) ::operator delete[](transformed_b, std::align_val_t(32));
    return result;
}

#pragma GCC pop_options

#endif

}  // namespace internal

template <class Mint>
std::vector<Mint> convolution_naive(const std::vector<Mint>& a, const std::vector<Mint>& b) {
    if (a.empty() || b.empty()) return {};
    std::vector<Mint> result(a.size() + b.size() - 1);
    if (a.size() < b.size()) {
        for (int i = 0; i < int(a.size()); i++) {
            for (int j = 0; j < int(b.size()); j++) result[i + j] += a[i] * b[j];
        }
    } else {
        for (int j = 0; j < int(b.size()); j++) {
            for (int i = 0; i < int(a.size()); i++) result[i + j] += a[i] * b[j];
        }
    }
    return result;
}

template <class Mint>
std::vector<Mint> convolution_ntt(const std::vector<Mint>& a, const std::vector<Mint>& b) {
    const int result_size = int(a.size() + b.size() - 1);
    int n = 1;
    while (n < result_size) n <<= 1;
    assert((Mint::mod() - 1) % uint32_t(n) == 0);

#ifdef M1UNE_FPS_HAS_X86_SIMD
    if constexpr (Mint::mod() == 998244353) {
        if (n >= 64 && __builtin_cpu_supports("avx2"))
            return internal::convolution_998244353_simd(a, b);
    }
#endif

    // Allocate the padded buffers directly.  Constructing from the inputs and
    // then resizing used to allocate and copy both large operands twice.
    const bool squaring = &a == &b;
    std::vector<Mint> fa(n);
    std::copy(a.begin(), a.end(), fa.begin());
    internal::ntt(fa, false);
    const Mint inverse_n = Mint(n).inv();
    if (squaring) {
        for (int i = 0; i < n; i++) fa[i] *= fa[i] * inverse_n;
    } else {
        std::vector<Mint> fb(n);
        std::copy(b.begin(), b.end(), fb.begin());
        internal::ntt(fb, false);
        for (int i = 0; i < n; i++) fa[i] *= fb[i] * inverse_n;
    }
    internal::ntt(fa, true, false);
    fa.resize(result_size);
    return fa;
}

namespace internal {

template <class Mint>
std::vector<Mint> convolution_998244353_blocked_scalar(const std::vector<Mint>& a,
                                                       const std::vector<Mint>& b,
                                                       int transform_size) {
    assert(Mint::mod() == 998244353);
    assert(transform_size >= 2 && (transform_size & (transform_size - 1)) == 0);
    assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);

    const int block_size = transform_size / 2;
    const int a_blocks = int((a.size() + block_size - 1) / block_size);
    const int b_blocks = int((b.size() + block_size - 1) / block_size);

    auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
        std::vector<std::vector<Mint>> blocks;
        blocks.reserve(block_count);
        for (int block = 0; block < block_count; block++) {
            const int begin = block * block_size;
            const int count = std::min(block_size, int(values.size()) - begin);
            std::vector<Mint> transformed(transform_size);
            std::copy_n(values.begin() + begin, count, transformed.begin());
            ntt(transformed, false);
            blocks.emplace_back(std::move(transformed));
        }
        return blocks;
    };

    std::vector<std::vector<Mint>> transformed_a = transform_blocks(a, a_blocks);
    std::vector<std::vector<Mint>> transformed_b = transform_blocks(b, b_blocks);
    const int result_size = int(a.size() + b.size() - 1);
    std::vector<Mint> result(result_size);
    std::vector<Mint> transformed_result(transform_size);
    for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
        std::fill(transformed_result.begin(), transformed_result.end(), Mint(0));
        const int first_a = std::max(0, diagonal - (b_blocks - 1));
        const int last_a = std::min(a_blocks - 1, diagonal);
        for (int a_block = first_a; a_block <= last_a; a_block++) {
            const int b_block = diagonal - a_block;
            for (int i = 0; i < transform_size; i++)
                transformed_result[i] +=
                    transformed_a[a_block][i] * transformed_b[b_block][i];
        }
        ntt(transformed_result, true);

        const int output_offset = diagonal * block_size;
        const int output_count = std::min(transform_size, result_size - output_offset);
        for (int i = 0; i < output_count; i++)
            result[output_offset + i] += transformed_result[i];
    }
    return result;
}

#ifdef M1UNE_FPS_HAS_X86_SIMD

class AlignedUint32Buffer {
   private:
    uint32_t* data_;

   public:
    explicit AlignedUint32Buffer(std::size_t size)
        : data_(static_cast<uint32_t*>(
              ::operator new[](sizeof(uint32_t) * size, std::align_val_t(32)))) {}

    AlignedUint32Buffer(const AlignedUint32Buffer&) = delete;
    AlignedUint32Buffer& operator=(const AlignedUint32Buffer&) = delete;

    AlignedUint32Buffer(AlignedUint32Buffer&& other) noexcept : data_(other.data_) {
        other.data_ = nullptr;
    }

    AlignedUint32Buffer& operator=(AlignedUint32Buffer&& other) noexcept {
        if (this == &other) return *this;
        ::operator delete[](data_, std::align_val_t(32));
        data_ = other.data_;
        other.data_ = nullptr;
        return *this;
    }

    ~AlignedUint32Buffer() {
        ::operator delete[](data_, std::align_val_t(32));
    }

    uint32_t* data() {
        return data_;
    }

    const uint32_t* data() const {
        return data_;
    }
};

template <class Mint>
__attribute__((target("avx2,bmi"), hot))
std::vector<Mint> convolution_998244353_blocked_simd(const std::vector<Mint>& a,
                                                     const std::vector<Mint>& b,
                                                     int transform_size) {
    assert(Mint::mod() == 998244353);
    assert(transform_size >= 64 && (transform_size & (transform_size - 1)) == 0);
    assert((Mint::mod() - 1) % uint32_t(transform_size) == 0);

    const int block_size = transform_size / 2;
    const int a_blocks = int((a.size() + block_size - 1) / block_size);
    const int b_blocks = int((b.size() + block_size - 1) / block_size);
    static constexpr fast998_v2::FNTT32_info transform(998244353);
    const std::size_t vector_size = std::size_t(transform_size) / 8;

    auto transform_blocks = [&](const std::vector<Mint>& values, int block_count) {
        std::vector<AlignedUint32Buffer> blocks;
        blocks.reserve(block_count);
        for (int block = 0; block < block_count; block++) {
            const int begin = block * block_size;
            const int count = std::min(block_size, int(values.size()) - begin);
            AlignedUint32Buffer transformed(transform_size);
            if constexpr (std::is_same_v<Mint, math::ModInt<998244353>>) {
                static_assert(sizeof(Mint) == sizeof(uint32_t) &&
                              std::is_trivially_copyable_v<Mint>);
                std::memcpy(transformed.data(), values.data() + begin,
                            sizeof(uint32_t) * count);
            } else {
                for (int i = 0; i < count; i++)
                    transformed.data()[i] = values[begin + i].val();
            }
            std::memset(transformed.data() + count, 0,
                        sizeof(uint32_t) * (transform_size - count));
            fast998_v2::vector_dif(reinterpret_cast<__m256i*>(transformed.data()),
                                   vector_size, &transform);
            blocks.emplace_back(std::move(transformed));
        }
        return blocks;
    };

    std::vector<AlignedUint32Buffer> transformed_a = transform_blocks(a, a_blocks);
    std::vector<AlignedUint32Buffer> transformed_b = transform_blocks(b, b_blocks);
    const int result_size = int(a.size() + b.size() - 1);
    std::vector<Mint> result(result_size);
    AlignedUint32Buffer transformed_result(transform_size);
    for (int diagonal = 0; diagonal < a_blocks + b_blocks - 1; diagonal++) {
        std::memset(transformed_result.data(), 0, sizeof(uint32_t) * transform_size);
        const int first_a = std::max(0, diagonal - (b_blocks - 1));
        const int last_a = std::min(a_blocks - 1, diagonal);
        for (int a_block = first_a; a_block <= last_a; a_block++) {
            const int b_block = diagonal - a_block;
            fast998_v2::vector_convolution_accumulate(
                reinterpret_cast<__m256i*>(transformed_result.data()),
                reinterpret_cast<const __m256i*>(transformed_a[a_block].data()),
                reinterpret_cast<const __m256i*>(transformed_b[b_block].data()),
                vector_size, &transform);
        }
        fast998_v2::vector_dit<true>(
            reinterpret_cast<__m256i*>(transformed_result.data()), vector_size,
            &transform);

        const int output_offset = diagonal * block_size;
        const int output_count = std::min(transform_size, result_size - output_offset);
        for (int i = 0; i < output_count; i++) {
            uint32_t value = result[output_offset + i].val() + transformed_result.data()[i];
            if (value >= Mint::mod()) value -= Mint::mod();
            result[output_offset + i] = Mint::raw(value);
        }
    }
    return result;
}

#endif

template <class Mint>
std::vector<Mint> convolution_998244353_blocked(const std::vector<Mint>& a,
                                                const std::vector<Mint>& b,
                                                int transform_size = 1 << 23) {
#ifdef M1UNE_FPS_HAS_X86_SIMD
    if (transform_size >= 64 && __builtin_cpu_supports("avx2"))
        return convolution_998244353_blocked_simd(a, b, transform_size);
#endif
    return convolution_998244353_blocked_scalar(a, b, transform_size);
}

}  // namespace internal

template <class Mint>
std::vector<Mint> convolution(const std::vector<Mint>& a, const std::vector<Mint>& b) {
    if (a.empty() || b.empty()) return {};
    if (std::min(a.size(), b.size()) <= 32) return convolution_naive(a, b);

    const int result_size = int(a.size() + b.size() - 1);
    int n = 1;
    while (n < result_size) n <<= 1;
    if constexpr (internal::has_static_modulus<Mint>::value) {
        if constexpr (Mint::mod() == 998244353) {
            if (n > (1 << 23))
                return internal::convolution_998244353_blocked(a, b);
        }
        if ((Mint::mod() - 1) % uint32_t(n) == 0) return convolution_ntt(a, b);
    }

    using Mint1 = math::ModInt<167772161>;
    using Mint2 = math::ModInt<469762049>;
    using Mint3 = math::ModInt<754974721>;
    assert(n <= (1 << 24));

    [[maybe_unused]] const unsigned __int128 coefficient_bound =
        static_cast<unsigned __int128>(std::min(a.size(), b.size())) * (Mint::mod() - 1) *
        (Mint::mod() - 1);
    [[maybe_unused]] const unsigned __int128 crt_modulus =
        static_cast<unsigned __int128>(Mint1::mod()) * Mint2::mod() * Mint3::mod();
    assert(coefficient_bound < crt_modulus);

    auto converted_convolution = [&]<class OtherMint>() {
        std::vector<OtherMint> converted_a(a.size());
        std::vector<OtherMint> converted_b(b.size());
        for (int i = 0; i < int(a.size()); i++) converted_a[i] = OtherMint(a[i].val());
        for (int i = 0; i < int(b.size()); i++) converted_b[i] = OtherMint(b[i].val());
        return convolution_ntt(converted_a, converted_b);
    };
    std::vector<Mint1> c1 = converted_convolution.template operator()<Mint1>();
    std::vector<Mint2> c2 = converted_convolution.template operator()<Mint2>();
    std::vector<Mint3> c3 = converted_convolution.template operator()<Mint3>();
    static const uint64_t inverse_mod1_mod2 = Mint2(Mint1::mod()).inv().val();
    static const uint64_t mod1_mod3 = Mint1::mod() % Mint3::mod();
    static const uint64_t mod1_mod2_mod3 =
        mod1_mod3 * (Mint2::mod() % Mint3::mod()) % Mint3::mod();
    static const uint64_t inverse_mod1_mod2_mod3 = Mint3(uint32_t(mod1_mod2_mod3)).inv().val();

    const uint64_t target_mod = Mint::mod();
    const uint64_t mod1_target = Mint1::mod() % target_mod;
    const uint64_t mod1_mod2_target = mod1_target * (Mint2::mod() % target_mod) % target_mod;
    std::vector<Mint> result(result_size);
    for (int i = 0; i < result_size; i++) {
        const uint64_t r1 = c1[i].val();
        const uint64_t r2 = c2[i].val();
        const uint64_t r3 = c3[i].val();
        const uint64_t first =
            (r2 + Mint2::mod() - r1 % Mint2::mod()) % Mint2::mod() * inverse_mod1_mod2 %
            Mint2::mod();
        const uint64_t combined_mod3 =
            (r1 % Mint3::mod() + mod1_mod3 * (first % Mint3::mod())) % Mint3::mod();
        const uint64_t second =
            (r3 + Mint3::mod() - combined_mod3) % Mint3::mod() * inverse_mod1_mod2_mod3 %
            Mint3::mod();

        uint64_t value = r1 % target_mod;
        value = (value + mod1_target * (first % target_mod)) % target_mod;
        value = (value + mod1_mod2_target * (second % target_mod)) % target_mod;
        result[i] = Mint::raw(uint32_t(value));
    }
    return result;
}

}  // namespace fps
}  // namespace m1une

#ifdef M1UNE_FPS_HAS_X86_SIMD
#undef M1UNE_FPS_HAS_X86_SIMD
#endif


#line 14 "graph/tree/distance_frequency.hpp"

namespace m1une {
namespace tree {

namespace distance_frequency_detail {

template <class Mint, class T>
std::vector<Mint> count_ordered_pairs(
    const m1une::graph::Graph<T>& tree,
    const CentroidDecomposition<T>& decomposition
) {
    const int size = tree.size();
    std::vector<Mint> count(static_cast<std::size_t>(size));
    std::vector<char> removed(std::size_t(size), false);
    std::vector<Mint> histogram;
    std::vector<std::pair<int, int>> stack;
    std::vector<int> parent(std::size_t(size), -1);

    for (int centroid : decomposition.order) {
        std::vector<Mint> total(1, Mint(1));
        for (const auto& edge : tree[centroid]) {
            if (!edge.alive || removed[std::size_t(edge.to)]) continue;

            histogram.clear();
            stack.clear();
            stack.emplace_back(edge.to, 1);
            parent[std::size_t(edge.to)] = centroid;
            while (!stack.empty()) {
                const auto [vertex, distance] = stack.back();
                stack.pop_back();
                if (int(histogram.size()) <= distance) {
                    histogram.resize(std::size_t(distance + 1));
                }
                histogram[std::size_t(distance)] += Mint(1);

                for (const auto& next : tree[vertex]) {
                    if (!next.alive || removed[std::size_t(next.to)]) continue;
                    if (next.to == parent[std::size_t(vertex)]) continue;
                    parent[std::size_t(next.to)] = vertex;
                    stack.emplace_back(next.to, distance + 1);
                }
            }

            if (total.size() < histogram.size()) {
                total.resize(histogram.size());
            }
            for (std::size_t distance = 0; distance < histogram.size(); distance++) {
                total[distance] += histogram[distance];
            }

            const std::vector<Mint> within_component =
                m1une::fps::convolution(histogram, histogram);
            const std::size_t limit = std::min(count.size(), within_component.size());
            for (std::size_t distance = 0; distance < limit; distance++) {
                count[distance] -= within_component[distance];
            }
        }

        const std::vector<Mint> through_centroid =
            m1une::fps::convolution(total, total);
        const std::size_t limit = std::min(count.size(), through_centroid.size());
        for (std::size_t distance = 0; distance < limit; distance++) {
            count[distance] += through_centroid[distance];
        }
        removed[std::size_t(centroid)] = true;
    }
    return count;
}

inline std::uint64_t combine_residues(std::uint32_t first, std::uint32_t second) {
    using First = m1une::math::ModInt<998244353>;
    using Second = m1une::math::ModInt<924844033>;
    static const std::uint64_t inverse = Second(First::mod()).inv().val();
    const std::uint64_t offset =
        (std::uint64_t(second) + Second::mod() - first % Second::mod()) %
        Second::mod();
    const std::uint64_t multiplier = offset * inverse % Second::mod();
    return std::uint64_t(first) + std::uint64_t(First::mod()) * multiplier;
}

}  // namespace distance_frequency_detail

template <class T>
std::vector<long long> tree_distance_frequency(
    const m1une::graph::Graph<T>& tree
) {
    const int size = tree.size();
    assert(tree.edge_count() == std::max(0, size - 1));
    if (size == 0) return {};

    const CentroidDecomposition<T> decomposition(tree);
    assert(decomposition.roots.size() == 1);

    using First = m1une::math::ModInt<998244353>;
    using Second = m1une::math::ModInt<924844033>;
    assert(
        std::uint64_t(size) * std::uint64_t(size - 1) <
        std::uint64_t(First::mod()) * Second::mod()
    );
    const std::vector<First> first =
        distance_frequency_detail::count_ordered_pairs<First>(
            tree,
            decomposition
        );
    const std::vector<Second> second =
        distance_frequency_detail::count_ordered_pairs<Second>(
            tree,
            decomposition
        );

    std::vector<long long> result(static_cast<std::size_t>(size));
    result[0] = size;
    for (int distance = 1; distance < size; distance++) {
        const std::uint64_t ordered =
            distance_frequency_detail::combine_residues(
                first[std::size_t(distance)].val(),
                second[std::size_t(distance)].val()
            );
        assert((ordered & 1) == 0);
        result[std::size_t(distance)] = static_cast<long long>(ordered / 2);
    }
    return result;
}

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/dsu_on_tree.hpp"



#line 7 "graph/tree/dsu_on_tree.hpp"

#line 9 "graph/tree/dsu_on_tree.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct DsuOnTree {
    int n;
    int root;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<int> subtree_size;
    std::vector<int> heavy_child;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> order;
    std::vector<std::vector<int>> children;

    DsuOnTree() : n(0), root(-1) {}

    explicit DsuOnTree(
        const m1une::graph::Graph<T>& graph,
        int root_vertex = 0
    ) {
        build(graph, root_vertex);
    }

    void build(
        const m1une::graph::Graph<T>& graph,
        int root_vertex = 0
    ) {
        n = graph.size();
        root = n == 0 ? -1 : root_vertex;
        parent.assign(n, -2);
        parent_edge.assign(n, -1);
        depth.assign(n, 0);
        subtree_size.assign(n, 1);
        heavy_child.assign(n, -1);
        tin.assign(n, -1);
        tout.assign(n, -1);
        order.clear();
        order.reserve(n);
        children.assign(n, {});
        if (n == 0) return;

        assert(0 <= root && root < n);
        std::vector<int> stack;
        stack.push_back(root);
        parent[root] = -1;
        while (!stack.empty()) {
            int vertex = stack.back();
            stack.pop_back();
            tin[vertex] = int(order.size());
            order.push_back(vertex);

            for (const auto& edge : graph[vertex]) {
                if (!edge.alive || parent[edge.to] != -2) continue;
                parent[edge.to] = vertex;
                parent_edge[edge.to] = edge.id;
                depth[edge.to] = depth[vertex] + 1;
                children[vertex].push_back(edge.to);
                stack.push_back(edge.to);
            }
        }
        assert(int(order.size()) == n);

        for (int index = n - 1; index >= 0; --index) {
            int vertex = order[index];
            for (int child : children[vertex]) {
                subtree_size[vertex] += subtree_size[child];
                if (
                    heavy_child[vertex] == -1 ||
                    subtree_size[heavy_child[vertex]] < subtree_size[child]
                ) {
                    heavy_child[vertex] = child;
                }
            }
            tout[vertex] = tin[vertex] + subtree_size[vertex];
        }
    }

    int size() const {
        return n;
    }

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

    std::pair<int, int> subtree_range(int vertex) const {
        assert(0 <= vertex && vertex < n);
        return {tin[vertex], tout[vertex]};
    }

    // Runs DSU on tree. `add(v)` inserts one vertex into the maintained state,
    // `remove(v)` erases it, and `answer(v)` observes the state for subtree(v).
    template <class Add, class Remove, class Answer>
    void run(Add add, Remove remove, Answer answer) const {
        if (n == 0) return;

        enum ActionType {
            Process,
            AddSubtree,
            AddVertex,
            AnswerVertex,
            RemoveSubtree,
        };
        struct Action {
            ActionType type;
            int vertex;
            bool keep;
        };

        std::vector<Action> actions;
        actions.reserve(3 * std::size_t(n));
        actions.push_back(Action{Process, root, true});

        while (!actions.empty()) {
            Action action = actions.back();
            actions.pop_back();
            int vertex = action.vertex;

            if (action.type == AddSubtree) {
                for (int index = tin[vertex]; index < tout[vertex]; ++index) {
                    add(order[index]);
                }
            } else if (action.type == AddVertex) {
                add(vertex);
            } else if (action.type == AnswerVertex) {
                answer(vertex);
            } else if (action.type == RemoveSubtree) {
                for (int index = tin[vertex]; index < tout[vertex]; ++index) {
                    remove(order[index]);
                }
            } else {
                if (!action.keep) {
                    actions.push_back(Action{
                        RemoveSubtree,
                        vertex,
                        false,
                    });
                }
                actions.push_back(Action{AnswerVertex, vertex, false});
                actions.push_back(Action{AddVertex, vertex, false});

                for (int child : children[vertex]) {
                    if (child != heavy_child[vertex]) {
                        actions.push_back(Action{
                            AddSubtree,
                            child,
                            false,
                        });
                    }
                }
                if (heavy_child[vertex] != -1) {
                    actions.push_back(Action{
                        Process,
                        heavy_child[vertex],
                        true,
                    });
                }
                for (int child : children[vertex]) {
                    if (child != heavy_child[vertex]) {
                        actions.push_back(Action{Process, child, false});
                    }
                }
            }
        }
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/euler_tour.hpp"



#line 8 "graph/tree/euler_tour.hpp"

#line 10 "graph/tree/euler_tour.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct EulerTour {
    using cost_type = T;
    using edge_type = m1une::graph::Edge<T>;

    int root;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<T> dist;
    std::vector<int> subtree_size;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> order;
    std::vector<std::vector<int>> children;

   private:
    int _n;

    void check_vertex(int v) const {
        assert(0 <= v && v < _n);
        assert(tin[v] != -1);
    }

   public:
    EulerTour() : root(-1), _n(0) {}
    explicit EulerTour(const m1une::graph::Graph<T>& g, int root_ = 0) {
        build(g, root_);
    }

    void build(const m1une::graph::Graph<T>& g, int root_ = 0) {
        _n = g.size();
        root = _n == 0 ? -1 : root_;
        parent.assign(_n, -2);
        parent_edge.assign(_n, -1);
        depth.assign(_n, 0);
        dist.assign(_n, T(0));
        subtree_size.assign(_n, 0);
        tin.assign(_n, -1);
        tout.assign(_n, -1);
        order.clear();
        order.reserve(_n);
        children.assign(_n, {});

        if (_n == 0) return;
        assert(0 <= root && root < _n);

        struct Frame {
            int v;
            int state;
        };

        std::vector<Frame> stack;
        stack.push_back({root, 0});
        parent[root] = -1;

        while (!stack.empty()) {
            Frame frame = stack.back();
            stack.pop_back();
            int v = frame.v;
            if (frame.state == 0) {
                tin[v] = int(order.size());
                order.push_back(v);
                stack.push_back({v, 1});
                const auto& adj = g[v];
                for (int i = int(adj.size()) - 1; i >= 0; --i) {
                    const auto& e = adj[i];
                    if (!e.alive) continue;
                    if (parent[e.to] != -2) continue;
                    parent[e.to] = v;
                    parent_edge[e.to] = e.id;
                    depth[e.to] = depth[v] + 1;
                    dist[e.to] = dist[v] + e.cost;
                    children[v].push_back(e.to);
                    stack.push_back({e.to, 0});
                }
                std::reverse(children[v].begin(), children[v].end());
            } else {
                subtree_size[v] = 1;
                for (int child : children[v]) subtree_size[v] += subtree_size[child];
                tout[v] = int(order.size());
            }
        }
    }

    int size() const {
        return _n;
    }

    int visited_size() const {
        return int(order.size());
    }

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

    bool is_ancestor(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        return tin[u] <= tin[v] && tout[v] <= tout[u];
    }

    bool in_subtree(int v, int u) const {
        return is_ancestor(u, v);
    }

    std::pair<int, int> subtree_range(int v, bool edge = false) const {
        check_vertex(v);
        return {tin[v] + (edge ? 1 : 0), tout[v]};
    }

    std::vector<int> subtree_vertices(int v) const {
        check_vertex(v);
        return std::vector<int>(order.begin() + tin[v], order.begin() + tout[v]);
    }

    template <class F>
    void for_each_subtree(int v, F f) const {
        auto [l, r] = subtree_range(v);
        for (int i = l; i < r; ++i) f(order[i]);
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/heavy_light_decomposition.hpp"



#line 8 "graph/tree/heavy_light_decomposition.hpp"

#line 10 "graph/tree/heavy_light_decomposition.hpp"

namespace m1une {
namespace tree {

struct HldPathSegment {
    int l;
    int r;
    bool reversed;
};

template <class T = int>
struct HeavyLightDecomposition {
    using cost_type = T;
    using edge_type = m1une::graph::Edge<T>;

    int root;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<T> dist;
    std::vector<int> subtree_size;
    std::vector<int> heavy;
    std::vector<int> head;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> order;

   private:
    int _n;

    void check_vertex(int v) const {
        assert(0 <= v && v < _n);
        assert(tin[v] != -1);
    }

    static void add_segment(std::vector<HldPathSegment>& result, int l, int r, bool reversed) {
        if (l < r) result.push_back({l, r, reversed});
    }

   public:
    HeavyLightDecomposition() : root(-1), _n(0) {}
    explicit HeavyLightDecomposition(const m1une::graph::Graph<T>& g, int root_ = 0) {
        build(g, root_);
    }

    void build(const m1une::graph::Graph<T>& g, int root_ = 0) {
        _n = g.size();
        root = _n == 0 ? -1 : root_;
        parent.assign(_n, -2);
        parent_edge.assign(_n, -1);
        depth.assign(_n, 0);
        dist.assign(_n, T(0));
        subtree_size.assign(_n, 1);
        heavy.assign(_n, -1);
        head.assign(_n, -1);
        tin.assign(_n, -1);
        tout.assign(_n, -1);
        order.clear();
        order.reserve(_n);
        if (_n == 0) return;
        assert(0 <= root && root < _n);

        std::vector<int> dfs_order;
        dfs_order.reserve(_n);
        std::vector<int> stack = {root};
        parent[root] = -1;
        while (!stack.empty()) {
            int v = stack.back();
            stack.pop_back();
            dfs_order.push_back(v);
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                if (parent[e.to] != -2) continue;
                parent[e.to] = v;
                parent_edge[e.to] = e.id;
                depth[e.to] = depth[v] + 1;
                dist[e.to] = dist[v] + e.cost;
                stack.push_back(e.to);
            }
        }

        for (int i = int(dfs_order.size()) - 1; i >= 0; i--) {
            int v = dfs_order[i];
            if (parent[v] == -1) continue;
            int p = parent[v];
            subtree_size[p] += subtree_size[v];
            if (heavy[p] == -1 || subtree_size[heavy[p]] < subtree_size[v]) heavy[p] = v;
        }

        order.assign(dfs_order.size(), -1);
        int timer = 0;
        std::vector<std::pair<int, int>> starts = {std::pair<int, int>{root, root}};
        while (!starts.empty()) {
            auto [start, h] = starts.back();
            starts.pop_back();
            for (int v = start; v != -1; v = heavy[v]) {
                head[v] = h;
                tin[v] = timer;
                order[timer++] = v;
                for (auto it = g[v].rbegin(); it != g[v].rend(); ++it) {
                    if (!it->alive) continue;
                    int to = it->to;
                    if (parent[to] != v || to == heavy[v]) continue;
                    starts.push_back({to, to});
                }
            }
        }
        for (int i = int(dfs_order.size()) - 1; i >= 0; i--) {
            int v = dfs_order[i];
            tout[v] = tin[v] + subtree_size[v];
        }
    }

    int size() const {
        return _n;
    }

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

    bool is_ancestor(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        return tin[u] <= tin[v] && tout[v] <= tout[u];
    }

    int lca(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        while (head[u] != head[v]) {
            if (depth[head[u]] < depth[head[v]]) std::swap(u, v);
            u = parent[head[u]];
        }
        return depth[u] < depth[v] ? u : v;
    }

    int dist_edges(int u, int v) const {
        int w = lca(u, v);
        return depth[u] + depth[v] - 2 * depth[w];
    }

    T dist_cost(int u, int v) const {
        int w = lca(u, v);
        return dist[u] + dist[v] - dist[w] - dist[w];
    }

    int kth_ancestor(int v, int k) const {
        check_vertex(v);
        assert(0 <= k);
        while (v != -1) {
            int h = head[v];
            int len = depth[v] - depth[h];
            if (k <= len) return order[tin[v] - k];
            k -= len + 1;
            v = parent[h];
        }
        return -1;
    }

    int jump(int from, int to, int k) const {
        check_vertex(from);
        check_vertex(to);
        assert(0 <= k);
        int w = lca(from, to);
        int up_len = depth[from] - depth[w];
        int down_len = depth[to] - depth[w];
        if (up_len + down_len < k) return -1;
        if (k <= up_len) return kth_ancestor(from, k);
        return kth_ancestor(to, down_len - (k - up_len));
    }

    std::pair<int, int> subtree_range(int v, bool edge = false) const {
        check_vertex(v);
        return {tin[v] + (edge ? 1 : 0), tout[v]};
    }

    std::vector<HldPathSegment> path_segments(int u, int v, bool edge = false) const {
        check_vertex(u);
        check_vertex(v);
        std::vector<HldPathSegment> result, down;
        while (head[u] != head[v]) {
            if (depth[head[u]] >= depth[head[v]]) {
                add_segment(result, tin[head[u]], tin[u] + 1, true);
                u = parent[head[u]];
            } else {
                add_segment(down, tin[head[v]], tin[v] + 1, false);
                v = parent[head[v]];
            }
        }

        if (depth[u] >= depth[v]) {
            add_segment(result, tin[v] + (edge ? 1 : 0), tin[u] + 1, true);
        } else {
            add_segment(down, tin[u] + (edge ? 1 : 0), tin[v] + 1, false);
        }
        std::reverse(down.begin(), down.end());
        result.insert(result.end(), down.begin(), down.end());
        return result;
    }

    template <class F>
    void for_each_path(int u, int v, F f, bool edge = false) const {
        for (auto seg : path_segments(u, v, edge)) f(seg.l, seg.r, seg.reversed);
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/mo_on_tree.hpp"



#line 7 "graph/tree/mo_on_tree.hpp"

#line 1 "algo/offline/mo.hpp"



#line 6 "algo/offline/mo.hpp"
#include <cmath>
#include <numeric>
#line 9 "algo/offline/mo.hpp"

namespace m1une {
namespace algo {

// Offline Mo's algorithm for half-open array ranges.
struct Mo {
    struct Query {
        int left;
        int right;
        int id;
    };

   private:
    int _n;
    std::vector<Query> _queries;

   public:
    Mo() : _n(0) {}

    explicit Mo(int n) : _n(n) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    int query_count() const {
        return int(_queries.size());
    }

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

    const std::vector<Query>& queries() const {
        return _queries;
    }

    void reserve(int query_capacity) {
        assert(0 <= query_capacity);
        _queries.reserve(query_capacity);
    }

    void clear() {
        _queries.clear();
    }

    // Adds [left, right) and returns its insertion-order ID.
    int add_query(int left, int right) {
        assert(0 <= left && left <= right && right <= _n);
        int id = query_count();
        _queries.push_back(Query{left, right, id});
        return id;
    }

    // Returns query IDs in Mo order. A non-positive block size selects one
    // automatically.
    std::vector<int> order(int block_size = 0) const {
        int query_size = query_count();
        std::vector<int> result(query_size);
        std::iota(result.begin(), result.end(), 0);
        if (query_size == 0) return result;

        if (block_size <= 0) {
            block_size = std::max(1, int(_n / std::sqrt(static_cast<double>(query_size))));
        }

        std::sort(result.begin(), result.end(), [&](int first, int second) {
            const Query& a = _queries[first];
            const Query& b = _queries[second];
            int first_block = a.left / block_size;
            int second_block = b.left / block_size;
            if (first_block != second_block) {
                return first_block < second_block;
            }
            if (first_block & 1) return a.right > b.right;
            return a.right < b.right;
        });
        return result;
    }

    // Maintains [left, right). Each movement callback receives the array index
    // being inserted or erased. `answer(query_id)` stores or reports a result.
    template <class AddLeft, class AddRight, class RemoveLeft, class RemoveRight, class Answer>
    void run(AddLeft add_left, AddRight add_right, RemoveLeft remove_left, RemoveRight remove_right, Answer answer,
             int block_size = 0) const {
        int left = 0;
        int right = 0;
        for (int query_index : order(block_size)) {
            const Query& query = _queries[query_index];
            while (query.left < left) add_left(--left);
            while (right < query.right) add_right(right++);
            while (left < query.left) remove_left(left++);
            while (query.right < right) remove_right(--right);
            answer(query.id);
        }
    }

    // Convenience overload for statistics whose update is independent of
    // which side moves.
    template <class Add, class Remove, class Answer>
    void run(Add add, Remove remove, Answer answer, int block_size = 0) const {
        run(add, add, remove, remove, answer, block_size);
    }
};

}  // namespace algo
}  // namespace m1une


#line 11 "graph/tree/mo_on_tree.hpp"

namespace m1une {
namespace tree {

// Offline Mo's algorithm for static paths in a tree.
template <class T = int>
struct MoOnTree {
    struct Query {
        int from;
        int to;
        int left;
        int right;
        int extra;
        int id;
        bool edge;
    };

    int root;
    std::vector<int> entry;
    std::vector<int> exit;
    std::vector<int> tour;

   private:
    int _n;
    HeavyLightDecomposition<T> _hld;
    m1une::algo::Mo _mo;
    std::vector<Query> _queries;

    void check_vertex(int vertex) const {
        assert(0 <= vertex && vertex < _n);
        assert(entry[vertex] != -1);
    }

    int add_path_query(int from, int to, bool edge) {
        check_vertex(from);
        check_vertex(to);
        assert(_queries.empty() || _queries.front().edge == edge);
        int original_from = from;
        int original_to = to;
        if (entry[from] > entry[to]) std::swap(from, to);

        int ancestor = _hld.lca(from, to);
        int left;
        int right = entry[to] + 1;
        int extra = -1;
        if (ancestor == from) {
            left = entry[from] + int(edge);
        } else {
            left = exit[from];
            if (!edge) extra = ancestor;
        }

        int id = _mo.add_query(left, right);
        _queries.push_back(Query{
            original_from,
            original_to,
            left,
            right,
            extra,
            id,
            edge,
        });
        return id;
    }

   public:
    MoOnTree() : root(-1), _n(0), _mo(0) {}

    explicit MoOnTree(
        const m1une::graph::Graph<T>& graph,
        int root_vertex = 0
    ) : root(-1), _n(0), _mo(0) {
        build(graph, root_vertex);
    }

    void build(
        const m1une::graph::Graph<T>& graph,
        int root_vertex = 0
    ) {
        _n = graph.size();
        root = _n == 0 ? -1 : root_vertex;
        entry.assign(_n, -1);
        exit.assign(_n, -1);
        tour.clear();
        tour.reserve(2 * _n);
        _queries.clear();
        _mo = m1une::algo::Mo(2 * _n);
        _hld.build(graph, root_vertex);
        if (_n == 0) return;

        assert(0 <= root && root < _n);
        for (int vertex = 0; vertex < _n; ++vertex) {
            assert(_hld.parent[vertex] != -2);
        }

        std::vector<std::vector<int>> children(_n);
        for (int vertex = 0; vertex < _n; ++vertex) {
            int parent = _hld.parent[vertex];
            if (parent != -1) children[parent].push_back(vertex);
        }

        struct Event {
            int vertex;
            bool leaving;
        };
        std::vector<Event> stack;
        stack.reserve(2 * _n);
        stack.push_back(Event{root, false});
        while (!stack.empty()) {
            Event event = stack.back();
            stack.pop_back();
            int vertex = event.vertex;
            if (event.leaving) {
                exit[vertex] = int(tour.size());
                tour.push_back(vertex);
                continue;
            }

            entry[vertex] = int(tour.size());
            tour.push_back(vertex);
            stack.push_back(Event{vertex, true});
            const auto& child_list = children[vertex];
            for (int index = int(child_list.size()) - 1; index >= 0; --index) {
                stack.push_back(Event{child_list[index], false});
            }
        }
        assert(int(tour.size()) == 2 * _n);
    }

    int size() const {
        return _n;
    }

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

    int query_count() const {
        return int(_queries.size());
    }

    const std::vector<Query>& queries() const {
        return _queries;
    }

    int parent(int vertex) const {
        check_vertex(vertex);
        return _hld.parent[vertex];
    }

    int parent_edge(int vertex) const {
        check_vertex(vertex);
        return _hld.parent_edge[vertex];
    }

    int depth(int vertex) const {
        check_vertex(vertex);
        return _hld.depth[vertex];
    }

    int lca(int first, int second) const {
        check_vertex(first);
        check_vertex(second);
        return _hld.lca(first, second);
    }

    void reserve(int query_capacity) {
        assert(0 <= query_capacity);
        _queries.reserve(query_capacity);
        _mo.reserve(query_capacity);
    }

    void clear() {
        _queries.clear();
        _mo.clear();
    }

    // Adds an inclusive vertex-path query and returns its insertion-order ID.
    // Vertex and edge queries cannot be mixed in one collection.
    int add_query(int from, int to) {
        return add_path_query(from, to, false);
    }

    // Adds an edge-path query. Each edge is represented by its child vertex.
    int add_edge_query(int from, int to) {
        return add_path_query(from, to, true);
    }

    std::vector<int> order(int block_size = 0) const {
        return _mo.order(block_size);
    }

    // `add(v)` and `remove(v)` maintain the current path. In edge mode, v
    // always represents the real edge parent_edge(v).
    template <class Add, class Remove, class Answer>
    void run(
        Add add,
        Remove remove,
        Answer answer,
        int block_size = 0
    ) const {
        bool edge_mode = !_queries.empty() && _queries.front().edge;
        std::vector<char> active(_n, false);
        auto toggle = [&](int tour_index) {
            int vertex = tour[tour_index];
            if (!edge_mode || vertex != root) {
                if (active[vertex]) {
                    remove(vertex);
                } else {
                    add(vertex);
                }
            }
            active[vertex] = !active[vertex];
        };

        _mo.run(
            toggle,
            toggle,
            [&](int query_id) {
                int extra = _queries[query_id].extra;
                if (extra != -1) {
                    assert(!active[extra]);
                    add(extra);
                }
                answer(query_id);
                if (extra != -1) remove(extra);
            },
            block_size
        );
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/range_contour_query.hpp"



#line 7 "graph/tree/range_contour_query.hpp"

#line 1 "graph/tree/rooted_tree.hpp"



#line 7 "graph/tree/rooted_tree.hpp"

#line 9 "graph/tree/rooted_tree.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct RootedTree {
    using cost_type = T;
    using edge_type = m1une::graph::Edge<T>;

    int root;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<T> dist;
    std::vector<int> subtree_size;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> order;
    std::vector<std::vector<int>> up;

   private:
    int _n;
    int _log;

    void check_vertex(int v) const {
        assert(0 <= v && v < _n);
        assert(tin[v] != -1);
    }

   public:
    RootedTree() : root(-1), _n(0), _log(0) {}
    explicit RootedTree(const m1une::graph::Graph<T>& g, int root_ = 0) {
        build(g, root_);
    }

    void build(const m1une::graph::Graph<T>& g, int root_ = 0) {
        _n = g.size();
        root = _n == 0 ? -1 : root_;
        _log = 1;
        while ((1U << _log) <= (unsigned int)(std::max(1, _n))) _log++;

        parent.assign(_n, -1);
        parent_edge.assign(_n, -1);
        depth.assign(_n, 0);
        dist.assign(_n, T(0));
        subtree_size.assign(_n, 0);
        tin.assign(_n, -1);
        tout.assign(_n, -1);
        order.clear();
        order.reserve(_n);
        up.assign(_log, std::vector<int>(_n, -1));

        if (_n == 0) return;
        assert(0 <= root && root < _n);

        struct Frame {
            int v;
            int state;
        };

        std::vector<char> visited(_n, false);
        std::vector<Frame> stack;
        stack.push_back({root, 0});
        visited[root] = true;
        int timer = 0;

        while (!stack.empty()) {
            Frame frame = stack.back();
            stack.pop_back();
            int v = frame.v;
            if (frame.state == 0) {
                tin[v] = timer++;
                order.push_back(v);
                up[0][v] = parent[v];
                for (int k = 1; k < _log; k++) {
                    int p = up[k - 1][v];
                    up[k][v] = p == -1 ? -1 : up[k - 1][p];
                }

                stack.push_back({v, 1});
                const auto& adj = g[v];
                for (int i = int(adj.size()) - 1; i >= 0; i--) {
                    const auto& e = adj[i];
                    if (!e.alive) continue;
                    if (visited[e.to]) continue;
                    visited[e.to] = true;
                    parent[e.to] = v;
                    parent_edge[e.to] = e.id;
                    depth[e.to] = depth[v] + 1;
                    dist[e.to] = dist[v] + e.cost;
                    stack.push_back({e.to, 0});
                }
            } else {
                subtree_size[v]++;
                if (parent[v] != -1) subtree_size[parent[v]] += subtree_size[v];
                tout[v] = timer;
            }
        }
    }

    int size() const {
        return _n;
    }

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

    int log() const {
        return _log;
    }

    bool is_ancestor(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        return tin[u] <= tin[v] && tout[v] <= tout[u];
    }

    bool in_subtree(int v, int u) const {
        return is_ancestor(u, v);
    }

    int kth_ancestor(int v, int k) const {
        check_vertex(v);
        assert(0 <= k);
        int bit = 0;
        while (k > 0 && v != -1) {
            if (k & 1) {
                if (_log <= bit) return -1;
                v = up[bit][v];
            }
            k >>= 1;
            bit++;
        }
        return v;
    }

    int lca(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        if (depth[u] < depth[v]) std::swap(u, v);
        u = kth_ancestor(u, depth[u] - depth[v]);
        if (u == v) return u;
        for (int k = _log - 1; k >= 0; k--) {
            if (up[k][u] != up[k][v]) {
                u = up[k][u];
                v = up[k][v];
            }
        }
        return parent[u];
    }

    int dist_edges(int u, int v) const {
        int w = lca(u, v);
        return depth[u] + depth[v] - 2 * depth[w];
    }

    T dist_cost(int u, int v) const {
        int w = lca(u, v);
        return dist[u] + dist[v] - dist[w] - dist[w];
    }

    int jump(int from, int to, int k) const {
        check_vertex(from);
        check_vertex(to);
        assert(0 <= k);
        int w = lca(from, to);
        int up_len = depth[from] - depth[w];
        int down_len = depth[to] - depth[w];
        if (up_len + down_len < k) return -1;
        if (k <= up_len) return kth_ancestor(from, k);
        return kth_ancestor(to, down_len - (k - up_len));
    }

    std::vector<int> path(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        int w = lca(u, v);
        std::vector<int> a, b;
        for (int x = u; x != w; x = parent[x]) a.push_back(x);
        a.push_back(w);
        for (int x = v; x != w; x = parent[x]) b.push_back(x);
        std::reverse(b.begin(), b.end());
        a.insert(a.end(), b.begin(), b.end());
        return a;
    }

    std::vector<int> path_edges(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        int w = lca(u, v);
        std::vector<int> a, b;
        for (int x = u; x != w; x = parent[x]) a.push_back(parent_edge[x]);
        for (int x = v; x != w; x = parent[x]) b.push_back(parent_edge[x]);
        std::reverse(b.begin(), b.end());
        a.insert(a.end(), b.begin(), b.end());
        return a;
    }

    std::pair<int, int> subtree_range(int v) const {
        check_vertex(v);
        return {tin[v], tout[v]};
    }

    std::vector<int> subtree_vertices(int v) const {
        check_vertex(v);
        return std::vector<int>(order.begin() + tin[v], order.begin() + tout[v]);
    }
};

}  // namespace tree
}  // namespace m1une


#line 13 "graph/tree/range_contour_query.hpp"

namespace m1une {
namespace tree {

namespace internal {

struct RangeContourPathEntry {
    int centroid;
    int distance;
    int subtree;
};

struct RangeContourLayout {
    int n = 0;
    std::vector<std::vector<RangeContourPathEntry>> path;
    std::vector<int> all_size;
    std::vector<int> subtree_size;

    template <class EdgeCost>
    void build(const m1une::graph::Graph<EdgeCost>& graph) {
        n = graph.size();
        path.assign(n, {});
        all_size.assign(n, 0);
        subtree_size.assign(n, 0);
        if (n == 0) return;

#ifndef NDEBUG
        std::vector<int> incidence(graph.edge_count(), 0);
        for (int vertex = 0; vertex < n; vertex++) {
            for (const auto& edge : graph[vertex]) {
                if (!edge.alive) continue;
                assert(0 <= edge.id && edge.id < graph.edge_count());
                incidence[edge.id]++;
            }
        }
        int active_edges = 0;
        for (int count : incidence) {
            if (count == 0) continue;
            assert(count == 2);
            active_edges++;
        }
        assert(active_edges == n - 1);
#endif

        RootedTree<EdgeCost> rooted(graph, 0);
        assert(int(rooted.order.size()) == n);
        CentroidDecomposition<EdgeCost> decomposition(graph);

        for (int vertex = 0; vertex < n; vertex++) {
            int previous = -1;
            for (
                int centroid = vertex;
                centroid != -1;
                centroid = decomposition.parent[centroid]
            ) {
                int distance = rooted.dist_edges(vertex, centroid);
                path[vertex].push_back(
                    RangeContourPathEntry{centroid, distance, previous}
                );
                all_size[centroid] = std::max(
                    all_size[centroid],
                    distance + 1
                );
                if (previous != -1) {
                    subtree_size[previous] = std::max(
                        subtree_size[previous],
                        distance + 1
                    );
                }
                previous = centroid;
            }
        }
    }
};

template <m1une::monoid::IsCommutativeGroup Group>
class RangeContourFenwick {
   public:
    using T = typename Group::value_type;

   private:
    int _n = 0;
    std::vector<T> _data;

    T prefix_product(int right) const {
        T result = Group::id();
        while (right > 0) {
            result = Group::op(result, _data[right]);
            right -= right & -right;
        }
        return result;
    }

   public:
    RangeContourFenwick() : _data(1, Group::id()) {}

    explicit RangeContourFenwick(int n)
        : _n(n), _data(n + 1, Group::id()) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    void apply(int index, const T& value) {
        assert(0 <= index && index < _n);
        for (index++; index <= _n; index += index & -index) {
            _data[index] = Group::op(_data[index], value);
        }
    }

    T product(int left, int right) const {
        left = std::max(left, 0);
        right = std::min(right, _n);
        if (right <= left) return Group::id();
        return Group::op(
            Group::inv(prefix_product(left)),
            prefix_product(right)
        );
    }

    void range_apply(int left, int right, const T& value) {
        left = std::max(left, 0);
        right = std::min(right, _n);
        if (right <= left) return;
        apply(left, value);
        if (right < _n) apply(right, Group::inv(value));
    }

    T get(int index) const {
        assert(0 <= index && index < _n);
        return prefix_product(index + 1);
    }
};

}  // namespace internal

template <m1une::monoid::IsCommutativeGroup Group>
class VertexApplyRangeContourProduct {
   public:
    using T = typename Group::value_type;

   private:
    internal::RangeContourLayout _layout;
    std::vector<T> _value;
    std::vector<internal::RangeContourFenwick<Group>> _all;
    std::vector<internal::RangeContourFenwick<Group>> _subtree;

    void check_vertex(int vertex) const {
        assert(0 <= vertex && vertex < size());
    }

   public:
    VertexApplyRangeContourProduct() = default;

    template <class EdgeCost>
    explicit VertexApplyRangeContourProduct(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<T>& initial = {}
    ) {
        build(graph, initial);
    }

    template <class EdgeCost>
    void build(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<T>& initial = {}
    ) {
        assert(initial.empty() || int(initial.size()) == graph.size());
        _layout.build(graph);
        const int n = _layout.n;
        _value.assign(n, Group::id());
        _all.assign(n, internal::RangeContourFenwick<Group>());
        _subtree.assign(n, internal::RangeContourFenwick<Group>());
        for (int index = 0; index < n; index++) {
            _all[index] =
                internal::RangeContourFenwick<Group>(_layout.all_size[index]);
            _subtree[index] =
                internal::RangeContourFenwick<Group>(
                    _layout.subtree_size[index]
                );
        }
        if (!initial.empty()) {
            for (int vertex = 0; vertex < n; vertex++) {
                apply(vertex, initial[vertex]);
            }
        }
    }

    int size() const {
        return _layout.n;
    }

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

    T get(int vertex) const {
        check_vertex(vertex);
        return _value[vertex];
    }

    void apply(int vertex, const T& value) {
        check_vertex(vertex);
        _value[vertex] = Group::op(_value[vertex], value);
        for (const auto& entry : _layout.path[vertex]) {
            _all[entry.centroid].apply(entry.distance, value);
            if (entry.subtree != -1) {
                _subtree[entry.subtree].apply(entry.distance, value);
            }
        }
    }

    void set(int vertex, const T& value) {
        check_vertex(vertex);
        apply(vertex, Group::op(Group::inv(_value[vertex]), value));
    }

    T prod(int vertex, int left_distance, int right_distance) const {
        check_vertex(vertex);
        assert(0 <= left_distance && left_distance <= right_distance);
        T result = Group::id();
        for (const auto& entry : _layout.path[vertex]) {
            int left = left_distance - entry.distance;
            int right = right_distance - entry.distance;
            result = Group::op(
                result,
                _all[entry.centroid].product(left, right)
            );
            if (entry.subtree != -1) {
                result = Group::op(
                    result,
                    Group::inv(
                        _subtree[entry.subtree].product(left, right)
                    )
                );
            }
        }
        return result;
    }
};

template <m1une::monoid::IsCommutativeGroup Group>
class VertexGetRangeContourApply {
   public:
    using T = typename Group::value_type;

   private:
    internal::RangeContourLayout _layout;
    std::vector<T> _base;
    std::vector<internal::RangeContourFenwick<Group>> _all;
    std::vector<internal::RangeContourFenwick<Group>> _subtree;

    void check_vertex(int vertex) const {
        assert(0 <= vertex && vertex < size());
    }

   public:
    VertexGetRangeContourApply() = default;

    template <class EdgeCost>
    explicit VertexGetRangeContourApply(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<T>& initial = {}
    ) {
        build(graph, initial);
    }

    template <class EdgeCost>
    void build(
        const m1une::graph::Graph<EdgeCost>& graph,
        const std::vector<T>& initial = {}
    ) {
        assert(initial.empty() || int(initial.size()) == graph.size());
        _layout.build(graph);
        const int n = _layout.n;
        _base = initial.empty() ? std::vector<T>(n, Group::id()) : initial;
        _all.assign(n, internal::RangeContourFenwick<Group>());
        _subtree.assign(n, internal::RangeContourFenwick<Group>());
        for (int index = 0; index < n; index++) {
            _all[index] =
                internal::RangeContourFenwick<Group>(_layout.all_size[index]);
            _subtree[index] =
                internal::RangeContourFenwick<Group>(
                    _layout.subtree_size[index]
                );
        }
    }

    int size() const {
        return _layout.n;
    }

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

    T get(int vertex) const {
        check_vertex(vertex);
        T result = _base[vertex];
        for (const auto& entry : _layout.path[vertex]) {
            result = Group::op(
                result,
                _all[entry.centroid].get(entry.distance)
            );
            if (entry.subtree != -1) {
                result = Group::op(
                    result,
                    Group::inv(
                        _subtree[entry.subtree].get(entry.distance)
                    )
                );
            }
        }
        return result;
    }

    void point_apply(int vertex, const T& value) {
        check_vertex(vertex);
        _base[vertex] = Group::op(_base[vertex], value);
    }

    void set(int vertex, const T& value) {
        check_vertex(vertex);
        _base[vertex] = Group::op(
            _base[vertex],
            Group::op(Group::inv(get(vertex)), value)
        );
    }

    void apply(
        int vertex,
        int left_distance,
        int right_distance,
        const T& value
    ) {
        check_vertex(vertex);
        assert(0 <= left_distance && left_distance <= right_distance);
        for (const auto& entry : _layout.path[vertex]) {
            int left = left_distance - entry.distance;
            int right = right_distance - entry.distance;
            _all[entry.centroid].range_apply(left, right, value);
            if (entry.subtree != -1) {
                _subtree[entry.subtree].range_apply(left, right, value);
            }
        }
    }
};

template <class T>
class VertexAddRangeContourSum
    : public VertexApplyRangeContourProduct<m1une::monoid::Add<T>> {
   private:
    using Base = VertexApplyRangeContourProduct<m1une::monoid::Add<T>>;

   public:
    using Base::Base;

    void add(int vertex, const T& delta) {
        Base::apply(vertex, delta);
    }

    T sum(int vertex, int left_distance, int right_distance) const {
        return Base::prod(vertex, left_distance, right_distance);
    }
};

template <class T>
class VertexGetRangeContourAdd
    : public VertexGetRangeContourApply<m1une::monoid::Add<T>> {
   private:
    using Base = VertexGetRangeContourApply<m1une::monoid::Add<T>>;

   public:
    using Base::Base;

    void add(int vertex, const T& delta) {
        Base::point_apply(vertex, delta);
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/rerooting_dp.hpp"



#line 5 "graph/tree/rerooting_dp.hpp"

#line 7 "graph/tree/rerooting_dp.hpp"

namespace m1une {
namespace tree {

template <class T, class DP, class Merge, class AddVertex, class AddEdge>
std::vector<DP> rerooting_dp(const m1une::graph::Graph<T>& g, DP id, Merge merge, AddVertex add_vertex,
                             AddEdge add_edge) {
    int n = g.size();
    std::vector<int> parent(n, -2), parent_edge(n, -1), order;
    order.reserve(n);
    for (int root = 0; root < n; root++) {
        if (parent[root] != -2) continue;
        parent[root] = -1;
        std::vector<int> stack = {root};
        while (!stack.empty()) {
            int v = stack.back();
            stack.pop_back();
            order.push_back(v);
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                if (parent[e.to] != -2) continue;
                parent[e.to] = v;
                parent_edge[e.to] = e.id;
                stack.push_back(e.to);
            }
        }
    }

    std::vector<DP> down(n, id), outside(n, id), answer(n, id);
    for (int i = n - 1; i >= 0; i--) {
        int v = order[i];
        DP acc = id;
        for (const auto& e : g[v]) {
            if (!e.alive) continue;
            if (parent[e.to] != v) continue;
            acc = merge(acc, add_edge(down[e.to], e));
        }
        down[v] = add_vertex(acc, v);
    }

    for (int v : order) {
        int d = int(g[v].size());
        std::vector<DP> contrib(d, id);
        for (int i = 0; i < d; i++) {
            const auto& e = g[v][i];
            if (!e.alive) continue;
            if (parent[e.to] == v) {
                contrib[i] = add_edge(down[e.to], e);
            } else if (parent[v] == e.to && parent_edge[v] == e.id) {
                contrib[i] = add_edge(outside[v], e);
            }
        }

        std::vector<DP> pref(d + 1, id), suff(d + 1, id);
        for (int i = 0; i < d; i++) pref[i + 1] = merge(pref[i], contrib[i]);
        for (int i = d - 1; i >= 0; i--) suff[i] = merge(contrib[i], suff[i + 1]);
        answer[v] = add_vertex(pref[d], v);

        for (int i = 0; i < d; i++) {
            const auto& e = g[v][i];
            if (!e.alive) continue;
            if (parent[e.to] != v) continue;
            outside[e.to] = add_vertex(merge(pref[i], suff[i + 1]), v);
        }
    }

    return answer;
}

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/rerooting_static_top_tree.hpp"



#line 6 "graph/tree/rerooting_static_top_tree.hpp"
#include <optional>
#line 10 "graph/tree/rerooting_static_top_tree.hpp"

#line 12 "graph/tree/rerooting_static_top_tree.hpp"

namespace m1une {
namespace tree {

namespace internal {

enum class RerootingStaticTopTreeNodeType {
    Compress,
    Rake,
    AddEdge,
    AddVertex,
};

enum class RerootingStaticTopTreeStepType {
    CompressLower,
    CompressUpper,
    AddEdge,
    RakeLeft,
    RakeRight,
    AddVertex,
};

}  // namespace internal

template <class T, class Vertex, class Path, class Point, class CompressDown, class CompressUp, class Rake,
          class AddEdgeDown, class AddEdgeUp, class AddVertex>
struct RerootingStaticTopTree {
    using cost_type = T;
    using vertex_type = Vertex;
    using path_type = Path;
    using point_type = Point;
    using edge_type = m1une::graph::Edge<T>;
    using node_type = internal::RerootingStaticTopTreeNodeType;
    using step_type = internal::RerootingStaticTopTreeStepType;

    struct Node {
        node_type type;
        int left = -1;
        int right = -1;
        int parent = -1;
        int vertex = -1;
        edge_type edge;
        int size = 0;
        int height = 1;
        std::optional<Path> path_down;
        std::optional<Path> path_up;
        std::optional<Point> point;
    };

    struct RerootingStep {
        step_type type;
        int node = -1;
        int sibling = -1;
        int vertex = -1;
        edge_type edge;
    };

   private:
    int _n;
    int _root;
    int _root_node;
    Point _point_id;
    CompressDown _compress_down;
    CompressUp _compress_up;
    Rake _rake;
    AddEdgeDown _add_edge_down;
    AddEdgeUp _add_edge_up;
    AddVertex _add_vertex;
    std::vector<Vertex> _values;
    std::vector<Node> _nodes;
    std::vector<int> _vertex_node;
    std::vector<int> _edge_node;
    std::vector<int> _parent;
    std::vector<int> _subtree_size;
    std::vector<int> _heavy;
    std::vector<edge_type> _heavy_edge;
    std::vector<std::vector<edge_type>> _children;

    static edge_type reversed_edge(edge_type e) {
        std::swap(e.from, e.to);
        return e;
    }

    const Path& node_path_down(int node) const {
        assert(0 <= node && node < int(_nodes.size()));
        assert(_nodes[node].path_down.has_value());
        return *_nodes[node].path_down;
    }

    const Path& node_path_up(int node) const {
        assert(0 <= node && node < int(_nodes.size()));
        assert(_nodes[node].path_up.has_value());
        return *_nodes[node].path_up;
    }

    const Point& node_point(int node) const {
        assert(0 <= node && node < int(_nodes.size()));
        assert(_nodes[node].point.has_value());
        return *_nodes[node].point;
    }

    void set_parent(int child, int parent) {
        if (child != -1) _nodes[child].parent = parent;
    }

    void recompute(int node) {
        auto& x = _nodes[node];
        if (x.type == node_type::Compress) {
            x.path_down = _compress_down(node_path_down(x.left), node_path_down(x.right), x.edge);
            x.path_up = _compress_up(node_path_up(x.right), node_path_up(x.left), reversed_edge(x.edge));
        } else if (x.type == node_type::Rake) {
            x.point = _rake(node_point(x.left), node_point(x.right));
        } else if (x.type == node_type::AddEdge) {
            x.point = _add_edge_down(node_path_down(x.left), x.edge);
        } else {
            const Point& side = x.left == -1 ? _point_id : node_point(x.left);
            Path path = _add_vertex(side, _values[x.vertex], x.vertex);
            x.path_down = path;
            x.path_up = std::move(path);
        }
    }

    int new_node(Node node) {
        int id = int(_nodes.size());
        _nodes.push_back(std::move(node));
        set_parent(_nodes[id].left, id);
        set_parent(_nodes[id].right, id);
        recompute(id);
        return id;
    }

    int new_compress(int left, int right, edge_type edge) {
        Node node;
        node.type = node_type::Compress;
        node.left = left;
        node.right = right;
        node.edge = edge;
        node.size = _nodes[left].size + _nodes[right].size;
        node.height = std::max(_nodes[left].height, _nodes[right].height) + 1;
        int id = new_node(std::move(node));
        if (0 <= edge.id && edge.id < int(_edge_node.size())) _edge_node[edge.id] = id;
        return id;
    }

    int new_rake(int left, int right) {
        Node node;
        node.type = node_type::Rake;
        node.left = left;
        node.right = right;
        node.size = _nodes[left].size + _nodes[right].size;
        node.height = std::max(_nodes[left].height, _nodes[right].height) + 1;
        return new_node(std::move(node));
    }

    int new_add_edge(int child, edge_type edge) {
        Node node;
        node.type = node_type::AddEdge;
        node.left = child;
        node.edge = edge;
        node.size = _nodes[child].size;
        node.height = _nodes[child].height + 1;
        int id = new_node(std::move(node));
        if (0 <= edge.id && edge.id < int(_edge_node.size())) _edge_node[edge.id] = id;
        return id;
    }

    int new_add_vertex(int side, int vertex) {
        Node node;
        node.type = node_type::AddVertex;
        node.left = side;
        node.vertex = vertex;
        node.size = 1 + (side == -1 ? 0 : _nodes[side].size);
        node.height = 1 + (side == -1 ? 0 : _nodes[side].height);
        int id = new_node(std::move(node));
        _vertex_node[vertex] = id;
        return id;
    }

    int weighted_split(const std::vector<int>& nodes, int l, int r) const {
        int total = 0;
        for (int i = l; i < r; i++) total += _nodes[nodes[i]].size;
        int left_sum = 0;
        for (int i = l; i + 1 < r; i++) {
            left_sum += _nodes[nodes[i]].size;
            if (2 * left_sum >= total) return i + 1;
        }
        return r - 1;
    }

    int build_rake(const std::vector<int>& nodes, int l, int r) {
        if (l == r) return -1;
        if (l + 1 == r) return nodes[l];
        int m = weighted_split(nodes, l, r);
        return new_rake(build_rake(nodes, l, m), build_rake(nodes, m, r));
    }

    int build_compress(const std::vector<int>& nodes, const std::vector<edge_type>& edges, int l, int r) {
        if (l + 1 == r) return nodes[l];
        int m = weighted_split(nodes, l, r);
        return new_compress(build_compress(nodes, edges, l, m), build_compress(nodes, edges, m, r), edges[m - 1]);
    }

    int build_vertex(int v) {
        std::vector<int> side_nodes;
        for (const auto& e : _children[v]) {
            if (e.to == _heavy[v]) continue;
            int child_path = build_path(e.to);
            side_nodes.push_back(new_add_edge(child_path, e));
        }
        return new_add_vertex(build_rake(side_nodes, 0, int(side_nodes.size())), v);
    }

    int build_path(int start) {
        std::vector<int> path_nodes;
        std::vector<edge_type> path_edges;
        for (int v = start; v != -1; v = _heavy[v]) {
            path_nodes.push_back(build_vertex(v));
            if (_heavy[v] != -1) path_edges.push_back(_heavy_edge[v]);
        }
        return build_compress(path_nodes, path_edges, 0, int(path_nodes.size()));
    }

    void recompute_up(int node) {
        while (node != -1) {
            recompute(node);
            node = _nodes[node].parent;
        }
    }

   public:
    RerootingStaticTopTree(const m1une::graph::Graph<T>& g, const std::vector<Vertex>& values, Point point_id,
                           CompressDown compress_down, CompressUp compress_up, Rake rake,
                           AddEdgeDown add_edge_down, AddEdgeUp add_edge_up, AddVertex add_vertex, int root = 0)
        : _n(g.size()),
          _root(_n == 0 ? -1 : root),
          _root_node(-1),
          _point_id(std::move(point_id)),
          _compress_down(std::move(compress_down)),
          _compress_up(std::move(compress_up)),
          _rake(std::move(rake)),
          _add_edge_down(std::move(add_edge_down)),
          _add_edge_up(std::move(add_edge_up)),
          _add_vertex(std::move(add_vertex)),
          _values(values) {
        build(g, root);
    }

    void build(const m1une::graph::Graph<T>& g, int root = 0) {
        _n = g.size();
        _root = _n == 0 ? -1 : root;
        assert(int(_values.size()) == _n);
        _nodes.clear();
        _vertex_node.assign(_n, -1);
        _edge_node.assign(g.edge_count(), -1);
        _parent.assign(_n, -2);
        _subtree_size.assign(_n, 1);
        _heavy.assign(_n, -1);
        _heavy_edge.assign(_n, edge_type());
        _children.assign(_n, {});
        _root_node = -1;

        if (_n == 0) return;
        assert(0 <= root && root < _n);
        assert(int(g.edges().size()) == _n - 1);

        std::vector<int> order;
        order.reserve(_n);
        std::vector<int> stack = {root};
        _parent[root] = -1;
        while (!stack.empty()) {
            int v = stack.back();
            stack.pop_back();
            order.push_back(v);
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                if (_parent[e.to] != -2) continue;
                _parent[e.to] = v;
                _children[v].push_back(e);
                stack.push_back(e.to);
            }
        }
        assert(int(order.size()) == _n);

        for (int i = int(order.size()) - 1; i >= 0; i--) {
            int v = order[i];
            for (const auto& e : _children[v]) {
                _subtree_size[v] += _subtree_size[e.to];
                if (_heavy[v] == -1 || _subtree_size[_heavy[v]] < _subtree_size[e.to]) {
                    _heavy[v] = e.to;
                    _heavy_edge[v] = e;
                }
            }
        }

        _root_node = build_path(root);
    }

    int size() const {
        return _n;
    }

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

    int root() const {
        return _root;
    }

    int root_node() const {
        return _root_node;
    }

    int node_count() const {
        return int(_nodes.size());
    }

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

    const std::vector<Node>& nodes() const {
        return _nodes;
    }

    const Node& node(int id) const {
        assert(0 <= id && id < int(_nodes.size()));
        return _nodes[id];
    }

    int parent_node(int id) const {
        return node(id).parent;
    }

    int vertex_node(int v) const {
        assert(0 <= v && v < _n);
        return _vertex_node[v];
    }

    int local_point_node(int v) const {
        int id = vertex_node(v);
        assert(_nodes[id].type == node_type::AddVertex);
        return _nodes[id].left;
    }

    const Point& local_point(int v) const {
        int point_node = local_point_node(v);
        return point_node == -1 ? _point_id : node_point(point_node);
    }

    const Vertex& get(int v) const {
        assert(0 <= v && v < _n);
        return _values[v];
    }

    const Vertex& operator[](int v) const {
        return get(v);
    }

    void set(int v, const Vertex& value) {
        assert(0 <= v && v < _n);
        assert(_vertex_node[v] != -1);
        _values[v] = value;
        recompute_up(_vertex_node[v]);
    }

    void set(int v, Vertex&& value) {
        assert(0 <= v && v < _n);
        assert(_vertex_node[v] != -1);
        _values[v] = std::move(value);
        recompute_up(_vertex_node[v]);
    }

    void set_edge_cost(int edge_id, T cost) {
        assert(0 <= edge_id && edge_id < int(_edge_node.size()));
        int node = _edge_node[edge_id];
        assert(node != -1);
        _nodes[node].edge.cost = cost;
        recompute_up(node);
    }

    const Path& path_down(int node_id) const {
        return node_path_down(node_id);
    }

    const Path& path_up(int node_id) const {
        return node_path_up(node_id);
    }

    const Point& point(int node_id) const {
        return node_point(node_id);
    }

    const Path& all_prod_down() const {
        assert(_root_node != -1);
        return path_down(_root_node);
    }

    const Path& all_prod_up() const {
        assert(_root_node != -1);
        return path_up(_root_node);
    }

    const Point& point_id() const {
        return _point_id;
    }

    template <class F>
    void for_each_rerooting_step(int v, F&& f) const {
        assert(0 <= v && v < _n);
        int cur = _vertex_node[v];
        assert(cur != -1);
        while (_nodes[cur].parent != -1) {
            int par = _nodes[cur].parent;
            const auto& p = _nodes[par];
            RerootingStep step;
            step.node = par;
            if (p.type == node_type::Compress) {
                step.edge = p.edge;
                if (p.left == cur) {
                    step.type = step_type::CompressLower;
                    step.sibling = p.right;
                } else {
                    assert(p.right == cur);
                    step.type = step_type::CompressUpper;
                    step.sibling = p.left;
                }
            } else if (p.type == node_type::Rake) {
                if (p.left == cur) {
                    step.type = step_type::RakeRight;
                    step.sibling = p.right;
                } else {
                    assert(p.right == cur);
                    step.type = step_type::RakeLeft;
                    step.sibling = p.left;
                }
            } else if (p.type == node_type::AddEdge) {
                assert(p.left == cur);
                step.type = step_type::AddEdge;
                step.edge = p.edge;
            } else {
                assert(p.type == node_type::AddVertex);
                assert(p.left == cur);
                step.type = step_type::AddVertex;
                step.vertex = p.vertex;
            }
            f(step);
            cur = par;
        }
    }

    std::vector<RerootingStep> rerooting_steps(int v) const {
        std::vector<RerootingStep> result;
        int cur = vertex_node(v);
        int depth = 0;
        while (_nodes[cur].parent != -1) {
            cur = _nodes[cur].parent;
            depth++;
        }
        result.reserve(depth);
        for_each_rerooting_step(v, [&](const RerootingStep& step) {
            result.push_back(step);
        });
        return result;
    }

    template <class Folder>
    auto fold_rerooting(int v, Folder folder) const {
        folder.start(v, _values[v], local_point(v));
        for_each_rerooting_step(v, [&](const RerootingStep& step) {
            if (step.type == step_type::CompressLower) {
                folder.compress_lower(path_down(step.sibling), step.edge);
            } else if (step.type == step_type::CompressUpper) {
                folder.compress_upper(path_up(step.sibling), reversed_edge(step.edge));
            } else if (step.type == step_type::AddEdge) {
                folder.add_edge(reversed_edge(step.edge));
            } else if (step.type == step_type::RakeLeft) {
                folder.rake_left(point(step.sibling));
            } else if (step.type == step_type::RakeRight) {
                folder.rake_right(point(step.sibling));
            } else {
                folder.add_vertex(step.vertex, _values[step.vertex]);
            }
        });
        return folder.result();
    }

    Path compress_down(const Path& upper, const Path& lower, edge_type edge) const {
        return _compress_down(upper, lower, edge);
    }

    Path compress_up(const Path& lower, const Path& upper, edge_type edge) const {
        return _compress_up(lower, upper, edge);
    }

    Point rake(const Point& left, const Point& right) const {
        return _rake(left, right);
    }

    Point add_edge_down(const Path& path, edge_type edge) const {
        return _add_edge_down(path, edge);
    }

    Point add_edge_up(const Path& path, edge_type edge) const {
        return _add_edge_up(path, edge);
    }

    Path add_vertex(const Point& side, const Vertex& value, int vertex) const {
        return _add_vertex(side, value, vertex);
    }

    static edge_type reverse_edge(edge_type edge) {
        return reversed_edge(edge);
    }
};

template <class T, class Vertex, class Point, class CompressDown, class CompressUp, class Rake, class AddEdgeDown,
          class AddEdgeUp, class AddVertex>
RerootingStaticTopTree(const m1une::graph::Graph<T>&, const std::vector<Vertex>&, Point, CompressDown, CompressUp,
                       Rake, AddEdgeDown, AddEdgeUp, AddVertex, int)
    -> RerootingStaticTopTree<T, Vertex, std::invoke_result_t<AddVertex, Point, Vertex, int>, Point, CompressDown,
                              CompressUp, Rake, AddEdgeDown, AddEdgeUp, AddVertex>;

template <class T, class Vertex, class Point, class CompressDown, class CompressUp, class Rake, class AddEdgeDown,
          class AddEdgeUp, class AddVertex>
RerootingStaticTopTree(const m1une::graph::Graph<T>&, const std::vector<Vertex>&, Point, CompressDown, CompressUp,
                       Rake, AddEdgeDown, AddEdgeUp, AddVertex)
    -> RerootingStaticTopTree<T, Vertex, std::invoke_result_t<AddVertex, Point, Vertex, int>, Point, CompressDown,
                              CompressUp, Rake, AddEdgeDown, AddEdgeUp, AddVertex>;

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/sparse_table_lca.hpp"



#line 9 "graph/tree/sparse_table_lca.hpp"

#line 1 "ds/range_query/sparse_table.hpp"



#include <bit>
#line 9 "ds/range_query/sparse_table.hpp"

#line 11 "ds/range_query/sparse_table.hpp"

namespace m1une {
namespace ds {

// A Sparse Table utilizing C++20 Concepts for type safety.
// It requires a Monoid struct that satisfies `m1une::monoid::IsMonoid`.
// [IMPORTANT] For O(1) range queries to work correctly, the monoid operation MUST be idempotent.
// i.e., Monoid::op(x, x) == x must hold (e.g., Min, Max, GCD, Bitwise AND/OR).
template <m1une::monoid::IsMonoid Monoid>
struct SparseTable {
    using T = typename Monoid::value_type;

   private:
    int _n;
    std::vector<std::vector<T>> _st;

   public:
    // Constructs an empty sparse table.
    SparseTable() : _n(0) {}

    // Constructs a sparse table from an existing vector in O(N log N) time.
    explicit SparseTable(const std::vector<T>& v) : _n(int(v.size())) {
        if (_n == 0) return;

        // Compute the maximum power of 2 needed
        int max_log = std::bit_width((unsigned int)_n);
        _st.assign(max_log, std::vector<T>(_n));

        // Initialize the base level
        for (int i = 0; i < _n; i++) {
            _st[0][i] = v[i];
        }

        // Build the sparse table
        for (int k = 1; k < max_log; k++) {
            for (int i = 0; i + (1 << k) <= _n; i++) {
                _st[k][i] = Monoid::op(_st[k - 1][i], _st[k - 1][i + (1 << (k - 1))]);
            }
        }
    }
    explicit SparseTable(std::vector<T>&& v) : _n(int(v.size())) {
        if (_n == 0) return;

        int max_log = std::bit_width((unsigned int)_n);
        _st.assign(max_log, std::vector<T>(_n));

        for (int i = 0; i < _n; i++) {
            _st[0][i] = std::move(v[i]);
        }

        for (int k = 1; k < max_log; k++) {
            for (int i = 0; i + (1 << k) <= _n; i++) {
                _st[k][i] = Monoid::op(_st[k - 1][i], _st[k - 1][i + (1 << (k - 1))]);
            }
        }
    }

    // Constructs a sparse table from a vector of a different type U.
    // It automatically adapts to the Monoid's initialization requirements:
    // 1. Monoid::make(val) if it exists.
    // 2. Monoid::make(val, index) if the monoid requires global indices.
    // 3. static_cast<T>(val) as a fallback for simple monoids.
    template <typename U>
    requires (!std::same_as<U, T>) && (
        requires(U x) { Monoid::make(x); } ||
        requires(U x, int i) { Monoid::make(x, i); } ||
        std::convertible_to<U, T>
    )
    explicit SparseTable(const std::vector<U>& v) : _n(int(v.size())) {
        if (_n == 0) return;

        int max_log = std::bit_width((unsigned int)_n);
        _st.assign(max_log, std::vector<T>(_n));

        // Compile-time branching based on the available make() signature
        for (int i = 0; i < _n; i++) {
            if constexpr (requires(U x) { Monoid::make(x); }) {
                _st[0][i] = Monoid::make(v[i]);
            } else if constexpr (requires(U x, int idx) { Monoid::make(x, idx); }) {
                _st[0][i] = Monoid::make(v[i], i);
            } else {
                _st[0][i] = static_cast<T>(v[i]);
            }
        }
        for (int k = 1; k < max_log; k++) {
            for (int i = 0; i + (1 << k) <= _n; i++) {
                _st[k][i] = Monoid::op(_st[k - 1][i], _st[k - 1][i + (1 << (k - 1))]);
            }
        }
    }

    // Returns the product (result of the monoid operation) in the range [l, r) in O(1) time.
    // Requires the monoid operation to be idempotent.
    T prod(int l, int r) const {
        assert(0 <= l && l <= r && r <= _n);
        if (l == r) return Monoid::id();

        // Calculate the largest power of 2 less than or equal to the interval length
        int k = std::bit_width((unsigned int)(r - l)) - 1;
        return Monoid::op(_st[k][l], _st[k][r - (1 << k)]);
    }
};

}  // namespace ds
}  // namespace m1une


#line 12 "graph/tree/sparse_table_lca.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct SparseTableLca {
    using cost_type = T;
    using edge_type = m1une::graph::Edge<T>;

    int root;
    std::vector<int> parent;
    std::vector<int> parent_edge;
    std::vector<int> depth;
    std::vector<T> dist;
    std::vector<int> subtree_size;
    std::vector<int> tin;
    std::vector<int> tout;
    std::vector<int> order;
    std::vector<int> first;
    std::vector<int> euler;

   private:
    struct RmqNode {
        int depth;
        int vertex;
    };

    struct RmqMonoid {
        using value_type = RmqNode;

        static value_type id() {
            return {std::numeric_limits<int>::max(), -1};
        }

        static value_type op(const value_type& a, const value_type& b) {
            if (a.depth != b.depth) return a.depth < b.depth ? a : b;
            return a.vertex < b.vertex ? a : b;
        }
    };

    int _n;
    m1une::ds::SparseTable<RmqMonoid> _st;

    void check_vertex(int v) const {
        assert(0 <= v && v < _n);
        assert(first[v] != -1);
    }

   public:
    SparseTableLca() : root(-1), _n(0) {}
    explicit SparseTableLca(const m1une::graph::Graph<T>& g, int root_ = 0) {
        build(g, root_);
    }

    void build(const m1une::graph::Graph<T>& g, int root_ = 0) {
        _n = g.size();
        root = _n == 0 ? -1 : root_;
        parent.assign(_n, -2);
        parent_edge.assign(_n, -1);
        depth.assign(_n, 0);
        dist.assign(_n, T(0));
        subtree_size.assign(_n, 0);
        tin.assign(_n, -1);
        tout.assign(_n, -1);
        order.clear();
        order.reserve(_n);
        first.assign(_n, -1);
        euler.clear();
        euler.reserve(std::max(0, 2 * _n - 1));
        _st = m1une::ds::SparseTable<RmqMonoid>();

        if (_n == 0) return;
        assert(0 <= root && root < _n);

        std::vector<int> it(_n, 0);
        std::vector<char> visited(_n, false);
        std::vector<int> stack = {root};
        visited[root] = true;
        parent[root] = -1;

        int timer = 0;
        tin[root] = timer++;
        order.push_back(root);
        first[root] = 0;
        euler.push_back(root);

        while (!stack.empty()) {
            int v = stack.back();
            if (it[v] < int(g[v].size())) {
                const auto& e = g[v][it[v]++];
                if (!e.alive) continue;
                if (visited[e.to]) continue;
                visited[e.to] = true;
                parent[e.to] = v;
                parent_edge[e.to] = e.id;
                depth[e.to] = depth[v] + 1;
                dist[e.to] = dist[v] + e.cost;
                tin[e.to] = timer++;
                order.push_back(e.to);
                first[e.to] = int(euler.size());
                euler.push_back(e.to);
                stack.push_back(e.to);
            } else {
                subtree_size[v]++;
                if (parent[v] != -1) subtree_size[parent[v]] += subtree_size[v];
                tout[v] = timer;
                stack.pop_back();
                if (!stack.empty()) euler.push_back(stack.back());
            }
        }

        std::vector<RmqNode> rmq;
        rmq.reserve(euler.size());
        for (int v : euler) rmq.push_back({depth[v], v});
        _st = m1une::ds::SparseTable<RmqMonoid>(std::move(rmq));
    }

    int size() const {
        return _n;
    }

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

    bool is_ancestor(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        return tin[u] <= tin[v] && tout[v] <= tout[u];
    }

    bool in_subtree(int v, int u) const {
        return is_ancestor(u, v);
    }

    int lca(int u, int v) const {
        check_vertex(u);
        check_vertex(v);
        int l = first[u], r = first[v];
        if (l > r) std::swap(l, r);
        return _st.prod(l, r + 1).vertex;
    }

    int dist_edges(int u, int v) const {
        int w = lca(u, v);
        return depth[u] + depth[v] - 2 * depth[w];
    }

    T dist_cost(int u, int v) const {
        int w = lca(u, v);
        return dist[u] + dist[v] - dist[w] - dist[w];
    }

    std::pair<int, int> subtree_range(int v) const {
        check_vertex(v);
        return {tin[v], tout[v]};
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/static_top_tree.hpp"



#line 10 "graph/tree/static_top_tree.hpp"

#line 12 "graph/tree/static_top_tree.hpp"

namespace m1une {
namespace tree {

namespace internal {

enum class StaticTopTreeNodeType {
    Compress,
    Rake,
    AddEdge,
    AddVertex,
};

}  // namespace internal

template <class T, class Vertex, class Path, class Point, class Compress, class Rake, class AddEdge,
          class AddVertex>
struct StaticTopTree {
    using cost_type = T;
    using vertex_type = Vertex;
    using path_type = Path;
    using point_type = Point;
    using edge_type = m1une::graph::Edge<T>;

   private:
    struct Node {
        internal::StaticTopTreeNodeType type;
        int left = -1;
        int right = -1;
        int parent = -1;
        int vertex = -1;
        edge_type edge;
        int size = 0;
        int height = 1;
        std::optional<Path> path;
        std::optional<Point> point;
    };

    int _n;
    int _root;
    int _root_node;
    Point _point_id;
    Compress _compress;
    Rake _rake;
    AddEdge _add_edge;
    AddVertex _add_vertex;
    std::vector<Vertex> _values;
    std::vector<Node> _nodes;
    std::vector<int> _vertex_node;
    std::vector<int> _edge_node;
    std::vector<int> _parent;
    std::vector<int> _subtree_size;
    std::vector<int> _heavy;
    std::vector<edge_type> _heavy_edge;
    std::vector<std::vector<edge_type>> _children;

    const Path& path_value(int node) const {
        assert(0 <= node && node < int(_nodes.size()));
        assert(_nodes[node].path.has_value());
        return *_nodes[node].path;
    }

    const Point& point_value(int node) const {
        assert(0 <= node && node < int(_nodes.size()));
        assert(_nodes[node].point.has_value());
        return *_nodes[node].point;
    }

    void set_parent(int child, int parent) {
        if (child != -1) _nodes[child].parent = parent;
    }

    void recompute(int node) {
        auto& x = _nodes[node];
        if (x.type == internal::StaticTopTreeNodeType::Compress) {
            x.path = _compress(path_value(x.left), path_value(x.right), x.edge);
        } else if (x.type == internal::StaticTopTreeNodeType::Rake) {
            x.point = _rake(point_value(x.left), point_value(x.right));
        } else if (x.type == internal::StaticTopTreeNodeType::AddEdge) {
            x.point = _add_edge(path_value(x.left), x.edge);
        } else {
            const Point& side = x.left == -1 ? _point_id : point_value(x.left);
            x.path = _add_vertex(side, _values[x.vertex], x.vertex);
        }
    }

    int new_node(Node node) {
        int id = int(_nodes.size());
        _nodes.push_back(std::move(node));
        set_parent(_nodes[id].left, id);
        set_parent(_nodes[id].right, id);
        recompute(id);
        return id;
    }

    int new_compress(int left, int right, edge_type edge) {
        Node node;
        node.type = internal::StaticTopTreeNodeType::Compress;
        node.left = left;
        node.right = right;
        node.edge = edge;
        node.size = _nodes[left].size + _nodes[right].size;
        node.height = std::max(_nodes[left].height, _nodes[right].height) + 1;
        int id = new_node(std::move(node));
        if (0 <= edge.id && edge.id < int(_edge_node.size())) _edge_node[edge.id] = id;
        return id;
    }

    int new_rake(int left, int right) {
        Node node;
        node.type = internal::StaticTopTreeNodeType::Rake;
        node.left = left;
        node.right = right;
        node.size = _nodes[left].size + _nodes[right].size;
        node.height = std::max(_nodes[left].height, _nodes[right].height) + 1;
        return new_node(std::move(node));
    }

    int new_add_edge(int child, edge_type edge) {
        Node node;
        node.type = internal::StaticTopTreeNodeType::AddEdge;
        node.left = child;
        node.edge = edge;
        node.size = _nodes[child].size;
        node.height = _nodes[child].height + 1;
        int id = new_node(std::move(node));
        if (0 <= edge.id && edge.id < int(_edge_node.size())) _edge_node[edge.id] = id;
        return id;
    }

    int new_add_vertex(int side, int vertex) {
        Node node;
        node.type = internal::StaticTopTreeNodeType::AddVertex;
        node.left = side;
        node.vertex = vertex;
        node.size = 1 + (side == -1 ? 0 : _nodes[side].size);
        node.height = 1 + (side == -1 ? 0 : _nodes[side].height);
        int id = new_node(std::move(node));
        _vertex_node[vertex] = id;
        return id;
    }

    int weighted_split(const std::vector<int>& nodes, int l, int r) const {
        int total = 0;
        for (int i = l; i < r; i++) total += _nodes[nodes[i]].size;
        int left_sum = 0;
        for (int i = l; i + 1 < r; i++) {
            left_sum += _nodes[nodes[i]].size;
            if (2 * left_sum >= total) return i + 1;
        }
        return r - 1;
    }

    int build_rake(const std::vector<int>& nodes, int l, int r) {
        if (l == r) return -1;
        if (l + 1 == r) return nodes[l];
        int m = weighted_split(nodes, l, r);
        return new_rake(build_rake(nodes, l, m), build_rake(nodes, m, r));
    }

    int build_compress(const std::vector<int>& nodes, const std::vector<edge_type>& edges, int l, int r) {
        if (l + 1 == r) return nodes[l];
        int m = weighted_split(nodes, l, r);
        return new_compress(build_compress(nodes, edges, l, m), build_compress(nodes, edges, m, r), edges[m - 1]);
    }

    int build_vertex(int v) {
        std::vector<int> side_nodes;
        for (const auto& e : _children[v]) {
            if (e.to == _heavy[v]) continue;
            int child_path = build_path(e.to);
            side_nodes.push_back(new_add_edge(child_path, e));
        }
        return new_add_vertex(build_rake(side_nodes, 0, int(side_nodes.size())), v);
    }

    int build_path(int start) {
        std::vector<int> path_nodes;
        std::vector<edge_type> path_edges;
        for (int v = start; v != -1; v = _heavy[v]) {
            path_nodes.push_back(build_vertex(v));
            if (_heavy[v] != -1) path_edges.push_back(_heavy_edge[v]);
        }
        return build_compress(path_nodes, path_edges, 0, int(path_nodes.size()));
    }

    void recompute_up(int node) {
        while (node != -1) {
            recompute(node);
            node = _nodes[node].parent;
        }
    }

   public:
    StaticTopTree(const m1une::graph::Graph<T>& g, const std::vector<Vertex>& values, Point point_id,
                  Compress compress, Rake rake, AddEdge add_edge, AddVertex add_vertex, int root = 0)
        : _n(g.size()),
          _root(_n == 0 ? -1 : root),
          _root_node(-1),
          _point_id(std::move(point_id)),
          _compress(std::move(compress)),
          _rake(std::move(rake)),
          _add_edge(std::move(add_edge)),
          _add_vertex(std::move(add_vertex)),
          _values(values) {
        build(g, root);
    }

    void build(const m1une::graph::Graph<T>& g, int root = 0) {
        _n = g.size();
        _root = _n == 0 ? -1 : root;
        assert(int(_values.size()) == _n);
        _nodes.clear();
        _vertex_node.assign(_n, -1);
        _edge_node.assign(g.edge_count(), -1);
        _parent.assign(_n, -2);
        _subtree_size.assign(_n, 1);
        _heavy.assign(_n, -1);
        _heavy_edge.assign(_n, edge_type());
        _children.assign(_n, {});
        _root_node = -1;

        if (_n == 0) return;
        assert(0 <= root && root < _n);
        assert(int(g.edges().size()) == _n - 1);

        std::vector<int> order;
        order.reserve(_n);
        std::vector<int> stack = {root};
        _parent[root] = -1;
        while (!stack.empty()) {
            int v = stack.back();
            stack.pop_back();
            order.push_back(v);
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                if (_parent[e.to] != -2) continue;
                _parent[e.to] = v;
                _children[v].push_back(e);
                stack.push_back(e.to);
            }
        }
        assert(int(order.size()) == _n);

        for (int i = int(order.size()) - 1; i >= 0; i--) {
            int v = order[i];
            for (const auto& e : _children[v]) {
                _subtree_size[v] += _subtree_size[e.to];
                if (_heavy[v] == -1 || _subtree_size[_heavy[v]] < _subtree_size[e.to]) {
                    _heavy[v] = e.to;
                    _heavy_edge[v] = e;
                }
            }
        }

        _root_node = build_path(root);
    }

    int size() const {
        return _n;
    }

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

    int root() const {
        return _root;
    }

    int node_count() const {
        return int(_nodes.size());
    }

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

    const Vertex& get(int v) const {
        assert(0 <= v && v < _n);
        return _values[v];
    }

    const Vertex& operator[](int v) const {
        return get(v);
    }

    void set(int v, const Vertex& value) {
        assert(0 <= v && v < _n);
        assert(_vertex_node[v] != -1);
        _values[v] = value;
        recompute_up(_vertex_node[v]);
    }

    void set(int v, Vertex&& value) {
        assert(0 <= v && v < _n);
        assert(_vertex_node[v] != -1);
        _values[v] = std::move(value);
        recompute_up(_vertex_node[v]);
    }

    void set_edge_cost(int edge_id, T cost) {
        assert(0 <= edge_id && edge_id < int(_edge_node.size()));
        int node = _edge_node[edge_id];
        assert(node != -1);
        _nodes[node].edge.cost = cost;
        recompute_up(node);
    }

    const Path& all_prod() const {
        assert(_root_node != -1);
        return path_value(_root_node);
    }

    const Path& query() const {
        return all_prod();
    }
};

template <class T, class Vertex, class Point, class Compress, class Rake, class AddEdge, class AddVertex>
StaticTopTree(const m1une::graph::Graph<T>&, const std::vector<Vertex>&, Point, Compress, Rake, AddEdge,
              AddVertex, int)
    -> StaticTopTree<T, Vertex, std::invoke_result_t<AddVertex, Point, Vertex, int>, Point, Compress, Rake,
                     AddEdge, AddVertex>;

template <class T, class Vertex, class Point, class Compress, class Rake, class AddEdge, class AddVertex>
StaticTopTree(const m1une::graph::Graph<T>&, const std::vector<Vertex>&, Point, Compress, Rake, AddEdge, AddVertex)
    -> StaticTopTree<T, Vertex, std::invoke_result_t<AddVertex, Point, Vertex, int>, Point, Compress, Rake,
                     AddEdge, AddVertex>;

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/tree.hpp"



#line 9 "graph/tree/tree.hpp"


#line 1 "graph/tree/tree_hash.hpp"



#line 9 "graph/tree/tree_hash.hpp"

#line 11 "graph/tree/tree_hash.hpp"

namespace m1une {
namespace tree {

using TreeHashValue = std::array<std::uint64_t, 2>;

class TreeHasher {
   private:
    static constexpr std::uint64_t mod = (std::uint64_t(1) << 61) - 1;
    std::uint64_t _seed;

    static std::uint64_t splitmix64(std::uint64_t x) {
        x += 0x9e3779b97f4a7c15ULL;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
        return x ^ (x >> 31);
    }

    static std::uint64_t mul_mod(std::uint64_t a, std::uint64_t b) {
        __uint128_t product = static_cast<__uint128_t>(a) * b;
        std::uint64_t result = std::uint64_t(product & mod) + std::uint64_t(product >> 61);
        if (mod <= result) result -= mod;
        return result;
    }

    static std::uint64_t add_mod(std::uint64_t a, std::uint64_t b) {
        std::uint64_t result = a + b;
        if (mod <= result) result -= mod;
        return result;
    }

    TreeHashValue salt(int height) const {
        std::uint64_t x = static_cast<std::uint64_t>(height);
        std::uint64_t first = splitmix64(_seed ^ (x + 0x243f6a8885a308d3ULL));
        std::uint64_t second = splitmix64(_seed ^ (x + 0x13198a2e03707344ULL));
        return {first % (mod - 1) + 1, second % (mod - 1) + 1};
    }

    template <class T>
    static std::vector<int> tree_centers(const m1une::graph::Graph<T>& g) {
        int n = g.size();
        if (n == 0) return {};

        std::vector<int> degree(n, 0);
        std::vector<int> queue;
        queue.reserve(n);
        long long active_arcs = 0;
        for (int v = 0; v < n; v++) {
            for (const auto& e : g[v]) {
                if (!e.alive) continue;
                degree[v]++;
                active_arcs++;
            }
            if (degree[v] <= 1) queue.push_back(v);
        }
        assert(active_arcs == 2LL * (n - 1));

        std::vector<char> removed(n, false);
        int remaining = n;
        int head = 0;
        while (2 < remaining) {
            int layer_end = int(queue.size());
            assert(head < layer_end);
            remaining -= layer_end - head;
            while (head < layer_end) {
                int v = queue[head++];
                removed[v] = true;
                for (const auto& e : g[v]) {
                    if (!e.alive || removed[e.to]) continue;
                    if (--degree[e.to] == 1) queue.push_back(e.to);
                }
            }
        }

        std::vector<int> centers;
        for (int v = 0; v < n; v++) {
            if (!removed[v]) centers.push_back(v);
        }
        assert(1 <= int(centers.size()) && int(centers.size()) <= 2);
        return centers;
    }

   public:
    explicit TreeHasher(std::uint64_t seed = 0x6a09e667f3bcc909ULL) : _seed(seed) {}

    std::uint64_t seed() const {
        return _seed;
    }

    template <class T>
    std::vector<TreeHashValue> hash_subtrees(const m1une::graph::Graph<T>& g, int root = 0) const {
        int n = g.size();
        if (n == 0) return {};
        assert(0 <= root && root < n);

        std::vector<int> parent(n, -1);
        std::vector<int> order;
        order.reserve(n);
        parent[root] = root;
        order.push_back(root);
        long long active_arcs = 0;
        for (int v = 0; v < n; v++) {
            for (const auto& e : g[v]) active_arcs += e.alive;
        }
        assert(active_arcs == 2LL * (n - 1));

        for (int i = 0; i < int(order.size()); i++) {
            int v = order[i];
            for (const auto& e : g[v]) {
                if (!e.alive || parent[e.to] != -1) continue;
                parent[e.to] = v;
                order.push_back(e.to);
            }
        }
        assert(int(order.size()) == n);

        std::vector<int> height(n, 0);
        std::vector<TreeHashValue> result(n, TreeHashValue{1, 1});
        for (int i = n - 1; i >= 0; i--) {
            int v = order[i];
            for (const auto& e : g[v]) {
                if (!e.alive || parent[e.to] != v) continue;
                height[v] = std::max(height[v], height[e.to] + 1);
            }
            TreeHashValue random = salt(height[v]);
            for (const auto& e : g[v]) {
                if (!e.alive || parent[e.to] != v) continue;
                result[v][0] = mul_mod(result[v][0], add_mod(result[e.to][0], random[0]));
                result[v][1] = mul_mod(result[v][1], add_mod(result[e.to][1], random[1]));
            }
        }
        return result;
    }

    template <class T>
    TreeHashValue hash_rooted(const m1une::graph::Graph<T>& g, int root = 0) const {
        if (g.empty()) return {0, 0};
        return hash_subtrees(g, root)[root];
    }

    template <class T>
    std::vector<TreeHashValue> hash_unrooted(const m1une::graph::Graph<T>& g) const {
        std::vector<int> centers = tree_centers(g);
        std::vector<TreeHashValue> result;
        result.reserve(centers.size());
        for (int center : centers) result.push_back(hash_rooted(g, center));
        std::sort(result.begin(), result.end());
        return result;
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/virtual_tree.hpp"



#line 8 "graph/tree/virtual_tree.hpp"

#line 11 "graph/tree/virtual_tree.hpp"

namespace m1une {
namespace tree {

template <class T = int>
struct VirtualTreeResult {
    std::vector<int> vertex;
    std::vector<int> parent;
    std::vector<int> parent_edge_count;
    std::vector<T> parent_cost;
    std::vector<std::vector<int>> children;
    std::vector<bool> is_key;

    int size() const {
        return int(vertex.size());
    }

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

    int edge_count() const {
        return vertex.empty() ? 0 : int(vertex.size()) - 1;
    }

    int root() const {
        return vertex.empty() ? -1 : 0;
    }

    int root_vertex() const {
        return vertex.empty() ? -1 : vertex[0];
    }
};

template <class T = int>
struct VirtualTree {
    using cost_type = T;
    using result_type = VirtualTreeResult<T>;

   private:
    SparseTableLca<T> _lca;
    std::vector<int> _key;
    std::vector<int> _vertices;
    std::vector<int> _stack;

   public:
    VirtualTree() = default;

    explicit VirtualTree(const m1une::graph::Graph<T>& graph, int root = 0) : _lca(graph, root) {}

    void build_lca(const m1une::graph::Graph<T>& graph, int root = 0) {
        _lca.build(graph, root);
    }

    int original_size() const {
        return _lca.size();
    }

    const SparseTableLca<T>& lca_data() const {
        return _lca;
    }

    result_type build(std::vector<int> key_vertices) {
        result_type result;
        if (key_vertices.empty()) return result;

        auto by_tin = [&](int u, int v) { return _lca.tin[u] < _lca.tin[v]; };
        for (int v : key_vertices) {
            assert(0 <= v && v < _lca.size());
            assert(_lca.tin[v] != -1);
        }
        std::sort(key_vertices.begin(), key_vertices.end(), by_tin);
        key_vertices.erase(std::unique(key_vertices.begin(), key_vertices.end()), key_vertices.end());

        _key = key_vertices;
        _vertices = key_vertices;
        _vertices.reserve(2 * _key.size());
        for (int i = 1; i < int(_key.size()); i++) {
            _vertices.push_back(_lca.lca(_key[i - 1], _key[i]));
        }
        std::sort(_vertices.begin(), _vertices.end(), by_tin);
        _vertices.erase(std::unique(_vertices.begin(), _vertices.end()), _vertices.end());

        int n = int(_vertices.size());
        result.vertex = _vertices;
        result.parent.assign(n, -1);
        result.parent_edge_count.assign(n, 0);
        result.parent_cost.assign(n, T(0));
        result.children.assign(n, {});
        result.is_key.assign(n, false);

        int key_index = 0;
        for (int i = 0; i < n; i++) {
            while (key_index < int(_key.size()) && _lca.tin[_key[key_index]] < _lca.tin[_vertices[i]]) {
                key_index++;
            }
            if (key_index < int(_key.size()) && _key[key_index] == _vertices[i]) result.is_key[i] = true;
        }

        _stack.clear();
        _stack.reserve(n);
        for (int i = 0; i < n; i++) {
            while (!_stack.empty() && !_lca.is_ancestor(_vertices[_stack.back()], _vertices[i])) {
                _stack.pop_back();
            }
            if (!_stack.empty()) {
                int p = _stack.back();
                result.parent[i] = p;
                result.parent_edge_count[i] = _lca.depth[_vertices[i]] - _lca.depth[_vertices[p]];
                result.parent_cost[i] = _lca.dist[_vertices[i]] - _lca.dist[_vertices[p]];
                result.children[p].push_back(i);
            }
            _stack.push_back(i);
        }
        return result;
    }
};

}  // namespace tree
}  // namespace m1une


#line 1 "graph/tree/zero_one_on_tree.hpp"



#line 5 "graph/tree/zero_one_on_tree.hpp"
#include <set>
#line 7 "graph/tree/zero_one_on_tree.hpp"

#line 9 "graph/tree/zero_one_on_tree.hpp"

namespace m1une {
namespace tree {

inline long long zero_one_on_tree(const std::vector<int>& parent,
                                  const std::vector<int>& value) {
    const int n = int(parent.size());
    assert(int(value.size()) == n);
    if (n == 0) return 0;

    int root = -1;
    std::vector<std::vector<int>> children(n);
    for (int v = 0; v < n; v++) {
        assert(value[v] == 0 || value[v] == 1);
        if (parent[v] == -1) {
            assert(root == -1);
            root = v;
        } else {
            assert(0 <= parent[v] && parent[v] < n && parent[v] != v);
            children[parent[v]].push_back(v);
        }
    }
    assert(root != -1);

    std::vector<int> stack(1, root);
    std::vector<char> visited(n, false);
    visited[root] = true;
    int visited_count = 0;
    while (!stack.empty()) {
        const int v = stack.back();
        stack.pop_back();
        visited_count++;
        for (int child : children[v]) {
            assert(!visited[child]);
            visited[child] = true;
            stack.push_back(child);
        }
    }
    assert(visited_count == n);

    struct Component {
        long long zeros;
        long long ones;
        int vertex;
    };
    struct Compare {
        bool operator()(const Component& lhs, const Component& rhs) const {
            const long long lhs_product = lhs.zeros * rhs.ones;
            const long long rhs_product = rhs.zeros * lhs.ones;
            if (lhs_product != rhs_product) return lhs_product < rhs_product;
            return lhs.vertex < rhs.vertex;
        }
    };

    std::vector<long long> zeros(n), ones(n);
    std::vector<int> dsu(n);
    std::set<Component, Compare> components;
    for (int v = 0; v < n; v++) {
        zeros[v] = value[v] == 0;
        ones[v] = value[v] == 1;
        dsu[v] = v;
        if (v != root) components.insert(Component{zeros[v], ones[v], v});
    }

    auto leader = [&](int v) {
        int result = v;
        while (dsu[result] != result) result = dsu[result];
        while (dsu[v] != v) {
            const int next = dsu[v];
            dsu[v] = result;
            v = next;
        }
        return result;
    };

    long long answer = 0;
    while (!components.empty()) {
        auto it = components.end();
        --it;
        const Component child = *it;
        components.erase(it);

        const int p = leader(parent[child.vertex]);
        if (p != root) {
            const int erased = int(components.erase(Component{zeros[p], ones[p], p}));
            assert(erased == 1);
        }

        answer += ones[p] * zeros[child.vertex];
        zeros[p] += zeros[child.vertex];
        ones[p] += ones[child.vertex];
        dsu[child.vertex] = p;

        if (p != root) components.insert(Component{zeros[p], ones[p], p});
    }
    return answer;
}

template <class T>
long long zero_one_on_tree(const m1une::graph::Graph<T>& graph,
                           const std::vector<int>& value, int root = 0) {
    const int n = graph.size();
    assert(int(value.size()) == n);
    if (n == 0) return 0;
    assert(0 <= root && root < n);
    assert(int(graph.edges().size()) == n - 1);

    RootedTree<T> rooted_tree(graph, root);
    assert(int(rooted_tree.order.size()) == n);
    return zero_one_on_tree(rooted_tree.parent, value);
}

}  // namespace tree
}  // namespace m1une


#line 23 "graph/tree/all.hpp"
Back to top page