Tree Path Product

tree/path_product.hpp · GitHub で見る

共通ページ

固定した木の頂点値を更新し、2 頂点を結ぶ経路の積を順序付きで求めます。結合則と単位元を持つ演算を指定します。可換性は不要です。

#include "lib/tree/path_product.hpp"
#include <cstdint>

std::int64_t op(std::int64_t a, std::int64_t b) { return a + b; }
std::int64_t e() { return 0; }

HLD tree(4);
tree.add_edge(0, 1);
tree.add_edge(1, 2);
tree.add_edge(1, 3);
tree.build();
TreePathProduct<std::int64_t, op, e> values(tree, {1, 2, 3, 4});
assert(values.prod(2, 3) == 9);
values.set(1, 10);
assert(values.prod(2, 3) == 17);

操作	内容・返り値	計算量
`TreePathProduct<S,op,e,reverse = nullptr>()`	空の木を作る	\(O(1)\)
`TreePathProduct(tree,values)`	構築済みの `HLD` と、元の頂点番号順の `vector<S>` から構築	\(O(n)\)
`build(tree,values)`	木と値を置き換えて再構築	\(O(n)\)
`size()` / `empty()`	頂点数 `int`／空か `bool`	\(O(1)\)
`get(v)`	頂点 `v` の値 `const S&`	\(O(1)\)
`set(v,value)`	頂点 `v` の値を変更	\(O(\log n)\)
`prod(u,v)`	頂点 `u` から `v` へ向かう順の経路積 `S`。両端を含む	\(O(\log n)\)
`apply(u,v,value,mapping)`	経路積を値 `value` に作用させ、同じ型で返す	\(O(\log n)\)

頂点番号は [0,n) です。values.size() は頂点数と一致させます。空の木も構築できますが、頂点を指定する操作は空でない木に対して呼び出してください。構築後に元の HLD を変更・破棄しても影響しません。木の辺を変更する場合は再構築します。get の参照は再構築・代入・破棄で無効になります。

S はコピー・代入可能な型とし、op(a,b) は S の結合演算、e() は単位元を返す関数です。演算は頂点列の左から右へ適用する順序で定義します。

apply の mapping(segment,value) は、積 segment を値へ作用させる関数です。mapping(e(),x) == x および mapping(op(a,b),x) == mapping(b,mapping(a,x)) を満たしてください。例えば、S が一次関数 a*x+b を表すなら、mapping はその一次関数の評価を行います。

第 4 テンプレート引数 reverse は省略できます。指定する場合は、集約値 S から頂点列を逆順にした集約値を返す関数とし、reverse(op(a,b)) == op(reverse(b),reverse(a))、reverse(e()) == e()、各頂点値 v について reverse(v) == v を満たしてください。順方向・逆方向の両方を持つ型などで指定できます。

頂点数を \(n\) として、構築は \(O(n)\) 時間・空間、頂点値更新・経路積・経路への作用は \(O(\log n)\) 時間です。値のコピーと演算・作用は定数時間とします。

依存ライブラリ

HLD (Heavy-Light Decomposition)

検証コード

test/hld_path_composite.test.cpp

コード

#line 1 "tree/path_product.hpp"

#line 1 "tree/hld.hpp"

#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <limits>
#include <utility>
#include <vector>

class HLD {
	int shape_ = 0, center_ = 0;
	struct Node { int position, head, parent, offset; };
	std::vector<Node> nodes;
	void check([[maybe_unused]] int v) const { assert(0 <= v && v < n && int(in.size()) == n); }

void initialize() {
		depth.assign(n, 0);
		size.assign(n, 1);
		heavy.assign(n, -1);
		head.resize(n), in.resize(n), out.resize(n), order.resize(n);
		nodes.resize(n);
		shape_ = center_ = 0;
	}
	void finish(int root) {
		in[root] = 0, out[root] = n, head[root] = root;
		for (int i = n - 2; i >= 0; --i) {
			int v = order[i], p = parent[v];
			depth[v] = depth[p] + 1;
			if (heavy[p] == v) in[v] = in[p] + 1, head[v] = head[p];
			else in[v] = out[p] -= size[v], head[v] = v;
			out[v] = in[v] + size[v];
		}
		for (int v = 0; v < n; ++v) out[v] = in[v] + size[v], order[in[v]] = v;
		parent[root] = -1;
		for (int v = 0; v < n; ++v) nodes[v] = {in[v], in[head[v]], parent[head[v]], depth[v] - in[v]};
	}

public:
	struct Segment {
		int l, r;
		bool reversed;
	};
	int n;
	std::vector<std::pair<int, int>> edges;
	std::vector<int> parent, depth, size, heavy, head, in, out, order;
	explicit HLD(int n = 0) : n(n) {
		assert(0 <= n && n <= std::numeric_limits<int>::max() / 2);
		edges.reserve(std::max(0, n - 1));
	}
	int add_edge(int u, int v) {
		assert(0 <= u && u < n && 0 <= v && v < n);
		edges.emplace_back(u, v);
		return int(edges.size()) - 1;
	}
	void build(int root = 0) {
		assert(n >= 0 && edges.size() == std::size_t(std::max(0, n - 1)));
		assert(!n || (0 <= root && root < n));
		parent.assign(n, -1);
		initialize();
		if (!n) return;
		std::vector<int> degree(n), neighbors(n);
		for (auto [u, v] : edges) {
			++degree[u], ++degree[v];
			neighbors[u] ^= v, neighbors[v] ^= u;
		}
		center_ = int(std::max_element(degree.begin(), degree.end()) - degree.begin());
		shape_ = degree[center_] <= 2 ? 1 : degree[center_] == n - 1 ? 2 : 0;
		degree[root] += 2;
		int count = 0;
		for (int v = 0; v < n; ++v) if (degree[v] == 1) order[count++] = v;
		for (int i = 0; i < count; ++i) {
			int v = order[i], p = neighbors[v];
			parent[v] = p;
			size[p] += size[v];
			if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
			neighbors[p] ^= v;
			if (--degree[p] == 1) order[count++] = p;
		}
		assert(count == n - 1);
		finish(root);
	}
	void build_from_parents(std::vector<int> parents) {
		assert(parents.size() <= std::size_t(std::numeric_limits<int>::max() / 2));
		n = static_cast<int>(parents.size());
		parent = std::move(parents);
		initialize();
		edges.clear();
		edges.reserve(std::max(0, n - 1));
		if (!n) return;
		std::fill(out.begin(), out.end(), 0);
		int root = -1;
		bool ordered = true;
		for (int v = 0; v < n; ++v) {
			int p = parent[v];
			assert(-1 <= p && p < n && p != v);
			if (p == -1) {
				assert(root == -1);
				root = v;
			} else {
				++out[p];
				edges.emplace_back(p, v);
				ordered &= p < v;
			}
		}
		assert(root != -1);
		int maximum = 0;
		for (int v = 0; v < n; ++v) {
			int degree = out[v] + (v != root);
			if (degree > maximum) maximum = degree, center_ = v;
		}
		shape_ = maximum <= 2 ? 1 : maximum == n - 1 ? 2 : 0;
		if (ordered) {
			for (int v = n - 1; v > 0; --v) {
				int p = parent[v];
				size[p] += size[v];
				if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
				order[n - 1 - v] = v;
			}
		} else {
			int count = 0;
			for (int v = 0; v < n; ++v) if (!out[v] && v != root) order[count++] = v;
			for (int i = 0; i < count; ++i) {
				int v = order[i], p = parent[v];
				size[p] += size[v];
				if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
				if (!--out[p] && p != root) order[count++] = p;
			}
			assert(count == n - 1);
		}
		finish(root);
	}
	int lca(int u, int v) const {
		check(u), check(v);
		if (shape_ == 1) return head[u] != head[v] ? order[0] : in[u] < in[v] ? u : v;
		if (shape_ == 2) return u == v ? u : u == order[0] || v == order[0] ? order[0] : center_;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) u = nodes[u].parent;
			else v = nodes[v].parent;
		}
		return nodes[u].position < nodes[v].position ? u : v;
	}
	int dist(int u, int v) const {
		int w = lca(u, v);
		return depth[u] + depth[v] - 2 * depth[w];
	}
	bool is_ancestor(int u, int v) const {
		check(u), check(v);
		return in[u] <= in[v] && in[v] < out[u];
	}
	int ancestor(int v, int k) const {
		check(v);
		if (shape_) {
			if (k < 0 || k > depth[v]) return -1;
			if (!k) return v;
			if (k == depth[v]) return order[0];
			return shape_ == 2 ? parent[v] : order[in[v] - k];
		}
		if (k < 0) return -1;
		int target = nodes[v].position + nodes[v].offset - k;
		if (target < 0) return -1;
		while (nodes[v].head + nodes[v].offset > target) v = nodes[v].parent;
		return order[target - nodes[v].offset];
	}
	int jump(int u, int v, int k) const {
		check(u), check(v);
		if (k < 0) return -1;
		if (!k) return u;
		if (shape_ == 1) {
			if (head[u] == head[v]) {
				int delta = in[v] - in[u];
				if (k > (delta < 0 ? -delta : delta)) return -1;
				return order[in[u] + (delta < 0 ? -k : k)];
			}
			int a = depth[u], b = depth[v];
			if (k > a + b) return -1;
			if (k == a) return order[0];
			return k < a ? order[in[u] - k] : order[in[v] - (a + b - k)];
		}
		if (shape_ == 2) {
			int d = (u != center_) + (v != center_);
			if (u == v || k > d) return -1;
			return k == d ? v : center_;
		}
		int from = u, to = v, target = nodes[u].position + nodes[u].offset - k;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) {
				if (nodes[u].head + nodes[u].offset <= target) return order[target - nodes[u].offset];
				u = nodes[u].parent;
			} else v = nodes[v].parent;
		}
		int w = nodes[u].position < nodes[v].position ? u : v;
		int d = nodes[w].position + nodes[w].offset;
		int length = nodes[from].position + nodes[from].offset + nodes[to].position + nodes[to].offset - 2 * d;
		if (k > length) return -1;
		if (target >= d) return order[target - nodes[u].offset];
		return ancestor(to, length - k);
	}
	std::pair<int, int> subtree(int v, bool edge = false) const {
		check(v);
		return {in[v] + edge, out[v]};
	}
	int edge_vertex(int id) const {
		assert(0 <= id && id < int(edges.size()));
		auto [u, v] = edges[id];
		check(u);
		return parent[u] == v ? u : v;
	}
	template <class F>
	void for_each_path(int u, int v, F &&f, bool edge = false) const {
		check(u), check(v);
		if (shape_ == 1) {
			int x = in[u], y = in[v];
			if (head[u] == head[v]) {
				int l = std::min(x, y) + edge, r = std::max(x, y) + 1;
				if (l < r) f(l, r, x > y);
			} else if (head[u] == order[0]) {
				if (edge < x + 1) f(int(edge), x + 1, true);
				f(in[head[v]], y + 1, false);
			} else {
				f(in[head[u]], x + 1, true);
				if (edge < y + 1) f(int(edge), y + 1, false);
			}
			return;
		}
		std::array<std::array<int, 2>, 32> right;
		int count = 0;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) {
				f(nodes[u].head, nodes[u].position + 1, true);
				u = nodes[u].parent;
			} else {
				assert(count < int(right.size()));
				right[count++] = {nodes[v].head, nodes[v].position + 1};
				v = nodes[v].parent;
			}
		}
		bool reversed = nodes[u].position > nodes[v].position;
		int l = std::min(nodes[u].position, nodes[v].position) + edge, r = std::max(nodes[u].position, nodes[v].position) + 1;
		if (l < r) f(l, r, reversed);
		while (count) {
			auto [a, b] = right[--count];
			f(a, b, false);
		}
	}
	std::vector<Segment> path(int u, int v, bool edge = false) const {
		std::vector<Segment> result;
		for_each_path(u, v, [&](int l, int r, bool reversed) { result.push_back({l, r, reversed}); }, edge);
		return result;
	}
};

#line 5 "tree/path_product.hpp"
#include <bit>
#include <cstdint>
#include <type_traits>

template<class S, auto op, auto e, auto reverse = nullptr> class TreePathProduct {
	static constexpr bool has_reverse = !std::is_same_v<decltype(reverse), std::nullptr_t>;
	struct Empty {};
	struct Back { S prod, prefix; };
	using BackStorage = std::conditional_t<has_reverse, Empty, Back>;
	struct alignas(64) Node {
		int left = 0, right = 0, parent = 0, head = 0, up = 0;
		unsigned char priority = 0;
		S value, prod, prefix;
		[[no_unique_address]] BackStorage back;
		static BackStorage make_back(const S& value) {
			if constexpr (has_reverse) return {};
			else return {value, value};
		}
		explicit Node(const S& v) : value(v), prod(v), prefix(v), back(make_back(v)) {}
	};
	std::vector<Node> nodes_;
	std::vector<int> index_;
	std::vector<S> dense_;
	std::vector<int> dense_heads_, dense_up_;
	unsigned base_ = 1;
	bool flat_ = false;
	int index(int v) const {
		assert(0 <= v && v < size());
		return index_[v];
	}
	void pull(int v) {
		auto& node = nodes_[v];
		node.prefix = node.left ? op(nodes_[node.left].prod, node.value) : node.value;
		node.prod = node.right ? op(node.prefix, nodes_[node.right].prod) : node.prefix;
		if constexpr (!has_reverse) {
			node.back.prefix = node.left ? op(node.value, nodes_[node.left].back.prod) : node.value;
			node.back.prod = node.right ? op(nodes_[node.right].back.prod, node.back.prefix) : node.back.prefix;
		}
	}
	S back_prod(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prod);
		else return nodes_[v].back.prod;
	}
	S back_prefix(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prefix);
		else return nodes_[v].back.prefix;
	}
	template<bool Reverse> S dense_interval(int l, int r) const {
		unsigned a = base_ + unsigned(l - 1), b = base_ + unsigned(r);
		S left = e(), right = e();
		while (a < b) {
			if (a & 1) left = op(left, dense_[a++]);
			if (b & 1) right = op(dense_[--b], right);
			a >>= 1, b >>= 1;
		}
		S result = op(left, right);
		if constexpr (Reverse && has_reverse) return reverse(result);
		else return result;
	}
	template<bool Reverse> S prefix(int v) const {
		S result = e();
		int last = v;
		for (; v; v = nodes_[v].parent) if (v <= last) {
			if constexpr (Reverse) result = op(result, back_prefix(v));
			else result = op(nodes_[v].prefix, result);
		}
		return result;
	}
	template<bool Reverse> S interval(int l, int r) const {
		int a = l, b = r;
		S left = e(), right = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) left = op(op(back_prod(node.right), node.value), left);
					else left = op(left, op(node.value, nodes_[node.right].prod));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) right = op(right, back_prefix(b));
					else right = op(nodes_[b].prefix, right);
				}
				b = nodes_[b].parent;
			}
		}
		if constexpr (Reverse) return op(right, op(nodes_[a].value, left));
		else return op(left, op(nodes_[a].value, right));
	}

template<bool Reverse, class X, class Apply> X apply_interval(int l, int r, X value, Apply& mapping) const {
		int a = l, b = r;
		S outside = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) outside = op(op(back_prod(node.right), node.value), outside);
					else value = mapping(nodes_[node.right].prod, mapping(node.value, std::move(value)));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) value = mapping(back_prefix(b), std::move(value));
					else outside = op(nodes_[b].prefix, outside);
				}
				b = nodes_[b].parent;
			}
		}
		return mapping(outside, mapping(nodes_[a].value, std::move(value)));
	}

template<class X, class Apply> X dense_apply(int u, int v, X value, Apply& mapping) const {
		std::array<std::pair<int, int>, 4> right;
		int count = 0;
		while (dense_heads_[u] != dense_heads_[v]) {
			if (dense_heads_[u] > dense_heads_[v]) {
				value = mapping(dense_interval<true>(dense_heads_[u], u), std::move(value));
				u = dense_up_[u];
			} else {
				assert(count < int(right.size()));
				right[count++] = {dense_heads_[v], v};
				v = dense_up_[v];
			}
		}
		if (u <= v) value = mapping(dense_interval<false>(u, v), std::move(value));
		else value = mapping(dense_interval<true>(v, u), std::move(value));
		while (count) {
			auto [l, r] = right[--count];
			value = mapping(dense_interval<false>(l, r), std::move(value));
		}
		return value;
	}

public:
	TreePathProduct() = default;
	TreePathProduct(const HLD& tree, const std::vector<S>& values) { build(tree, values); }
	void build(const HLD& tree, const std::vector<S>& values) {
		assert(values.size() == std::size_t(tree.n) && tree.order.size() == values.size());
		index_.resize(values.size());
		nodes_.clear();
		for (int v : tree.order) index_[v] = tree.in[v] + 1;
		flat_ = false;
		dense_.clear();
		if constexpr (has_reverse) {
			std::vector<unsigned char> levels(values.size());
			unsigned maximum = 0;
			for (int v : tree.order) {
				if (tree.parent[v] != -1) levels[v] = levels[tree.parent[v]] + (tree.head[v] == v);
				maximum = std::max(maximum, unsigned(levels[v]));
			}
			flat_ = maximum <= 4;
		}
		if (flat_) {
			base_ = std::bit_ceil(unsigned(size()) ? unsigned(size()) : 1U);
			dense_.assign(std::size_t(base_) * 2, e());
			dense_heads_.assign(values.size() + 1, 0), dense_up_.assign(values.size() + 1, 0);
			for (int v : tree.order) {
				int id = index_[v], h = tree.head[v];
				dense_heads_[id] = index_[h];
				dense_up_[id] = tree.parent[h] == -1 ? 0 : index_[tree.parent[h]];
				dense_[base_ + unsigned(id - 1)] = values[v];
			}
			for (unsigned i = base_ - 1; i; --i) dense_[i] = op(dense_[i * 2], dense_[i * 2 + 1]);
			return;
		}
		dense_heads_.clear(), dense_up_.clear();
		nodes_.reserve(values.size() + 1);
		nodes_.emplace_back(e());
		for (int v : tree.order) nodes_.emplace_back(values[v]);
		std::vector<int> stack, order;
		order.reserve(values.size());
		for (int v : tree.order) {
			if (tree.head[v] != v) continue;
			stack.clear();
			unsigned position = 0;
			int head = index_[v], up = tree.parent[v] == -1 ? 0 : index_[tree.parent[v]];
			for (int u = v; u != -1; u = tree.heavy[u]) {
				int id = index_[u], weight = tree.size[u] - (tree.heavy[u] == -1 ? 0 : tree.size[tree.heavy[u]]);
				auto& node = nodes_[id];
				node.head = head, node.up = up;
				node.priority = std::bit_width(position ^ (position + unsigned(weight)));
				position += unsigned(weight);
				int last = 0;
				while (!stack.empty() && nodes_[stack.back()].priority < node.priority) {
					last = stack.back();
					stack.pop_back();
				}
				node.left = last;
				if (last) nodes_[last].parent = id;
				if (!stack.empty()) nodes_[stack.back()].right = id, node.parent = stack.back();
				stack.push_back(id);
			}
			order.push_back(stack[0]);
		}
		for (std::size_t i = 0; i < order.size(); ++i) {
			const auto& node = nodes_[order[i]];
			if (node.left) order.push_back(node.left);
			if (node.right) order.push_back(node.right);
		}
		for (auto it = order.rbegin(); it != order.rend(); ++it) pull(*it);
	}
	int size() const { return int(index_.size()); }
	bool empty() const { return index_.empty(); }
	const S& get(int v) const {
		int id = index(v);
		return flat_ ? dense_[base_ + unsigned(id - 1)] : nodes_[id].value;
	}
	void set(int v, const S& value) {
		int id = index(v);
		if (flat_) {
			unsigned k = base_ + unsigned(id - 1);
			dense_[k] = value;
			for (k >>= 1; k; k >>= 1) dense_[k] = op(dense_[k * 2], dense_[k * 2 + 1]);
			return;
		}
		nodes_[id].value = value;
		for (; id; id = nodes_[id].parent) pull(id);
	}
	S prod(int u, int v) const {
		u = index(u), v = index(v);
		if (flat_) {
			auto append = [](const S& part, S value) { return op(std::move(value), part); };
			return dense_apply(u, v, e(), append);
		}
		S left = e(), right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				left = op(left, prefix<true>(u));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		S middle = u <= v ? interval<false>(u, v) : interval<true>(v, u);
		return op(left, op(middle, right));
	}
	template<class X, class Apply> X apply(int u, int v, X value, Apply mapping) const {
		u = index(u), v = index(v);
		if (flat_) return dense_apply(u, v, std::move(value), mapping);
		S right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				for (int id = u; id; id = nodes_[id].parent)
					if (id <= u) value = mapping(back_prefix(id), std::move(value));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		if (u <= v) value = apply_interval<false>(u, v, std::move(value), mapping);
		else value = apply_interval<true>(v, u, std::move(value), mapping);
		return mapping(right, std::move(value));
	}

};

#ifndef TREE_PATH_PRODUCT_HPP
#define TREE_PATH_PRODUCT_HPP

#include "hld.hpp"
#include <bit>
#include <cstdint>
#include <type_traits>

template<class S, auto op, auto e, auto reverse = nullptr> class TreePathProduct {
	static constexpr bool has_reverse = !std::is_same_v<decltype(reverse), std::nullptr_t>;
	struct Empty {};
	struct Back { S prod, prefix; };
	using BackStorage = std::conditional_t<has_reverse, Empty, Back>;
	struct alignas(64) Node {
		int left = 0, right = 0, parent = 0, head = 0, up = 0;
		unsigned char priority = 0;
		S value, prod, prefix;
		[[no_unique_address]] BackStorage back;
		static BackStorage make_back(const S& value) {
			if constexpr (has_reverse) return {};
			else return {value, value};
		}
		explicit Node(const S& v) : value(v), prod(v), prefix(v), back(make_back(v)) {}
	};
	std::vector<Node> nodes_;
	std::vector<int> index_;
	std::vector<S> dense_;
	std::vector<int> dense_heads_, dense_up_;
	unsigned base_ = 1;
	bool flat_ = false;
	int index(int v) const {
		assert(0 <= v && v < size());
		return index_[v];
	}
	void pull(int v) {
		auto& node = nodes_[v];
		node.prefix = node.left ? op(nodes_[node.left].prod, node.value) : node.value;
		node.prod = node.right ? op(node.prefix, nodes_[node.right].prod) : node.prefix;
		if constexpr (!has_reverse) {
			node.back.prefix = node.left ? op(node.value, nodes_[node.left].back.prod) : node.value;
			node.back.prod = node.right ? op(nodes_[node.right].back.prod, node.back.prefix) : node.back.prefix;
		}
	}
	S back_prod(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prod);
		else return nodes_[v].back.prod;
	}
	S back_prefix(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prefix);
		else return nodes_[v].back.prefix;
	}
	template<bool Reverse> S dense_interval(int l, int r) const {
		unsigned a = base_ + unsigned(l - 1), b = base_ + unsigned(r);
		S left = e(), right = e();
		while (a < b) {
			if (a & 1) left = op(left, dense_[a++]);
			if (b & 1) right = op(dense_[--b], right);
			a >>= 1, b >>= 1;
		}
		S result = op(left, right);
		if constexpr (Reverse && has_reverse) return reverse(result);
		else return result;
	}
	template<bool Reverse> S prefix(int v) const {
		S result = e();
		int last = v;
		for (; v; v = nodes_[v].parent) if (v <= last) {
			if constexpr (Reverse) result = op(result, back_prefix(v));
			else result = op(nodes_[v].prefix, result);
		}
		return result;
	}
	template<bool Reverse> S interval(int l, int r) const {
		int a = l, b = r;
		S left = e(), right = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) left = op(op(back_prod(node.right), node.value), left);
					else left = op(left, op(node.value, nodes_[node.right].prod));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) right = op(right, back_prefix(b));
					else right = op(nodes_[b].prefix, right);
				}
				b = nodes_[b].parent;
			}
		}
		if constexpr (Reverse) return op(right, op(nodes_[a].value, left));
		else return op(left, op(nodes_[a].value, right));
	}

	template<bool Reverse, class X, class Apply> X apply_interval(int l, int r, X value, Apply& mapping) const {
		int a = l, b = r;
		S outside = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) outside = op(op(back_prod(node.right), node.value), outside);
					else value = mapping(nodes_[node.right].prod, mapping(node.value, std::move(value)));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) value = mapping(back_prefix(b), std::move(value));
					else outside = op(nodes_[b].prefix, outside);
				}
				b = nodes_[b].parent;
			}
		}
		return mapping(outside, mapping(nodes_[a].value, std::move(value)));
	}

	template<class X, class Apply> X dense_apply(int u, int v, X value, Apply& mapping) const {
		std::array<std::pair<int, int>, 4> right;
		int count = 0;
		while (dense_heads_[u] != dense_heads_[v]) {
			if (dense_heads_[u] > dense_heads_[v]) {
				value = mapping(dense_interval<true>(dense_heads_[u], u), std::move(value));
				u = dense_up_[u];
			} else {
				assert(count < int(right.size()));
				right[count++] = {dense_heads_[v], v};
				v = dense_up_[v];
			}
		}
		if (u <= v) value = mapping(dense_interval<false>(u, v), std::move(value));
		else value = mapping(dense_interval<true>(v, u), std::move(value));
		while (count) {
			auto [l, r] = right[--count];
			value = mapping(dense_interval<false>(l, r), std::move(value));
		}
		return value;
	}

  public:
	TreePathProduct() = default;
	TreePathProduct(const HLD& tree, const std::vector<S>& values) { build(tree, values); }
	void build(const HLD& tree, const std::vector<S>& values) {
		assert(values.size() == std::size_t(tree.n) && tree.order.size() == values.size());
		index_.resize(values.size());
		nodes_.clear();
		for (int v : tree.order) index_[v] = tree.in[v] + 1;
		flat_ = false;
		dense_.clear();
		if constexpr (has_reverse) {
			std::vector<unsigned char> levels(values.size());
			unsigned maximum = 0;
			for (int v : tree.order) {
				if (tree.parent[v] != -1) levels[v] = levels[tree.parent[v]] + (tree.head[v] == v);
				maximum = std::max(maximum, unsigned(levels[v]));
			}
			flat_ = maximum <= 4;
		}
		if (flat_) {
			base_ = std::bit_ceil(unsigned(size()) ? unsigned(size()) : 1U);
			dense_.assign(std::size_t(base_) * 2, e());
			dense_heads_.assign(values.size() + 1, 0), dense_up_.assign(values.size() + 1, 0);
			for (int v : tree.order) {
				int id = index_[v], h = tree.head[v];
				dense_heads_[id] = index_[h];
				dense_up_[id] = tree.parent[h] == -1 ? 0 : index_[tree.parent[h]];
				dense_[base_ + unsigned(id - 1)] = values[v];
			}
			for (unsigned i = base_ - 1; i; --i) dense_[i] = op(dense_[i * 2], dense_[i * 2 + 1]);
			return;
		}
		dense_heads_.clear(), dense_up_.clear();
		nodes_.reserve(values.size() + 1);
		nodes_.emplace_back(e());
		for (int v : tree.order) nodes_.emplace_back(values[v]);
		std::vector<int> stack, order;
		order.reserve(values.size());
		for (int v : tree.order) {
			if (tree.head[v] != v) continue;
			stack.clear();
			unsigned position = 0;
			int head = index_[v], up = tree.parent[v] == -1 ? 0 : index_[tree.parent[v]];
			for (int u = v; u != -1; u = tree.heavy[u]) {
				int id = index_[u], weight = tree.size[u] - (tree.heavy[u] == -1 ? 0 : tree.size[tree.heavy[u]]);
				auto& node = nodes_[id];
				node.head = head, node.up = up;
				node.priority = std::bit_width(position ^ (position + unsigned(weight)));
				position += unsigned(weight);
				int last = 0;
				while (!stack.empty() && nodes_[stack.back()].priority < node.priority) {
					last = stack.back();
					stack.pop_back();
				}
				node.left = last;
				if (last) nodes_[last].parent = id;
				if (!stack.empty()) nodes_[stack.back()].right = id, node.parent = stack.back();
				stack.push_back(id);
			}
			order.push_back(stack[0]);
		}
		for (std::size_t i = 0; i < order.size(); ++i) {
			const auto& node = nodes_[order[i]];
			if (node.left) order.push_back(node.left);
			if (node.right) order.push_back(node.right);
		}
		for (auto it = order.rbegin(); it != order.rend(); ++it) pull(*it);
	}
	int size() const { return int(index_.size()); }
	bool empty() const { return index_.empty(); }
	const S& get(int v) const {
		int id = index(v);
		return flat_ ? dense_[base_ + unsigned(id - 1)] : nodes_[id].value;
	}
	void set(int v, const S& value) {
		int id = index(v);
		if (flat_) {
			unsigned k = base_ + unsigned(id - 1);
			dense_[k] = value;
			for (k >>= 1; k; k >>= 1) dense_[k] = op(dense_[k * 2], dense_[k * 2 + 1]);
			return;
		}
		nodes_[id].value = value;
		for (; id; id = nodes_[id].parent) pull(id);
	}
	S prod(int u, int v) const {
		u = index(u), v = index(v);
		if (flat_) {
			auto append = [](const S& part, S value) { return op(std::move(value), part); };
			return dense_apply(u, v, e(), append);
		}
		S left = e(), right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				left = op(left, prefix<true>(u));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		S middle = u <= v ? interval<false>(u, v) : interval<true>(v, u);
		return op(left, op(middle, right));
	}
	template<class X, class Apply> X apply(int u, int v, X value, Apply mapping) const {
		u = index(u), v = index(v);
		if (flat_) return dense_apply(u, v, std::move(value), mapping);
		S right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				for (int id = u; id; id = nodes_[id].parent)
					if (id <= u) value = mapping(back_prefix(id), std::move(value));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		if (u <= v) value = apply_interval<false>(u, v, std::move(value), mapping);
		else value = apply_interval<true>(v, u, std::move(value), mapping);
		return mapping(right, std::move(value));
	}

};

#endif

#line 1 "tree/path_product.hpp"



#line 1 "tree/hld.hpp"



#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <limits>
#include <utility>
#include <vector>

class HLD {
	int shape_ = 0, center_ = 0;
	struct Node { int position, head, parent, offset; };
	std::vector<Node> nodes;
	void check([[maybe_unused]] int v) const { assert(0 <= v && v < n && int(in.size()) == n); }

	void initialize() {
		depth.assign(n, 0);
		size.assign(n, 1);
		heavy.assign(n, -1);
		head.resize(n), in.resize(n), out.resize(n), order.resize(n);
		nodes.resize(n);
		shape_ = center_ = 0;
	}
	void finish(int root) {
		in[root] = 0, out[root] = n, head[root] = root;
		for (int i = n - 2; i >= 0; --i) {
			int v = order[i], p = parent[v];
			depth[v] = depth[p] + 1;
			if (heavy[p] == v) in[v] = in[p] + 1, head[v] = head[p];
			else in[v] = out[p] -= size[v], head[v] = v;
			out[v] = in[v] + size[v];
		}
		for (int v = 0; v < n; ++v) out[v] = in[v] + size[v], order[in[v]] = v;
		parent[root] = -1;
		for (int v = 0; v < n; ++v) nodes[v] = {in[v], in[head[v]], parent[head[v]], depth[v] - in[v]};
	}

  public:
	struct Segment {
		int l, r;
		bool reversed;
	};
	int n;
	std::vector<std::pair<int, int>> edges;
	std::vector<int> parent, depth, size, heavy, head, in, out, order;
	explicit HLD(int n = 0) : n(n) {
		assert(0 <= n && n <= std::numeric_limits<int>::max() / 2);
		edges.reserve(std::max(0, n - 1));
	}
	int add_edge(int u, int v) {
		assert(0 <= u && u < n && 0 <= v && v < n);
		edges.emplace_back(u, v);
		return int(edges.size()) - 1;
	}
	void build(int root = 0) {
		assert(n >= 0 && edges.size() == std::size_t(std::max(0, n - 1)));
		assert(!n || (0 <= root && root < n));
		parent.assign(n, -1);
		initialize();
		if (!n) return;
		std::vector<int> degree(n), neighbors(n);
		for (auto [u, v] : edges) {
			++degree[u], ++degree[v];
			neighbors[u] ^= v, neighbors[v] ^= u;
		}
		center_ = int(std::max_element(degree.begin(), degree.end()) - degree.begin());
		shape_ = degree[center_] <= 2 ? 1 : degree[center_] == n - 1 ? 2 : 0;
		degree[root] += 2;
		int count = 0;
		for (int v = 0; v < n; ++v) if (degree[v] == 1) order[count++] = v;
		for (int i = 0; i < count; ++i) {
			int v = order[i], p = neighbors[v];
			parent[v] = p;
			size[p] += size[v];
			if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
			neighbors[p] ^= v;
			if (--degree[p] == 1) order[count++] = p;
		}
		assert(count == n - 1);
		finish(root);
	}
	void build_from_parents(std::vector<int> parents) {
		assert(parents.size() <= std::size_t(std::numeric_limits<int>::max() / 2));
		n = static_cast<int>(parents.size());
		parent = std::move(parents);
		initialize();
		edges.clear();
		edges.reserve(std::max(0, n - 1));
		if (!n) return;
		std::fill(out.begin(), out.end(), 0);
		int root = -1;
		bool ordered = true;
		for (int v = 0; v < n; ++v) {
			int p = parent[v];
			assert(-1 <= p && p < n && p != v);
			if (p == -1) {
				assert(root == -1);
				root = v;
			} else {
				++out[p];
				edges.emplace_back(p, v);
				ordered &= p < v;
			}
		}
		assert(root != -1);
		int maximum = 0;
		for (int v = 0; v < n; ++v) {
			int degree = out[v] + (v != root);
			if (degree > maximum) maximum = degree, center_ = v;
		}
		shape_ = maximum <= 2 ? 1 : maximum == n - 1 ? 2 : 0;
		if (ordered) {
			for (int v = n - 1; v > 0; --v) {
				int p = parent[v];
				size[p] += size[v];
				if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
				order[n - 1 - v] = v;
			}
		} else {
			int count = 0;
			for (int v = 0; v < n; ++v) if (!out[v] && v != root) order[count++] = v;
			for (int i = 0; i < count; ++i) {
				int v = order[i], p = parent[v];
				size[p] += size[v];
				if (heavy[p] == -1 || size[heavy[p]] < size[v]) heavy[p] = v;
				if (!--out[p] && p != root) order[count++] = p;
			}
			assert(count == n - 1);
		}
		finish(root);
	}
	int lca(int u, int v) const {
		check(u), check(v);
		if (shape_ == 1) return head[u] != head[v] ? order[0] : in[u] < in[v] ? u : v;
		if (shape_ == 2) return u == v ? u : u == order[0] || v == order[0] ? order[0] : center_;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) u = nodes[u].parent;
			else v = nodes[v].parent;
		}
		return nodes[u].position < nodes[v].position ? u : v;
	}
	int dist(int u, int v) const {
		int w = lca(u, v);
		return depth[u] + depth[v] - 2 * depth[w];
	}
	bool is_ancestor(int u, int v) const {
		check(u), check(v);
		return in[u] <= in[v] && in[v] < out[u];
	}
	int ancestor(int v, int k) const {
		check(v);
		if (shape_) {
			if (k < 0 || k > depth[v]) return -1;
			if (!k) return v;
			if (k == depth[v]) return order[0];
			return shape_ == 2 ? parent[v] : order[in[v] - k];
		}
		if (k < 0) return -1;
		int target = nodes[v].position + nodes[v].offset - k;
		if (target < 0) return -1;
		while (nodes[v].head + nodes[v].offset > target) v = nodes[v].parent;
		return order[target - nodes[v].offset];
	}
	int jump(int u, int v, int k) const {
		check(u), check(v);
		if (k < 0) return -1;
		if (!k) return u;
		if (shape_ == 1) {
			if (head[u] == head[v]) {
				int delta = in[v] - in[u];
				if (k > (delta < 0 ? -delta : delta)) return -1;
				return order[in[u] + (delta < 0 ? -k : k)];
			}
			int a = depth[u], b = depth[v];
			if (k > a + b) return -1;
			if (k == a) return order[0];
			return k < a ? order[in[u] - k] : order[in[v] - (a + b - k)];
		}
		if (shape_ == 2) {
			int d = (u != center_) + (v != center_);
			if (u == v || k > d) return -1;
			return k == d ? v : center_;
		}
		int from = u, to = v, target = nodes[u].position + nodes[u].offset - k;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) {
				if (nodes[u].head + nodes[u].offset <= target) return order[target - nodes[u].offset];
				u = nodes[u].parent;
			} else v = nodes[v].parent;
		}
		int w = nodes[u].position < nodes[v].position ? u : v;
		int d = nodes[w].position + nodes[w].offset;
		int length = nodes[from].position + nodes[from].offset + nodes[to].position + nodes[to].offset - 2 * d;
		if (k > length) return -1;
		if (target >= d) return order[target - nodes[u].offset];
		return ancestor(to, length - k);
	}
	std::pair<int, int> subtree(int v, bool edge = false) const {
		check(v);
		return {in[v] + edge, out[v]};
	}
	int edge_vertex(int id) const {
		assert(0 <= id && id < int(edges.size()));
		auto [u, v] = edges[id];
		check(u);
		return parent[u] == v ? u : v;
	}
	template <class F>
	void for_each_path(int u, int v, F &&f, bool edge = false) const {
		check(u), check(v);
		if (shape_ == 1) {
			int x = in[u], y = in[v];
			if (head[u] == head[v]) {
				int l = std::min(x, y) + edge, r = std::max(x, y) + 1;
				if (l < r) f(l, r, x > y);
			} else if (head[u] == order[0]) {
				if (edge < x + 1) f(int(edge), x + 1, true);
				f(in[head[v]], y + 1, false);
			} else {
				f(in[head[u]], x + 1, true);
				if (edge < y + 1) f(int(edge), y + 1, false);
			}
			return;
		}
		std::array<std::array<int, 2>, 32> right;
		int count = 0;
		while (nodes[u].head != nodes[v].head) {
			if (nodes[u].position > nodes[v].position) {
				f(nodes[u].head, nodes[u].position + 1, true);
				u = nodes[u].parent;
			} else {
				assert(count < int(right.size()));
				right[count++] = {nodes[v].head, nodes[v].position + 1};
				v = nodes[v].parent;
			}
		}
		bool reversed = nodes[u].position > nodes[v].position;
		int l = std::min(nodes[u].position, nodes[v].position) + edge, r = std::max(nodes[u].position, nodes[v].position) + 1;
		if (l < r) f(l, r, reversed);
		while (count) {
			auto [a, b] = right[--count];
			f(a, b, false);
		}
	}
	std::vector<Segment> path(int u, int v, bool edge = false) const {
		std::vector<Segment> result;
		for_each_path(u, v, [&](int l, int r, bool reversed) { result.push_back({l, r, reversed}); }, edge);
		return result;
	}
};


#line 5 "tree/path_product.hpp"
#include <bit>
#include <cstdint>
#include <type_traits>

template<class S, auto op, auto e, auto reverse = nullptr> class TreePathProduct {
	static constexpr bool has_reverse = !std::is_same_v<decltype(reverse), std::nullptr_t>;
	struct Empty {};
	struct Back { S prod, prefix; };
	using BackStorage = std::conditional_t<has_reverse, Empty, Back>;
	struct alignas(64) Node {
		int left = 0, right = 0, parent = 0, head = 0, up = 0;
		unsigned char priority = 0;
		S value, prod, prefix;
		[[no_unique_address]] BackStorage back;
		static BackStorage make_back(const S& value) {
			if constexpr (has_reverse) return {};
			else return {value, value};
		}
		explicit Node(const S& v) : value(v), prod(v), prefix(v), back(make_back(v)) {}
	};
	std::vector<Node> nodes_;
	std::vector<int> index_;
	std::vector<S> dense_;
	std::vector<int> dense_heads_, dense_up_;
	unsigned base_ = 1;
	bool flat_ = false;
	int index(int v) const {
		assert(0 <= v && v < size());
		return index_[v];
	}
	void pull(int v) {
		auto& node = nodes_[v];
		node.prefix = node.left ? op(nodes_[node.left].prod, node.value) : node.value;
		node.prod = node.right ? op(node.prefix, nodes_[node.right].prod) : node.prefix;
		if constexpr (!has_reverse) {
			node.back.prefix = node.left ? op(node.value, nodes_[node.left].back.prod) : node.value;
			node.back.prod = node.right ? op(nodes_[node.right].back.prod, node.back.prefix) : node.back.prefix;
		}
	}
	S back_prod(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prod);
		else return nodes_[v].back.prod;
	}
	S back_prefix(int v) const {
		if constexpr (has_reverse) return reverse(nodes_[v].prefix);
		else return nodes_[v].back.prefix;
	}
	template<bool Reverse> S dense_interval(int l, int r) const {
		unsigned a = base_ + unsigned(l - 1), b = base_ + unsigned(r);
		S left = e(), right = e();
		while (a < b) {
			if (a & 1) left = op(left, dense_[a++]);
			if (b & 1) right = op(dense_[--b], right);
			a >>= 1, b >>= 1;
		}
		S result = op(left, right);
		if constexpr (Reverse && has_reverse) return reverse(result);
		else return result;
	}
	template<bool Reverse> S prefix(int v) const {
		S result = e();
		int last = v;
		for (; v; v = nodes_[v].parent) if (v <= last) {
			if constexpr (Reverse) result = op(result, back_prefix(v));
			else result = op(nodes_[v].prefix, result);
		}
		return result;
	}
	template<bool Reverse> S interval(int l, int r) const {
		int a = l, b = r;
		S left = e(), right = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) left = op(op(back_prod(node.right), node.value), left);
					else left = op(left, op(node.value, nodes_[node.right].prod));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) right = op(right, back_prefix(b));
					else right = op(nodes_[b].prefix, right);
				}
				b = nodes_[b].parent;
			}
		}
		if constexpr (Reverse) return op(right, op(nodes_[a].value, left));
		else return op(left, op(nodes_[a].value, right));
	}

	template<bool Reverse, class X, class Apply> X apply_interval(int l, int r, X value, Apply& mapping) const {
		int a = l, b = r;
		S outside = e();
		while (a != b) {
			if (nodes_[a].priority < nodes_[b].priority) {
				if (l <= a) {
					const auto& node = nodes_[a];
					if constexpr (Reverse) outside = op(op(back_prod(node.right), node.value), outside);
					else value = mapping(nodes_[node.right].prod, mapping(node.value, std::move(value)));
				}
				a = nodes_[a].parent;
			} else {
				if (b <= r) {
					if constexpr (Reverse) value = mapping(back_prefix(b), std::move(value));
					else outside = op(nodes_[b].prefix, outside);
				}
				b = nodes_[b].parent;
			}
		}
		return mapping(outside, mapping(nodes_[a].value, std::move(value)));
	}

	template<class X, class Apply> X dense_apply(int u, int v, X value, Apply& mapping) const {
		std::array<std::pair<int, int>, 4> right;
		int count = 0;
		while (dense_heads_[u] != dense_heads_[v]) {
			if (dense_heads_[u] > dense_heads_[v]) {
				value = mapping(dense_interval<true>(dense_heads_[u], u), std::move(value));
				u = dense_up_[u];
			} else {
				assert(count < int(right.size()));
				right[count++] = {dense_heads_[v], v};
				v = dense_up_[v];
			}
		}
		if (u <= v) value = mapping(dense_interval<false>(u, v), std::move(value));
		else value = mapping(dense_interval<true>(v, u), std::move(value));
		while (count) {
			auto [l, r] = right[--count];
			value = mapping(dense_interval<false>(l, r), std::move(value));
		}
		return value;
	}

  public:
	TreePathProduct() = default;
	TreePathProduct(const HLD& tree, const std::vector<S>& values) { build(tree, values); }
	void build(const HLD& tree, const std::vector<S>& values) {
		assert(values.size() == std::size_t(tree.n) && tree.order.size() == values.size());
		index_.resize(values.size());
		nodes_.clear();
		for (int v : tree.order) index_[v] = tree.in[v] + 1;
		flat_ = false;
		dense_.clear();
		if constexpr (has_reverse) {
			std::vector<unsigned char> levels(values.size());
			unsigned maximum = 0;
			for (int v : tree.order) {
				if (tree.parent[v] != -1) levels[v] = levels[tree.parent[v]] + (tree.head[v] == v);
				maximum = std::max(maximum, unsigned(levels[v]));
			}
			flat_ = maximum <= 4;
		}
		if (flat_) {
			base_ = std::bit_ceil(unsigned(size()) ? unsigned(size()) : 1U);
			dense_.assign(std::size_t(base_) * 2, e());
			dense_heads_.assign(values.size() + 1, 0), dense_up_.assign(values.size() + 1, 0);
			for (int v : tree.order) {
				int id = index_[v], h = tree.head[v];
				dense_heads_[id] = index_[h];
				dense_up_[id] = tree.parent[h] == -1 ? 0 : index_[tree.parent[h]];
				dense_[base_ + unsigned(id - 1)] = values[v];
			}
			for (unsigned i = base_ - 1; i; --i) dense_[i] = op(dense_[i * 2], dense_[i * 2 + 1]);
			return;
		}
		dense_heads_.clear(), dense_up_.clear();
		nodes_.reserve(values.size() + 1);
		nodes_.emplace_back(e());
		for (int v : tree.order) nodes_.emplace_back(values[v]);
		std::vector<int> stack, order;
		order.reserve(values.size());
		for (int v : tree.order) {
			if (tree.head[v] != v) continue;
			stack.clear();
			unsigned position = 0;
			int head = index_[v], up = tree.parent[v] == -1 ? 0 : index_[tree.parent[v]];
			for (int u = v; u != -1; u = tree.heavy[u]) {
				int id = index_[u], weight = tree.size[u] - (tree.heavy[u] == -1 ? 0 : tree.size[tree.heavy[u]]);
				auto& node = nodes_[id];
				node.head = head, node.up = up;
				node.priority = std::bit_width(position ^ (position + unsigned(weight)));
				position += unsigned(weight);
				int last = 0;
				while (!stack.empty() && nodes_[stack.back()].priority < node.priority) {
					last = stack.back();
					stack.pop_back();
				}
				node.left = last;
				if (last) nodes_[last].parent = id;
				if (!stack.empty()) nodes_[stack.back()].right = id, node.parent = stack.back();
				stack.push_back(id);
			}
			order.push_back(stack[0]);
		}
		for (std::size_t i = 0; i < order.size(); ++i) {
			const auto& node = nodes_[order[i]];
			if (node.left) order.push_back(node.left);
			if (node.right) order.push_back(node.right);
		}
		for (auto it = order.rbegin(); it != order.rend(); ++it) pull(*it);
	}
	int size() const { return int(index_.size()); }
	bool empty() const { return index_.empty(); }
	const S& get(int v) const {
		int id = index(v);
		return flat_ ? dense_[base_ + unsigned(id - 1)] : nodes_[id].value;
	}
	void set(int v, const S& value) {
		int id = index(v);
		if (flat_) {
			unsigned k = base_ + unsigned(id - 1);
			dense_[k] = value;
			for (k >>= 1; k; k >>= 1) dense_[k] = op(dense_[k * 2], dense_[k * 2 + 1]);
			return;
		}
		nodes_[id].value = value;
		for (; id; id = nodes_[id].parent) pull(id);
	}
	S prod(int u, int v) const {
		u = index(u), v = index(v);
		if (flat_) {
			auto append = [](const S& part, S value) { return op(std::move(value), part); };
			return dense_apply(u, v, e(), append);
		}
		S left = e(), right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				left = op(left, prefix<true>(u));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		S middle = u <= v ? interval<false>(u, v) : interval<true>(v, u);
		return op(left, op(middle, right));
	}
	template<class X, class Apply> X apply(int u, int v, X value, Apply mapping) const {
		u = index(u), v = index(v);
		if (flat_) return dense_apply(u, v, std::move(value), mapping);
		S right = e();
		while (nodes_[u].head != nodes_[v].head) {
			if (nodes_[u].head > nodes_[v].head) {
				for (int id = u; id; id = nodes_[id].parent)
					if (id <= u) value = mapping(back_prefix(id), std::move(value));
				u = nodes_[u].up;
			} else {
				right = op(prefix<false>(v), right);
				v = nodes_[v].up;
			}
		}
		if (u <= v) value = apply_interval<false>(u, v, std::move(value), mapping);
		else value = apply_interval<true>(v, u, std::move(value), mapping);
		return mapping(right, std::move(value));
	}

};