noya2_Library

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

View the Project on GitHub noya2ruler/noya2_Library

:heavy_check_mark: Heavy Light Decomposition
(tree/heavy_light_decomposition.hpp)

Heavy Light Decomposition

HLD のライブラリです。時空間にかなり気を遣って実装していますが、一般的な実装では通らないがこのライブラリを使うと通るような問題はないかもしれません。HLD 自体にあまり関係のないメンバ関数については書いていません（標準入力からの木の入力、CSR の構築など）。

(constructor)

hld_tree (); // (1)
hld_tree (int n, int root = 0); // (2)
hld_tree (const std::vector<int> &par); // (3)
hld_tree (std::vector<int> &&par); // (4)

備考

add_edge

void add_edge(int u, int v);

計算量

備考

leader

int leader(int v) const ;

計算量

la

int la(int v, int d) const ;

計算量

lca

int lca(int u, int v) const ;

計算量

dist

int dist(int u, int v) const ;

計算量

jump

int jump(int from, int to, int d) const ;

計算量

parent

int parent(int v) const ;

計算量

index

int index(int vertex) const ; // (1)
int index_edge(int vertex1, int vertex2) const ; // (2)

計算量

備考

subtree_size

int subtree_size(int v) const ;

計算量

subtree_l

int subtree_l(int v) const ;

計算量

subtree_r

int subtree_r(int v) const ;

計算量

備考

is_in_subtree

bool is_in_subtree(int r, int v) const ;

計算量

dist_table

std::vector<int> dist_table(int s) const ;

計算量

備考

diameter

std::tuple<int, int, int> diameter() const ;

計算量

path

std::vector<int> path(int from, int to) const ;

計算量

備考

path_query

template<bool vertex = true>
void path_query(int u, int v, auto f) const ;

計算量

備考

よくある使い方を示す。

非可換クエリ

モノイド $(S,\oplus,e)$ が非可換のときを扱う。

列に向きがあるため、一般のモノイド $(S,\oplus,e)$ に対して、列を逆向きにするという作用を備えたモノイド $(S\times S, \oplus’, (e,e))$ を扱うことにすると良い。ここで $\oplus’$ は

\[\oplus'((a,a'),(b,b'))=(\oplus (a,b),\oplus (b',a'))\]

で定められる。逆向きにする作用は $(a, b) \mapsto (b, a)$ である。

以下では S,op,e が $(S\times S, \oplus’, (e,e))$ に、 reversed が逆向きにする作用に対応している。

l, r の大小関係で向きの判定を行い、左右からの累積を持っておくと良い。

hld_tree g;
segtree<S,op,e> seg;
S sml = e(), smr = e();
auto f = [&](int l, int r){
    if (l < r){
        smr = op(seg.prod(l,r), smr);
    }
    else {
        sml = op(sml, reversed(seg.prod(r,l)));
    }
};
g.path_query(u, v, f);
S ans = op(sml, smr);

なお、順向きのものだけ取り出すには $(a,b)\mapsto a$ を適用すると良い。

可換クエリ

簡単のためパス上の辺に書かれた整数の総和を扱うことにする。

可換なので、左右や逆向きの区別は必要ない。この場合でも l > r となるものが与えられることに注意せよ。

hld_tree g;
fenwick_tree<int> fen;
int ans = 0;
auto f = [&](int l, int r){
    if (l > r) std::swap(l, r);
    ans += fen.sum(l, r);
};
g.path_query<false>(u, v, f);
// ans is calculated

なお、頂点重みにするには、単に <false> を <true> に書き換えればよい。あるいは <false> を消してしまっても良い（デフォルトで vertex = true である）。

compressed_tree

std::pair<std::vector<int>, std::vector<int>> compressed_tree(std::vector<int> vs) const ;

計算量

備考

省メモリで爆速な HLD のための工夫

木の頂点数を $n$ とする。

本ライブラリの利点

本ライブラリの欠点

本ライブラリが目指したもの

非常に高速に動作する HLD を用いて、木の基本的なクエリを処理することを目指した。

多くの場合は例えば親の取得やパスの復元がボトルネックになることはない。必要な空間に対して線形の時間で取得するものは、高速にするメリットがあまりない（それでも、時空間に気を遣って実装してある）。一方、LCA, LA, jump, dist などは部分問題によく登場する割に、二分探索内でクエリを飛ばすなどして、呼び出しにかかる $\log$ が $2$ つ目の $\log$ になることも多い。したがって、そのようなクエリに対しては高速な実装を提供したい、という動機があった。

また、HLD の構築は、線形時空間にしては重い。depth, sub_size, in_time, out_time, heavy_path_top, adjacent_list, euler_tour, ... （変数名で意図がわかるようにしているつもり） 持っておいて損がないようなものはたくさんあるが、無駄が多い。Library Checker 上の実装やインターネット上に存在する実装を見ても、平均して $8n$ くらいの空間を使用しているようだ（要出典）。空間の重さは時間の重さに直結する。同等の計算回数で達成できるなら、空間を抑えた方がキャッシュが効きやすく、高速に動作する。配列へのランダムアクセスの回数を抑えることも効果的である。

非再帰 euler tour

根付き木 $T$ 上の euler tour を考える。根を $r$ とする。

euler tour で得る配列 down は、euler tour で頂点 v を訪れる時刻が down[v] であることを表す。教科書的な方法は以下の通り。子配列を child[v] とする。

std::vector<std::vector<int>> child;
std::vector<int> down(n);
int t = 0;
auto dfs = [&](auto self, int v) -> void {
    down[v] = t++;
    for (int u : child[v]){
        self(self, u);
    }
};
dfs(dfs,r);

この方法は深さ $O(n)$ の再帰呼び出しが行われる上、二次元配列 child へのアクセスも含み、低速である。（それでも、通常は無視できるほど高速に動作する）

この方法は top-down に時刻を決定するため、再帰呼び出しを回避しやすい bottom-up への書き換えが難しい。

この問題を巧妙に解消しているのが次の実装である。 $T$ の親配列を par 、任意の topological order を ord とする。ここで topological order とは、任意の頂点がその親よりも後に現れるような、すべての頂点の順列である。また、 std::exchange(int &a, int b) は a に b を代入し、元の a を返す関数である。

std::vector<int> down(n);
for (int v : ord | std::views::drop(1) | std::views::reverse){
    down[par[v]] += down[v] + 1;
}
for (int v : ord | std::views::drop(1)){
    down[v] = std::exchange(down[par[v]], down[par[v]] - down[v] - 1);
}

$1$ つ目の for 文の直後では、 down には部分木サイズ $-1$ が格納されている。これを、次のように読み替える： v から euler tour を始めたとき、最後に訪れる頂点の時刻が down[v] である。

$2$ つ目の for 文では、次の解釈に従って値を変更する： par[v] から euler tour を始めるとして、 par[v] の子のうち v には最後に潜るとし、par[v] から v を切り離す。最後の時刻を親から引き継ぎ、親の持つ最後の時刻から自身の部分木サイズを減じるのである。

結局は bottom-up に部分木サイズを求め、top-down に euler tour の時刻を求めているだけである。ただし、euler tour の順に頂点 v を訪れているわけではないことに注意せよ。

さて、よくある根付き木の入力として、親配列が par[i] < i を満たすように与えられることがある。このとき、 ord として $(0,1,\dots ,n-1)$ を取ることができ、 par へのアクセスがシーケンシャルになると同時に ord の保持の必要がなくなる。このときは非常に高速に動作する。

std::vector<int> down(n);
for (int v = n - 1; v >= 1; v--){
    down[par[v]] += down[v] + 1;
}
for (int v = 1; v < n; v++){
    down[v] = std::exchange(down[par[v]], down[par[v]] - down[v] - 1);
}

辺配列から親配列と topological order を取得

辺が親へ向かう有向辺であるなら、葉から切り離していく要領で簡単に計算できる。問題になるのは、よくある木の入力として与えられるような、辺の端点に順序のない場合である。

教科書的な方法は、隣接リストを作り、再帰呼び出しを行なって euler tour する方法である。しかしこれは低速なので、別の方法を用いる。

有向辺のときと同様に次数を管理しながら実装できる。葉をキューに入れて実装するのが普通だが、キューに入る順番がそのまま topological order の逆順になるため、キューを配列で模倣して実装する。

辺配列を std::vector<std::pair<int, int>> edges とする。

std::vector<int> deg(n), par(n), ord(n);
for (auto [u, v] : edges){
    deg[u]++;
    deg[v]++;
    par[u] ^= v;
    par[v] ^= u;
}
int back = 0;
for (int v = 0; v < n; v++){
    if (deg[v] == 1){
        ord[back++] = v;
    }
}
for (int front = 0; front < n-1; front++){
    int v = ord[front];
    par[par[v]] ^= v;
    if (--deg[par[v]] == 1){
        ord[back++] = par[v];
    }
}
int root = ord[n - 1];
par[root] = -1;
std::ranges::reverse(ord);

この実装は $n=1$ のときも正しく動作する。

重要な点は par を、隣接する頂点の番号の総 XOR にしていることである。葉を切り離していく際に deg に加えて par も更新することで、 v が葉になって以降は par[v] が親の番号となる。

この実装は木の根が後から定まるが、根を前もって決定したい場合はキューに最後まで入れないようにすれば良い。

また、 edges の保持は、現実的には全く必要ないこともわかる。辺 u,v の追加のタイミングで deg[u],deg[v],par[u],par[v] を書き換えれば良いし、その後必要になることはない。（辺配列は親配列から簡単に得ることができる。）

ここからさらに euler tour までしたい場合は、 for (int front ...) の内部で部分木サイズまで計算することができ、追加の for 文は $1$ 回で済む。

省メモリのための戦略

本ライブラリは次の $4$ つの長さ $n$ の std::vector<int> を保持する。

HLD の主な目的を、パスクエリを処理すること、LCA を求めることとする。このときに必要なのは down, nxt のみである。

int leader(int v){
    return nxt[v] < 0 ? v : nxt[v];
}
int lca(int u, int v){
    while (leader(u) != leader(v)){
        if (down[u] > down[v]) std::swap(u, v);
        v = ~nxt[leader(v)];
    }
    return down[u] < down[v] ? u : v;
}
void path_query(int u, int v, auto f){
    while (leader(u) != leader(v)){
        if (down[u] < down[v]){
            f(down[leader(v)], down[v] + 1);
            v = ~nxt[leader(v)];
        }
        else {
            f(down[u] + 1, down[leader(u)]);
            u = ~nxt[leader(u)];
        }
    }
    if (down[u] < down[v]){
        f(down[u], down[v] + 1);
    }
    else {
        f(down[u] + 1, down[v]);
    }
}

down, nxt のみで処理できる主なクエリは次のとおり。

これだけでも十分に見えるが、HLD の前計算（構築）はより多くの空間を必要とする。

親配列から HLD を行って down, nxt を得るには、追加で長さ $n$ の std::vector<int> が必要に見える。（相当な時間を費やしたが、追加の空間を $o(n)$ にする方法は思いつかなかった。）素直な方法では、追加で使用した配列は sub に相当するものを格納することになる。よって、はじめから保持しておき、より広いクエリに対処することにした。

down, nxt, sub で新たに処理できる主なクエリは次のとおり。

辺配列からの構築を行う際には、さらに追加で長さ $n$ の std::vector<int> が必要に見える。追加で LA, jump クエリを処理することを考えると、 down の逆順列を持つのが都合が良い。そこで tour も保持することにした。

down, nxt, sub, tour で新たに処理できる主なクエリは次のとおり。

教科書的な HLD の実装では親配列をそのまま保持することが多い（要出典）。親配列から HLD を行う場合はなおさらである。本ライブラリは親配列を保持しないが、次のように実装できる。

int parent(int v){
    return nxt[v] < 0 ? ~nxt[v] : tour[down[v] - 1];
}

heavy path 上では down の値が連続していることに注目すると、正当性が確認できる。

爆速な LCA

本ライブラリの実装は次のとおり。

int lca(int u, int v) const {
    int du = down[u], dv = down[v];
    if (du > dv){
        std::swap(du, dv);
        std::swap(u, v);
    }
    if (dv < du + sub[u]){
        return u;
    }
    while (du < dv){
        v = ~nxt[leader(v)];
        dv = down[v];
    }
    return v;
}

この実装の正当性の確認の前に、LCA クエリの基本的な処理方法について確認する。

nxt2 を nxt のよくある定義に置き換えたバージョンとして次のように定義する。

down, nxt2 に加えて親配列 par と深さ配列 dep を保持したときの教科書的な LCA の実装は次のとおり：

int lca(int u, int v){
    while (nxt2[u] != nxt2[v]){
        if (down[u] > down[v]) std::swap(u, v);
        v = par[nxt2[v]];
    }
    return dep[u] < dep[v] ? u : v;
}

この実装は、アクセスしている配列の種類が多すぎる。

次の点に注目せよ。

nxt2 の情報を完全に保持しながら par へのアクセスに応えたものが本ライブラリの nxt である。 nxt2[v] は leader(v) として得ることができる。

ここまでの工夫で得られる実装が 省メモリのための戦略 に載せたものである。

さらに次の点に注目せよ。

次の疑似コードによる最も基礎的な実装の正当性を確認せよ。

int lca(int u, int v){
    if u が v の祖先 : return u
    while v が u の祖先でない {
        v = par[v]
    }
    return v
}

これを模倣したものが本ライブラリの実装である。

実は、教科書的な実装において swap が発生するのは高々 $2$ 回である。down[u] < down[v] とすれば、まず v が lca(u,v) まで登り、次に u が v に代わって登っていただけだ。これを用いれば if, swap を取り除くことができ、while の内部はずっと単純になる。また、 u が登らなくて良くなったため、通常の実装と単純に比較して、約半分程度の計算で済むようになる。

その代わりに sub を必要とすることに注意せよ。はじめに祖先の判定を行わなかった場合、v は lca(u,v) を通り過ぎてしまう。

LCA クエリを処理するだけなら、 nxt を次の nxt3 に置き換えれば良い。

なお、この工夫は heavy path の端点を失うため、パスクエリとは共存しない。

爆速な HLD

euler tour の時刻 down を非再帰で取得する方法を紹介した。本ライブラリでは、実装の都合上、少し異なる方法で行っている：

のではなく、

としている。

先に紹介したものは、子に潜る順番に制約がなかった。しかし、HLD では最もサイズの大きい部分木を持つ子にはじめに潜る必要がある1。よって、部分木サイズの最大値を保持する必要がある。

親配列 nxt から down, nxt, sub, tour を行うために for 文を $4$ 回まわしている。

  1. 部分木サイズ sub と、子の部分木サイズの最大値 down の計算
  2. 自身が親の heavy child であるかの判定
  3. sub を用いた down と nxt の構築
  4. tour を down の逆順列であるとして計算

辺配列からは for 文を $7$ 回まわしている。

  1. 次数 down と隣接頂点の総 XOR nxt の計算
  2. 次数 $1$ の頂点をキュー tour に入れる
  3. 次数 $1$ の頂点を切り離しながら親配列 nxt と topological order tour を取得
  4. 以降は親配列からの構築と同じ

ここにはもう少し工夫の余地があると思っているが、特に思いつかなかった。

Verified with

Code

#pragma once

#include <vector>
#include <algorithm>
#include <stack>
#include <ranges>
#include <cassert>
#include <utility>

// #include "data_structure/csr.hpp"

namespace noya2 {

struct hld_tree {
    int n, root;
    std::vector<int> down, nxt, sub, tour;
	// noya2::internal::csr<int> childs;

    // default constructor (nop)
    hld_tree () {}

    // tree with _n node
    // after construct, call input_edges / input_parents / add_edge _n - 1 times
    hld_tree (int _n, int _root = 0) : n(_n), root(_root), down(n), nxt(n), sub(n, 1), tour(n) {
        if (n == 1){
            nxt[0] = -1;
            down[0] = -1;
            build_from_parents();
        }
    }

    // par[i] < i, par[0] == -1
    hld_tree (const std::vector<int> &par) : n(par.size()), root(0), down(n, -1), nxt(par), sub(n, 1), tour(n){
        build_from_parents();
    }

    // par[i] < i, par[0] == -1
    hld_tree (std::vector<int> &&par) : n(par.size()), root(0), down(n, -1), sub(n, 1), tour(n) {
        nxt.swap(par);
        build_from_parents();
    }

    // distinct unweighted undirected n - 1 edges of tree 
    hld_tree (const std::vector<std::pair<int, int>> &es, int _root = 0) : n(es.size() + 1), root(_root), down(n), nxt(n), sub(n, 1), tour(n) {
        for (auto &[u, v] : es){
            down[u]++;
            down[v]++;
            nxt[u] ^= v;
            nxt[v] ^= u;
        }
        build_from_edges();
    }

    // input parents from cin
    template<int indexed = 1>
    void input_parents(){
        // using std::cin;
        nxt[0] = -1;
        for (int u = 1; u < n; u++){
            cin >> nxt[u];
            nxt[u] -= indexed;
        }
        build_from_parents();
    }

    // input n - 1 edges from cin
    template<int indexed = 1>
    void input_edges(){
        // using std::cin;
        for (int i = 1; i < n; i++){
            int u, v; cin >> u >> v;
            u -= indexed;
            v -= indexed;
            down[u]++;
            down[v]++;
            nxt[u] ^= v;
            nxt[v] ^= u;
        }
        build_from_edges();
    }

    void add_edge(int u, int v){
        down[u]++;
        down[v]++;
        nxt[u] ^= v;
        nxt[v] ^= u;
        // use tour[0] as counter
        if (++tour[0] == n - 1){
            build_from_edges();
        }
    }

    size_t size() const {
        return n;
    }

    // top vertex of heavy path which contains v
    int leader(int v) const {
        return nxt[v] < 0 ? v : nxt[v];
    }

    // level ancestor
    // ret is ancestor of v, dist(ret, v) == d
    // if d > depth(v), return -1
    int la(int v, int d) const {
        while (v != -1){
            int u = leader(v);
            if (down[v] - d >= down[u]){
                v = tour[down[v] - d];
                break;
            }
            d -= down[v] - down[u] + 1;
            v = (u == root ? -1 : ~nxt[u]);
        }
        return v;
    }

    // lowest common ancestor of u and v
    int lca(int u, int v) const {
        int du = down[u], dv = down[v];
        if (du > dv){
            std::swap(du, dv);
            std::swap(u, v);
        }
        if (dv < du + sub[u]){
            return u;
        }
        while (du < dv){
            v = ~nxt[leader(v)];
            dv = down[v];
        }
        return v;
    }

    // distance from u to v
    int dist(int u, int v) const {
        int _dist = 0;
        while (leader(u) != leader(v)){
            if (down[u] > down[v]) std::swap(u, v);
            _dist += down[v] - down[leader(v)] + 1;
            v = ~nxt[leader(v)];
        }
        _dist += std::abs(down[u] - down[v]);
        return _dist;
    }

    // d times move from to its neighbor (direction of to)
    // if d > dist(from, to), return -1
    int jump(int from, int to, int d) const {
        int _from = from, _to = to;
        int dist_from_lca = 0, dist_to_lca = 0;
        while (leader(_from) != leader(_to)){
            if (down[_from] > down[_to]){
                dist_from_lca += down[_from] - down[leader(_from)] + 1;
                _from = ~nxt[leader(_from)];
            }
            else {
                dist_to_lca += down[_to] - down[leader(_to)] + 1;
                _to = ~nxt[leader(_to)];
            }
        }
        if (down[_from] > down[_to]){
            dist_from_lca += down[_from] - down[_to];
        }
        else {
            dist_to_lca += down[_to] - down[_from];
        }
        if (d <= dist_from_lca){
            return la(from, d);
        }
        d -= dist_from_lca;
        if (d <= dist_to_lca){
            return la(to, dist_to_lca - d);
        }
        return -1;
    }

    // parent of v (if v is root, return -1)
    int parent(int v) const {
        if (v == root) return -1;
        return (nxt[v] < 0 ? ~nxt[v] : tour[down[v] - 1]);
    }

    // visiting time in euler tour
    // usage : seg.set(index(v), X[v])
    int index(int vertex) const {
        return down[vertex];
    }
    // usage : seg.set(index_edge(e.u, e.v), e.val)
    int index(int vertex1, int vertex2) const {
        return std::max(down[vertex1], down[vertex2]);
    }

    // subtree size of v
    int subtree_size(int v) const {
        return sub[v];
    }

    // prod in subtree v : seg.prod(subtree_l(v), subtree_r(v))
    int subtree_l(int v) const {
        return down[v];
    }
    int subtree_r(int v) const {
        return down[v] + sub[v];
    }

    // v is in subtree r
    bool is_in_subtree(int r, int v) const {
        return subtree_l(r) <= subtree_l(v) && subtree_r(v) <= subtree_r(r);
    }
    
    // distance table from s
    std::vector<int> dist_table(int s) const {
        std::vector<int> table(n, -1);
        table[s] = 0;
        while (s != root){
            table[parent(s)] = table[s] + 1;
            s = parent(s);
        }
        for (int v : tour){
            if (table[v] == -1){
                table[v] = table[parent(v)] + 1;
            }
        }
        return table;
    }

    // dist, v1, v2
    std::tuple<int, int, int> diameter() const {
        std::vector<int> dep = dist_table(root);
        int v1 = std::ranges::max_element(dep) - dep.begin();
        std::vector<int> fromv1 = dist_table(v1);
        int v2 = std::ranges::max_element(fromv1) - fromv1.begin();
        return {fromv1[v2], v1, v2};
    }

    // vertex array {from, ..., to}
    std::vector<int> path(int from, int to) const {
        int d = dist(from, to);
        std::vector<int> _path(d + 1);
        int front = 0, back = d;
        while (from != to){
            if (down[from] > down[to]){
                _path[front++] = from;
                from = parent(from);
            }
            else {
                _path[back--] = to;
                to = parent(to);
            }
        }
        _path[front] = from;
        return _path;
    }

    // path decomposition and query (vertex weighted)
    // if l < r, decsending order tour[l, r)
    // if l > r, acsending order tour(l, r]
    template<bool vertex = true>
    void path_query(int u, int v, auto f) const {
        while (leader(u) != leader(v)){
            if (down[u] < down[v]){
                f(down[leader(v)], down[v] + 1);
                v = ~nxt[leader(v)];
            }
            else {
                f(down[u] + 1, down[leader(u)]);
                u = ~nxt[leader(u)];
            }
        }
        if constexpr (vertex){
            if (down[u] < down[v]){
                f(down[u], down[v] + 1);
            }
            else {
                f(down[u] + 1, down[v]);
            }
        }
        else {
            if (down[u] != down[v]){
                f(down[u] + 1, down[v] + 1);
            }
        }
    }

    // {parent, mapping} : cptree i is correspond to tree mapping[i]. parent[i] is parent of i in cptree.
    // parent[i] < i, parent[0] == -1
	std::pair<std::vector<int>, std::vector<int>> compressed_tree(std::vector<int> vs) const {
        if (vs.empty()){
            return {{},{}};
        }
        auto comp = [&](int l, int r){
            return down[l] < down[r];
        };
		std::ranges::sort(vs, comp);
		int sz = vs.size(); vs.reserve(2*sz);
        for (int i = 0; i < sz-1; i++){
            vs.emplace_back(lca(vs[i], vs[i+1]));
        }
        std::sort(vs.begin() + sz, vs.end(), comp);
        std::ranges::inplace_merge(vs, vs.begin() + sz, comp);
        auto del = std::ranges::unique(vs);
        vs.erase(del.begin(), del.end());
        sz = vs.size();
        std::stack<int> st;
        std::vector<int> par(sz);
        par[0] = -1;
        st.push(0);
        for (int i = 1; i < sz; i++){
            while (!is_in_subtree(vs[st.top()], vs[i])) st.pop();
            par[i] = st.top();
            st.push(i);
        }
        return {par, vs};
	}

/*  CSR

	// build csr for using operator()
	void build_csr(){
		childs = noya2::internal::csr<int>(n, n - 1);
		for (int v = 0; v < n; v++){
			if (v == root) continue;
			childs.add(parent(v), v);
		}
		childs.build();
	}
	const auto operator()(int v) const {
		return childs[v];
	}
	auto operator()(int v){
		return childs[v];
	}
*/

    // hld_tree g;
    // euler tour order : `for (int v : g)`
    // with range_adaptor : `for (int v : g | std::views::reverse)`
    // bottom-up DP : `for (int v : g | std::views::drop(1) | std::views::reverse){ update dp[g.parent(v)] by dp[v] }`
    auto begin() const {
        return tour.begin();
    }
    auto end() const {
        return tour.end();
    }

  private:
    // nxt[v] : parent of v, nxt[0] == -1
    void build_from_parents(){
        for (int u = n - 1; u >= 1; u--){
            int v = nxt[u];
            sub[v] += sub[u];
            down[v] = std::max(down[v], sub[u]);
        }
        for (int u = n - 1; u >= 1; u--){
            int v = nxt[u];
            if (down[v] == sub[u]){
                sub[u] = ~sub[u];
                down[v] = ~down[v];
            }
        }

        sub[0] = ~down[0] + 1;
        down[0] = 0;
        for (int u = 1; u < n; u++){
            int v = nxt[u];
            int nsub = ~down[u] + 1;
            if (sub[u] < 0){
                down[u] = down[v] + 1;
                nxt[u] = (nxt[v] < 0 ? v : nxt[v]);
            }
            else {
                down[u] = down[v] + sub[v];
                sub[v] += sub[u];
                nxt[u] = ~v;
            }
            sub[u] = nsub;
        }

        for (int u = 0; u < n; u++){
            tour[down[u]] = u;
        }
    }

    // down[v] : degree of v
    // nxt[v] : xor prod of neighbor of v
    void build_from_edges(){
        // use tour as queue
        int back = 0;
        for (int u = 0; u < n; u++){
            if (u != root && down[u] == 1){
                tour[back++] = u;
            }
        }
        for (int front = 0; front < n - 1; front++){
            int u = tour[front];
            down[u] = -1;
            int v = nxt[u]; // parent of v
            nxt[v] ^= u;
            if (--down[v] == 1 && v != root){
                tour[back++] = v;
            }
        }
        // check : now, tour is reverse of topological order

        tour.pop_back();

        // check : now, down[*] <= 1
        for (int u : tour){
            int v = nxt[u];
            // subtree size (initialized (1,1,...,1))
            sub[v] += sub[u];
            // heaviest subtree of its child
            down[v] = std::max(down[v], sub[u]);
        }
        for (int u : tour){
            int v = nxt[u];
            // whether u is not the top of heavy path
            if (down[v] == sub[u]){
                sub[u] = ~sub[u];
                down[v] = ~down[v];
            }
        }

        // after appearing v as u (or v == root), 
        // down[v] is the visiting time of euler tour
        // nxt[v] is the lowest vertex of heavy path which contains v
        //   (if v itself, nxt[v] is ~(parent of v))
        // sub[v] + down[v] is the light child's starting time of euler tour
        // note : heavy child's visiting time of euler tour is (the time of its parent) + 1
        sub[root] = ~down[root] + 1;
        down[root] = 0;
        nxt[root] = -1;
        for (int u : tour | std::views::reverse){
            int v = nxt[u];
            int nsub = ~down[u] + 1;
            // heavy child
            if (sub[u] < 0){
                down[u] = down[v] + 1;
                nxt[u] = (nxt[v] < 0 ? v : nxt[v]);
            }
            // light child
            else {
                down[u] = down[v] + sub[v];
                sub[v] += sub[u];
                nxt[u] = ~v;
            }
            sub[u] = nsub;
        }

        // tour is inverse permutation of down
        tour.push_back(0);
        for (int u = 0; u < n; u++){
            tour[down[u]] = u;
        }
    }
};

} // namespace noya2
#line 2 "tree/heavy_light_decomposition.hpp"

#include <vector>
#include <algorithm>
#include <stack>
#include <ranges>
#include <cassert>
#include <utility>

// #include "data_structure/csr.hpp"

namespace noya2 {

struct hld_tree {
    int n, root;
    std::vector<int> down, nxt, sub, tour;
	// noya2::internal::csr<int> childs;

    // default constructor (nop)
    hld_tree () {}

    // tree with _n node
    // after construct, call input_edges / input_parents / add_edge _n - 1 times
    hld_tree (int _n, int _root = 0) : n(_n), root(_root), down(n), nxt(n), sub(n, 1), tour(n) {
        if (n == 1){
            nxt[0] = -1;
            down[0] = -1;
            build_from_parents();
        }
    }

    // par[i] < i, par[0] == -1
    hld_tree (const std::vector<int> &par) : n(par.size()), root(0), down(n, -1), nxt(par), sub(n, 1), tour(n){
        build_from_parents();
    }

    // par[i] < i, par[0] == -1
    hld_tree (std::vector<int> &&par) : n(par.size()), root(0), down(n, -1), sub(n, 1), tour(n) {
        nxt.swap(par);
        build_from_parents();
    }

    // distinct unweighted undirected n - 1 edges of tree 
    hld_tree (const std::vector<std::pair<int, int>> &es, int _root = 0) : n(es.size() + 1), root(_root), down(n), nxt(n), sub(n, 1), tour(n) {
        for (auto &[u, v] : es){
            down[u]++;
            down[v]++;
            nxt[u] ^= v;
            nxt[v] ^= u;
        }
        build_from_edges();
    }

    // input parents from cin
    template<int indexed = 1>
    void input_parents(){
        // using std::cin;
        nxt[0] = -1;
        for (int u = 1; u < n; u++){
            cin >> nxt[u];
            nxt[u] -= indexed;
        }
        build_from_parents();
    }

    // input n - 1 edges from cin
    template<int indexed = 1>
    void input_edges(){
        // using std::cin;
        for (int i = 1; i < n; i++){
            int u, v; cin >> u >> v;
            u -= indexed;
            v -= indexed;
            down[u]++;
            down[v]++;
            nxt[u] ^= v;
            nxt[v] ^= u;
        }
        build_from_edges();
    }

    void add_edge(int u, int v){
        down[u]++;
        down[v]++;
        nxt[u] ^= v;
        nxt[v] ^= u;
        // use tour[0] as counter
        if (++tour[0] == n - 1){
            build_from_edges();
        }
    }

    size_t size() const {
        return n;
    }

    // top vertex of heavy path which contains v
    int leader(int v) const {
        return nxt[v] < 0 ? v : nxt[v];
    }

    // level ancestor
    // ret is ancestor of v, dist(ret, v) == d
    // if d > depth(v), return -1
    int la(int v, int d) const {
        while (v != -1){
            int u = leader(v);
            if (down[v] - d >= down[u]){
                v = tour[down[v] - d];
                break;
            }
            d -= down[v] - down[u] + 1;
            v = (u == root ? -1 : ~nxt[u]);
        }
        return v;
    }

    // lowest common ancestor of u and v
    int lca(int u, int v) const {
        int du = down[u], dv = down[v];
        if (du > dv){
            std::swap(du, dv);
            std::swap(u, v);
        }
        if (dv < du + sub[u]){
            return u;
        }
        while (du < dv){
            v = ~nxt[leader(v)];
            dv = down[v];
        }
        return v;
    }

    // distance from u to v
    int dist(int u, int v) const {
        int _dist = 0;
        while (leader(u) != leader(v)){
            if (down[u] > down[v]) std::swap(u, v);
            _dist += down[v] - down[leader(v)] + 1;
            v = ~nxt[leader(v)];
        }
        _dist += std::abs(down[u] - down[v]);
        return _dist;
    }

    // d times move from to its neighbor (direction of to)
    // if d > dist(from, to), return -1
    int jump(int from, int to, int d) const {
        int _from = from, _to = to;
        int dist_from_lca = 0, dist_to_lca = 0;
        while (leader(_from) != leader(_to)){
            if (down[_from] > down[_to]){
                dist_from_lca += down[_from] - down[leader(_from)] + 1;
                _from = ~nxt[leader(_from)];
            }
            else {
                dist_to_lca += down[_to] - down[leader(_to)] + 1;
                _to = ~nxt[leader(_to)];
            }
        }
        if (down[_from] > down[_to]){
            dist_from_lca += down[_from] - down[_to];
        }
        else {
            dist_to_lca += down[_to] - down[_from];
        }
        if (d <= dist_from_lca){
            return la(from, d);
        }
        d -= dist_from_lca;
        if (d <= dist_to_lca){
            return la(to, dist_to_lca - d);
        }
        return -1;
    }

    // parent of v (if v is root, return -1)
    int parent(int v) const {
        if (v == root) return -1;
        return (nxt[v] < 0 ? ~nxt[v] : tour[down[v] - 1]);
    }

    // visiting time in euler tour
    // usage : seg.set(index(v), X[v])
    int index(int vertex) const {
        return down[vertex];
    }
    // usage : seg.set(index_edge(e.u, e.v), e.val)
    int index(int vertex1, int vertex2) const {
        return std::max(down[vertex1], down[vertex2]);
    }

    // subtree size of v
    int subtree_size(int v) const {
        return sub[v];
    }

    // prod in subtree v : seg.prod(subtree_l(v), subtree_r(v))
    int subtree_l(int v) const {
        return down[v];
    }
    int subtree_r(int v) const {
        return down[v] + sub[v];
    }

    // v is in subtree r
    bool is_in_subtree(int r, int v) const {
        return subtree_l(r) <= subtree_l(v) && subtree_r(v) <= subtree_r(r);
    }
    
    // distance table from s
    std::vector<int> dist_table(int s) const {
        std::vector<int> table(n, -1);
        table[s] = 0;
        while (s != root){
            table[parent(s)] = table[s] + 1;
            s = parent(s);
        }
        for (int v : tour){
            if (table[v] == -1){
                table[v] = table[parent(v)] + 1;
            }
        }
        return table;
    }

    // dist, v1, v2
    std::tuple<int, int, int> diameter() const {
        std::vector<int> dep = dist_table(root);
        int v1 = std::ranges::max_element(dep) - dep.begin();
        std::vector<int> fromv1 = dist_table(v1);
        int v2 = std::ranges::max_element(fromv1) - fromv1.begin();
        return {fromv1[v2], v1, v2};
    }

    // vertex array {from, ..., to}
    std::vector<int> path(int from, int to) const {
        int d = dist(from, to);
        std::vector<int> _path(d + 1);
        int front = 0, back = d;
        while (from != to){
            if (down[from] > down[to]){
                _path[front++] = from;
                from = parent(from);
            }
            else {
                _path[back--] = to;
                to = parent(to);
            }
        }
        _path[front] = from;
        return _path;
    }

    // path decomposition and query (vertex weighted)
    // if l < r, decsending order tour[l, r)
    // if l > r, acsending order tour(l, r]
    template<bool vertex = true>
    void path_query(int u, int v, auto f) const {
        while (leader(u) != leader(v)){
            if (down[u] < down[v]){
                f(down[leader(v)], down[v] + 1);
                v = ~nxt[leader(v)];
            }
            else {
                f(down[u] + 1, down[leader(u)]);
                u = ~nxt[leader(u)];
            }
        }
        if constexpr (vertex){
            if (down[u] < down[v]){
                f(down[u], down[v] + 1);
            }
            else {
                f(down[u] + 1, down[v]);
            }
        }
        else {
            if (down[u] != down[v]){
                f(down[u] + 1, down[v] + 1);
            }
        }
    }

    // {parent, mapping} : cptree i is correspond to tree mapping[i]. parent[i] is parent of i in cptree.
    // parent[i] < i, parent[0] == -1
	std::pair<std::vector<int>, std::vector<int>> compressed_tree(std::vector<int> vs) const {
        if (vs.empty()){
            return {{},{}};
        }
        auto comp = [&](int l, int r){
            return down[l] < down[r];
        };
		std::ranges::sort(vs, comp);
		int sz = vs.size(); vs.reserve(2*sz);
        for (int i = 0; i < sz-1; i++){
            vs.emplace_back(lca(vs[i], vs[i+1]));
        }
        std::sort(vs.begin() + sz, vs.end(), comp);
        std::ranges::inplace_merge(vs, vs.begin() + sz, comp);
        auto del = std::ranges::unique(vs);
        vs.erase(del.begin(), del.end());
        sz = vs.size();
        std::stack<int> st;
        std::vector<int> par(sz);
        par[0] = -1;
        st.push(0);
        for (int i = 1; i < sz; i++){
            while (!is_in_subtree(vs[st.top()], vs[i])) st.pop();
            par[i] = st.top();
            st.push(i);
        }
        return {par, vs};
	}

/*  CSR

	// build csr for using operator()
	void build_csr(){
		childs = noya2::internal::csr<int>(n, n - 1);
		for (int v = 0; v < n; v++){
			if (v == root) continue;
			childs.add(parent(v), v);
		}
		childs.build();
	}
	const auto operator()(int v) const {
		return childs[v];
	}
	auto operator()(int v){
		return childs[v];
	}
*/

    // hld_tree g;
    // euler tour order : `for (int v : g)`
    // with range_adaptor : `for (int v : g | std::views::reverse)`
    // bottom-up DP : `for (int v : g | std::views::drop(1) | std::views::reverse){ update dp[g.parent(v)] by dp[v] }`
    auto begin() const {
        return tour.begin();
    }
    auto end() const {
        return tour.end();
    }

  private:
    // nxt[v] : parent of v, nxt[0] == -1
    void build_from_parents(){
        for (int u = n - 1; u >= 1; u--){
            int v = nxt[u];
            sub[v] += sub[u];
            down[v] = std::max(down[v], sub[u]);
        }
        for (int u = n - 1; u >= 1; u--){
            int v = nxt[u];
            if (down[v] == sub[u]){
                sub[u] = ~sub[u];
                down[v] = ~down[v];
            }
        }

        sub[0] = ~down[0] + 1;
        down[0] = 0;
        for (int u = 1; u < n; u++){
            int v = nxt[u];
            int nsub = ~down[u] + 1;
            if (sub[u] < 0){
                down[u] = down[v] + 1;
                nxt[u] = (nxt[v] < 0 ? v : nxt[v]);
            }
            else {
                down[u] = down[v] + sub[v];
                sub[v] += sub[u];
                nxt[u] = ~v;
            }
            sub[u] = nsub;
        }

        for (int u = 0; u < n; u++){
            tour[down[u]] = u;
        }
    }

    // down[v] : degree of v
    // nxt[v] : xor prod of neighbor of v
    void build_from_edges(){
        // use tour as queue
        int back = 0;
        for (int u = 0; u < n; u++){
            if (u != root && down[u] == 1){
                tour[back++] = u;
            }
        }
        for (int front = 0; front < n - 1; front++){
            int u = tour[front];
            down[u] = -1;
            int v = nxt[u]; // parent of v
            nxt[v] ^= u;
            if (--down[v] == 1 && v != root){
                tour[back++] = v;
            }
        }
        // check : now, tour is reverse of topological order

        tour.pop_back();

        // check : now, down[*] <= 1
        for (int u : tour){
            int v = nxt[u];
            // subtree size (initialized (1,1,...,1))
            sub[v] += sub[u];
            // heaviest subtree of its child
            down[v] = std::max(down[v], sub[u]);
        }
        for (int u : tour){
            int v = nxt[u];
            // whether u is not the top of heavy path
            if (down[v] == sub[u]){
                sub[u] = ~sub[u];
                down[v] = ~down[v];
            }
        }

        // after appearing v as u (or v == root), 
        // down[v] is the visiting time of euler tour
        // nxt[v] is the lowest vertex of heavy path which contains v
        //   (if v itself, nxt[v] is ~(parent of v))
        // sub[v] + down[v] is the light child's starting time of euler tour
        // note : heavy child's visiting time of euler tour is (the time of its parent) + 1
        sub[root] = ~down[root] + 1;
        down[root] = 0;
        nxt[root] = -1;
        for (int u : tour | std::views::reverse){
            int v = nxt[u];
            int nsub = ~down[u] + 1;
            // heavy child
            if (sub[u] < 0){
                down[u] = down[v] + 1;
                nxt[u] = (nxt[v] < 0 ? v : nxt[v]);
            }
            // light child
            else {
                down[u] = down[v] + sub[v];
                sub[v] += sub[u];
                nxt[u] = ~v;
            }
            sub[u] = nsub;
        }

        // tour is inverse permutation of down
        tour.push_back(0);
        for (int u = 0; u < n; u++){
            tour[down[u]] = u;
        }
    }
};

} // namespace noya2
Back to top page