tpl_mo_on_trees.H

Go to the documentation of this file.

 
/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
 
# ifndef TPL_MO_ON_TREES_H
# define TPL_MO_ON_TREES_H
 
# include <ah-graph-concepts.H>
 
# include <algorithm>
# include <cmath>
# include <concepts>
# include <vector>
# include <tpl_array.H>
# include <tpl_dynMapOhash.H>
# include <tpl_dynList.H>
# include <tpl_dynListStack.H>
# include <tpl_mo_algorithm.H>
# include <tpl_tree_node.H>
# include <ah-errors.H>
 
namespace Aleph
{
  template <class GT>
  concept MoTreeGraph = requires(const GT & g,
                                 typename GT::Node * p,
                                 typename GT::Arc * a)
    {
      typename GT::Node;
      typename GT::Arc;
      typename GT::Node::Node_Type;
      typename GT::Node_Arc_Iterator;
      { g.vsize() } -> std::convertible_to<size_t>;
      { g.get_connected_node(a, p) } -> std::same_as<typename GT::Node *>;
      { typename GT::Node_Arc_Iterator(p) };
      { g.for_each_node([](typename GT::Node *) {}) };
    };
 
  // ================================================================
  // Gen_Mo_On_Trees
  // ================================================================
 
  template <AlephGraph GT, class Policy>
    requires MoTreeGraph<GT> and MoPolicy<Policy, typename GT::Node::Node_Type>
  class Gen_Mo_On_Trees
  {
  public:
    using Node = typename GT::Node;
    using Arc = typename GT::Arc;
    using T = typename Node::Node_Type;
    using answer_type = typename Policy::answer_type;
 
  private:
    static constexpr size_t NONE = ~static_cast<size_t>(0);
 
    const GT & g_;
    Node * root_;
    mutable Policy pol_;
    size_t n_ = 0;
 
    // Node <-> ID mapping (external, zero impact on graph)
    Array<Node*> id_to_node_;
    MapOLhash<Node*, size_t> node_to_id_;
    Array<T> node_values_;
 
    // Single-occurrence Euler tour (for subtree queries, size n)
    Array<size_t> tin_;
    Array<size_t> tout_;
    Array<T> flat_sub_;
 
    // Double-occurrence Euler tour (for path queries, size 2n)
    Array<size_t> first_;
    Array<size_t> last_;
    Array<size_t> flat_node_;
 
    // Binary lifting for LCA
    size_t log_n_ = 0;
    Array<size_t> depth_;
    Array<size_t> parent_;
    Array<size_t> up_;  // flattened [log_n_][n_]: up_(k * n_ + v)
 
    // ── Build helpers ──────────────────────────────────────────────
 
    void build()
    {
      n_ = g_.vsize();
      if (n_ == 0)
        return;
 
      // 1. Assign contiguous IDs to nodes (pre-allocate OLhash for n_ entries)
      id_to_node_ = Array<Node*>::create(n_);
      node_values_ = Array<T>::create(n_);
      {
        MapOLhash<Node*, size_t> tmp(n_);
        node_to_id_.swap(tmp);
      }
 
      {
        size_t next_id = 0;
        g_.for_each_node([&](Node * p)
          {
            node_to_id_.insert(p, next_id);
            id_to_node_(next_id) = p;
            node_values_(next_id) = p->get_info();
            ++next_id;
          });
      }
 
      ah_domain_error_if(node_to_id_.search(root_) == nullptr)
        << "Gen_Mo_On_Trees: root node not found in graph";
 
      // 2. Pre-collect adjacency lists as node IDs (Array for O(1) access)
      auto adj = Array<Array<size_t>>::create(n_);
      g_.for_each_node([&](Node * p)
        {
          const size_t pid = node_to_id_.find(p);
          DynList<size_t> neighbors;
          for (auto it = typename GT::Node_Arc_Iterator(p);
               it.has_curr(); it.next_ne())
            {
              auto * a = it.get_curr();
              auto * neighbor = g_.get_connected_node(a, p);
              neighbors.append(node_to_id_.find(neighbor));
            }
          adj(pid) = Array<size_t>(neighbors);
        });
 
      // 3. Iterative DFS to build Euler tours + parent/depth
      tin_       = Array<size_t>::create(n_);
      tout_      = Array<size_t>::create(n_);
      flat_sub_  = Array<T>::create(n_);
      first_     = Array<size_t>::create(n_);
      last_      = Array<size_t>::create(n_);
      flat_node_ = Array<size_t>::create(2 * n_);
      depth_     = Array<size_t>::create(n_);
      parent_    = Array<size_t>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        parent_(i) = NONE;
 
      auto visited = Array<bool>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        visited(i) = false;
 
      size_t sub_timer  = 0;
      size_t path_timer = 0;
 
      struct Frame { size_t id; size_t child_idx; };
      DynListStack<Frame> stk;
 
      const size_t root_id = node_to_id_.find(root_);
      visited(root_id) = true;
      depth_(root_id)  = 0;
      parent_(root_id) = NONE;
 
      tin_(root_id) = sub_timer;
      flat_sub_(sub_timer) = node_values_(root_id);
      ++sub_timer;
 
      first_(root_id) = path_timer;
      flat_node_(path_timer) = root_id;
      ++path_timer;
 
      stk.push({root_id, 0});
 
      while (not stk.is_empty())
        {
          auto & fr = stk.top();
          bool pushed = false;
 
          while (fr.child_idx < adj(fr.id).size())
            {
              const size_t nid = adj(fr.id)(fr.child_idx++);
              if (visited(nid))
                continue;
 
              visited(nid) = true;
              parent_(nid) = fr.id;
              depth_(nid)  = depth_(fr.id) + 1;
 
              tin_(nid) = sub_timer;
              flat_sub_(sub_timer) = node_values_(nid);
              ++sub_timer;
 
              first_(nid) = path_timer;
              flat_node_(path_timer) = nid;
              ++path_timer;
 
              stk.push({nid, 0});
              pushed = true;
              break;
            }
 
          if (not pushed)
            {
              tout_(fr.id) = sub_timer - 1;
 
              last_(fr.id) = path_timer;
              flat_node_(path_timer) = fr.id;
              ++path_timer;
 
              (void) stk.pop();
            }
        }
 
      ah_domain_error_if(sub_timer != n_)
        << "Gen_Mo_On_Trees: graph is not connected (not a tree)";
 
      // 4. Binary lifting table for LCA
      log_n_ = 1;
      while ((static_cast<size_t>(1) << log_n_) < n_)
        ++log_n_;
      ++log_n_;
 
      up_ = Array<size_t>::create(log_n_ * n_);
      for (size_t i = 0; i < log_n_ * n_; ++i)
        up_(i) = NONE;
 
      for (size_t v = 0; v < n_; ++v)
        up_(0 * n_ + v) = parent_(v);
 
      for (size_t k = 1; k < log_n_; ++k)
        for (size_t v = 0; v < n_; ++v)
          up_(k * n_ + v) = (up_((k-1) * n_ + v) == NONE)
            ? NONE
            : up_((k-1) * n_ + up_((k-1) * n_ + v));
    }
 
    // LCA via binary lifting – O(log N)
    size_t lca(size_t u, size_t v) const
    {
      if (depth_(u) < depth_(v))
        std::swap(u, v);
 
      const size_t diff = depth_(u) - depth_(v);
      for (size_t k = 0; k < log_n_; ++k)
        if ((diff >> k) & 1)
          u = up_(k * n_ + u);
 
      if (u == v)
        return u;
 
      for (int k = static_cast<int>(log_n_) - 1; k >= 0; --k)
        if (up_(k * n_ + u) != up_(k * n_ + v))
          {
            u = up_(k * n_ + u);
            v = up_(k * n_ + v);
          }
 
      return up_(0 * n_ + u);
    }
 
    // Standard Mo sweep (used for subtree queries)
    Array<answer_type> mo_sweep(const Array<T> & data,
                                Array<Mo_Query> queries,
                                size_t n) const
    {
      const size_t q = queries.size();
      const size_t block = std::max<size_t>(
        1, static_cast<size_t>(std::sqrt(static_cast<double>(n))));
 
      std::sort(&queries(0), &queries(0) + q,
                [block](const Mo_Query & a, const Mo_Query & b)
                  {
                    const size_t ba = a.l / block;
                    const size_t bb = b.l / block;
                    if (ba != bb)
                      return ba < bb;
                    return (ba & 1) ? (a.r > b.r) : (a.r < b.r);
                  });
 
      pol_.init(data, n);
 
      auto answers = Array<answer_type>::create(q);
 
      size_t cur_l = queries(0).l;
      size_t cur_r = queries(0).l;
      pol_.add(data, cur_l);
 
      while (cur_r < queries(0).r)
        pol_.add(data, ++cur_r);
 
      answers(queries(0).id) = pol_.answer();
 
      for (size_t i = 1; i < q; ++i)
        {
          const size_t ql = queries(i).l;
          const size_t qr = queries(i).r;
 
          while (cur_r < qr)
            pol_.add(data, ++cur_r);
 
          while (cur_l > ql)
            pol_.add(data, --cur_l);
 
          while (cur_r > qr)
            pol_.remove(data, cur_r--);
 
          while (cur_l < ql)
            pol_.remove(data, cur_l++);
 
          answers(queries(i).id) = pol_.answer();
        }
 
      return answers;
    }
 
  public:
 
    Gen_Mo_On_Trees(const GT & g, Node * root, Policy p = Policy())
      : g_(g), root_(root), pol_(std::move(p))
    {
      build();
    }
 
    [[nodiscard]] size_t size() const noexcept { return n_; }
 
    [[nodiscard]] bool is_empty() const noexcept { return n_ == 0; }
 
    // ── Subtree queries ──────────────────────────────────────────
 
    Array<answer_type> subtree_solve(
        const Array<Node*> & query_roots) const
    {
      const size_t q = query_roots.size();
      if (q == 0)
        return Array<answer_type>();
 
      auto queries = Array<Mo_Query>::create(q);
      for (size_t i = 0; i < q; ++i)
        {
          auto * p = node_to_id_.search(query_roots(i));
          ah_domain_error_if(p == nullptr)
            << "subtree_solve: query " << i << " node not in tree";
          const size_t id = p->second;
          queries(i) = {tin_(id), tout_(id), i};
        }
 
      return mo_sweep(flat_sub_, std::move(queries), n_);
    }
 
    Array<answer_type> subtree_solve(
        std::initializer_list<Node*> il) const
    {
      const Array<Node*> roots(il);
      return subtree_solve(roots);
    }
 
    // ── Path queries ─────────────────────────────────────────────
 
    Array<answer_type> path_solve(
        const Array<std::pair<Node*, Node*>> & query_pairs) const
    {
      const size_t q = query_pairs.size();
      if (q == 0)
        return Array<answer_type>();
 
      struct Path_Query
      {
        size_t l, r, id, lca_id;
      };
 
      auto pqueries = Array<Path_Query>::create(q);
      for (size_t i = 0; i < q; ++i)
        {
          auto * pu = node_to_id_.search(query_pairs(i).first);
          auto * pv = node_to_id_.search(query_pairs(i).second);
          ah_domain_error_if(pu == nullptr or pv == nullptr)
            << "path_solve: query " << i << " node not in tree";
 
          size_t u = pu->second;
          size_t v = pv->second;
 
          if (first_(u) > first_(v))
            std::swap(u, v);
 
          if (const size_t ancestor = lca(u, v); ancestor == u)
            pqueries(i) = {first_(u), first_(v), i, NONE};
          else
            pqueries(i) = {last_(u), first_(v), i, ancestor};
        }
 
      // Snake sort
      const size_t tour_sz = 2 * n_;
      const size_t block = std::max<size_t>(
        1, static_cast<size_t>(
          std::sqrt(static_cast<double>(tour_sz))));
 
      std::sort(&pqueries(0), &pqueries(0) + q,
                [block](const Path_Query & a, const Path_Query & b)
                  {
                    const size_t ba = a.l / block;
                    const size_t bb = b.l / block;
                    if (ba != bb)
                      return ba < bb;
                    return (ba & 1) ? (a.r > b.r) : (a.r < b.r);
                  });
 
      pol_.init(node_values_, n_);
 
      auto answers = Array<answer_type>::create(q);
      auto active = Array<bool>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        active(i) = false;
 
      // Toggle: flip node activity and call policy add/remove
      auto toggle = [&](const size_t pos)
        {
          if (const size_t nid = flat_node_(pos); active(nid))
            {
              pol_.remove(node_values_, nid);
              active(nid) = false;
            }
          else
            {
              pol_.add(node_values_, nid);
              active(nid) = true;
            }
        };
 
      // Initialize window to first query
      size_t cur_l = pqueries(0).l;
      size_t cur_r = pqueries(0).l;
      toggle(cur_l);
 
      while (cur_r < pqueries(0).r)
        toggle(++cur_r);
 
      if (pqueries(0).lca_id != NONE)
        {
          pol_.add(node_values_, pqueries(0).lca_id);
          answers(pqueries(0).id) = pol_.answer();
          pol_.remove(node_values_, pqueries(0).lca_id);
        }
      else
        answers(pqueries(0).id) = pol_.answer();
 
      // Sweep remaining queries
      for (size_t i = 1; i < q; ++i)
        {
          const size_t ql = pqueries(i).l;
          const size_t qr = pqueries(i).r;
 
          while (cur_r < qr) toggle(++cur_r);
          while (cur_l > ql) toggle(--cur_l);
          while (cur_r > qr) toggle(cur_r--);
          while (cur_l < ql) toggle(cur_l++);
 
          if (pqueries(i).lca_id != NONE)
            {
              pol_.add(node_values_, pqueries(i).lca_id);
              answers(pqueries(i).id) = pol_.answer();
              pol_.remove(node_values_, pqueries(i).lca_id);
            }
          else
            answers(pqueries(i).id) = pol_.answer();
        }
 
      return answers;
    }
 
    Array<answer_type> path_solve(
        std::initializer_list<std::pair<Node*, Node*>> il) const
    {
      const Array<std::pair<Node*, Node*>> pairs(il);
      return path_solve(pairs);
    }
  };
 
  // ================================================================
  // Gen_Mo_On_Tree_Node  —  Mo on N-ary trees (Tree_Node<T>)
  // ================================================================
 
  template <typename T, class Policy>
    requires MoPolicy<Policy, T>
  class Gen_Mo_On_Tree_Node
  {
  public:
    using Node = Tree_Node<T>;
    using answer_type = typename Policy::answer_type;
 
  private:
    static constexpr size_t NONE = ~size_t(0);
 
    Node * root_;
    mutable Policy pol_;
    size_t n_ = 0;
 
    // Node <-> ID mapping
    Array<Node*> id_to_node_;
    MapOLhash<Node*, size_t> node_to_id_;
    Array<T> node_values_;
 
    // Single-occurrence Euler tour (subtree queries, size n)
    Array<size_t> tin_;
    Array<size_t> tout_;
    Array<T> flat_sub_;
 
    // Double-occurrence Euler tour (path queries, size 2n)
    Array<size_t> first_;
    Array<size_t> last_;
    Array<size_t> flat_node_;
 
    // Binary lifting for LCA
    size_t log_n_ = 0;
    Array<size_t> depth_;
    Array<size_t> parent_;
    Array<size_t> up_;  // flattened [log_n_][n_]: up_(k * n_ + v)
 
    // ── Build ──────────────────────────────────────────────────────
 
    void build()
    {
      // 1. Collect nodes iteratively in preorder (no recursion).
      n_ = 0;
      if (root_ == nullptr)
        return;
 
      DynList<Node*> preorder;
      {
        DynListStack<Node*> stk;
        stk.push(root_);
 
        while (not stk.is_empty())
          {
            Node * curr = stk.pop();
            preorder.append(curr);
            ++n_;
 
            // Push children right-to-left so stack pops left-to-right.
            DynListStack<Node*> rev_children;
            for (Node * c = curr->get_left_child(); c != nullptr;
                 c = c->get_right_sibling())
              rev_children.push(c);
 
            while (not rev_children.is_empty())
              stk.push(rev_children.pop());
          }
      }
 
      if (n_ == 0)
        return;
 
      // 2. Assign IDs (pre-allocate OLhash for n_ entries)
      id_to_node_  = Array<Node*>::create(n_);
      node_values_ = Array<T>::create(n_);
      {
        MapOLhash<Node*, size_t> tmp(n_);
        node_to_id_.swap(tmp);
      }
 
      // Pre-collect adjacency as Array<Array<size_t>> for O(1) indexed access
      auto adj = Array<Array<size_t>>::create(n_);
      {
        size_t next_id = 0;
        for (auto it = preorder.get_it(); it.has_curr(); it.next_ne())
          {
            Node * mp = it.get_curr();
            node_to_id_.insert(mp, next_id);
            id_to_node_(next_id) = mp;
            node_values_(next_id) = mp->get_key();
            ++next_id;
          }
 
        // Build adjacency: children plus parent if it belongs to this rooted tree.
        for (size_t pid = 0; pid < n_; ++pid)
          {
            Node * mp = id_to_node_(pid);
            DynList<size_t> neighbors;
 
            for (Node * c = mp->get_left_child(); c != nullptr;
                 c = c->get_right_sibling())
              neighbors.append(node_to_id_.find(c));
 
            Node * par = mp->get_parent();
            if (par != nullptr)
              if (auto * pp = node_to_id_.search(par); pp != nullptr)
                neighbors.append(pp->second);
 
            adj(pid) = Array<size_t>(neighbors);
          }
      }
 
      // 3. Iterative DFS for Euler tours + parent/depth
      tin_       = Array<size_t>::create(n_);
      tout_      = Array<size_t>::create(n_);
      flat_sub_  = Array<T>::create(n_);
      first_     = Array<size_t>::create(n_);
      last_      = Array<size_t>::create(n_);
      flat_node_ = Array<size_t>::create(2 * n_);
      depth_     = Array<size_t>::create(n_);
      parent_    = Array<size_t>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        parent_(i) = NONE;
 
      auto visited = Array<bool>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        visited(i) = false;
 
      size_t sub_timer  = 0;
      size_t path_timer = 0;
 
      struct Frame { size_t id; size_t child_idx; };
      DynListStack<Frame> stk;
 
      const size_t root_id = node_to_id_.find(root_);
      visited(root_id) = true;
      depth_(root_id)  = 0;
      parent_(root_id) = NONE;
 
      tin_(root_id) = sub_timer;
      flat_sub_(sub_timer) = node_values_(root_id);
      ++sub_timer;
 
      first_(root_id) = path_timer;
      flat_node_(path_timer) = root_id;
      ++path_timer;
 
      stk.push({root_id, 0});
 
      while (not stk.is_empty())
        {
          auto & fr = stk.top();
          bool pushed = false;
 
          while (fr.child_idx < adj(fr.id).size())
            {
              const size_t nid = adj(fr.id)(fr.child_idx++);
              if (visited(nid))
                continue;
 
              visited(nid) = true;
              parent_(nid) = fr.id;
              depth_(nid)  = depth_(fr.id) + 1;
 
              tin_(nid) = sub_timer;
              flat_sub_(sub_timer) = node_values_(nid);
              ++sub_timer;
 
              first_(nid) = path_timer;
              flat_node_(path_timer) = nid;
              ++path_timer;
 
              stk.push({nid, 0});
              pushed = true;
              break;
            }
 
          if (not pushed)
            {
              tout_(fr.id) = sub_timer - 1;
 
              last_(fr.id) = path_timer;
              flat_node_(path_timer) = fr.id;
              ++path_timer;
 
              (void) stk.pop();
            }
        }
 
      ah_domain_error_if(sub_timer != n_)
        << "Gen_Mo_On_Tree_Node: tree traversal inconsistency";
 
      // 4. Binary lifting table for LCA
      log_n_ = 1;
      while ((static_cast<size_t>(1) << log_n_) < n_)
        ++log_n_;
      ++log_n_;
 
      up_ = Array<size_t>::create(log_n_ * n_);
      for (size_t i = 0; i < log_n_ * n_; ++i)
        up_(i) = NONE;
 
      for (size_t v = 0; v < n_; ++v)
        up_(0 * n_ + v) = parent_(v);
 
      for (size_t k = 1; k < log_n_; ++k)
        for (size_t v = 0; v < n_; ++v)
          up_(k * n_ + v) = (up_((k-1) * n_ + v) == NONE)
            ? NONE
            : up_((k-1) * n_ + up_((k-1) * n_ + v));
    }
 
    // LCA via binary lifting – O(log N)
    size_t lca(size_t u, size_t v) const
    {
      if (depth_(u) < depth_(v))
        std::swap(u, v);
 
      size_t diff = depth_(u) - depth_(v);
      for (size_t k = 0; k < log_n_; ++k)
        if ((diff >> k) & 1)
          u = up_(k * n_ + u);
 
      if (u == v)
        return u;
 
      for (int k = static_cast<int>(log_n_) - 1; k >= 0; --k)
        if (up_(k * n_ + u) != up_(k * n_ + v))
          {
            u = up_(k * n_ + u);
            v = up_(k * n_ + v);
          }
 
      return up_(0 * n_ + u);
    }
 
    // Standard Mo sweep
    Array<answer_type> mo_sweep(const Array<T> & data,
                                Array<Mo_Query> queries,
                                size_t nn) const
    {
      const size_t q = queries.size();
      const size_t block = std::max<size_t>(
        1, static_cast<size_t>(std::sqrt(static_cast<double>(nn))));
 
      std::sort(&queries(0), &queries(0) + q,
                [block](const Mo_Query & a, const Mo_Query & b)
                  {
                    const size_t ba = a.l / block;
                    const size_t bb = b.l / block;
                    if (ba != bb)
                      return ba < bb;
                    return (ba & 1) ? (a.r > b.r) : (a.r < b.r);
                  });
 
      pol_.init(data, nn);
 
      auto answers = Array<answer_type>::create(q);
 
      size_t cur_l = queries(0).l;
      size_t cur_r = queries(0).l;
      pol_.add(data, cur_l);
 
      while (cur_r < queries(0).r)
        pol_.add(data, ++cur_r);
 
      answers(queries(0).id) = pol_.answer();
 
      for (size_t i = 1; i < q; ++i)
        {
          const size_t ql = queries(i).l;
          const size_t qr = queries(i).r;
 
          while (cur_r < qr) pol_.add(data, ++cur_r);
          while (cur_l > ql) pol_.add(data, --cur_l);
          while (cur_r > qr) pol_.remove(data, cur_r--);
          while (cur_l < ql) pol_.remove(data, cur_l++);
 
          answers(queries(i).id) = pol_.answer();
        }
 
      return answers;
    }
 
  public:
 
    Gen_Mo_On_Tree_Node(Node * root, Policy p = Policy())
      : root_(root), pol_(std::move(p))
    {
      ah_domain_error_if(root == nullptr)
        << "Gen_Mo_On_Tree_Node: root is null";
      build();
    }
 
    [[nodiscard]] size_t size() const noexcept { return n_; }
 
    [[nodiscard]] bool is_empty() const noexcept { return n_ == 0; }
 
    // ── Subtree queries ──────────────────────────────────────────
 
    Array<answer_type> subtree_solve(
        const Array<Node*> & query_roots) const
    {
      const size_t q = query_roots.size();
      if (q == 0)
        return Array<answer_type>();
 
      auto queries = Array<Mo_Query>::create(q);
      for (size_t i = 0; i < q; ++i)
        {
          auto * p = node_to_id_.search(query_roots(i));
          ah_domain_error_if(p == nullptr)
            << "subtree_solve: query " << i << " node not in tree";
          const size_t id = p->second;
          queries(i) = {tin_(id), tout_(id), i};
        }
 
      return mo_sweep(flat_sub_, std::move(queries), n_);
    }
 
    Array<answer_type> subtree_solve(
        std::initializer_list<Node*> il) const
    {
      const Array<Node*> roots(il);
      return subtree_solve(roots);
    }
 
    // ── Path queries ─────────────────────────────────────────────
 
    Array<answer_type> path_solve(
        const Array<std::pair<Node*, Node*>> & query_pairs) const
    {
      const size_t q = query_pairs.size();
      if (q == 0)
        return Array<answer_type>();
 
      struct Path_Query
      {
        size_t l, r, id, lca_id;
      };
 
      auto pqueries = Array<Path_Query>::create(q);
      for (size_t i = 0; i < q; ++i)
        {
          auto * pu = node_to_id_.search(query_pairs(i).first);
          auto * pv = node_to_id_.search(query_pairs(i).second);
          ah_domain_error_if(pu == nullptr or pv == nullptr)
            << "path_solve: query " << i << " node not in tree";
 
          size_t u = pu->second;
          size_t v = pv->second;
 
          if (first_(u) > first_(v))
            std::swap(u, v);
 
          const size_t ancestor = lca(u, v);
 
          if (ancestor == u)
            pqueries(i) = {first_(u), first_(v), i, NONE};
          else
            pqueries(i) = {last_(u), first_(v), i, ancestor};
        }
 
      // Snake sort
      const size_t tour_sz = 2 * n_;
      const size_t block = std::max<size_t>(
        1, static_cast<size_t>(
          std::sqrt(static_cast<double>(tour_sz))));
 
      std::sort(&pqueries(0), &pqueries(0) + q,
                [block](const Path_Query & a, const Path_Query & b)
                  {
                    const size_t ba = a.l / block;
                    const size_t bb = b.l / block;
                    if (ba != bb)
                      return ba < bb;
                    return (ba & 1) ? (a.r > b.r) : (a.r < b.r);
                  });
 
      pol_.init(node_values_, n_);
 
      auto answers = Array<answer_type>::create(q);
      auto active = Array<bool>::create(n_);
      for (size_t i = 0; i < n_; ++i)
        active(i) = false;
 
      auto toggle = [&](const size_t pos)
        {
          if (const size_t nid = flat_node_(pos); active(nid))
            {
              pol_.remove(node_values_, nid);
              active(nid) = false;
            }
          else
            {
              pol_.add(node_values_, nid);
              active(nid) = true;
            }
        };
 
      size_t cur_l = pqueries(0).l;
      size_t cur_r = pqueries(0).l;
      toggle(cur_l);
 
      while (cur_r < pqueries(0).r)
        toggle(++cur_r);
 
      if (pqueries(0).lca_id != NONE)
        {
          pol_.add(node_values_, pqueries(0).lca_id);
          answers(pqueries(0).id) = pol_.answer();
          pol_.remove(node_values_, pqueries(0).lca_id);
        }
      else
        answers(pqueries(0).id) = pol_.answer();
 
      for (size_t i = 1; i < q; ++i)
        {
          const size_t ql = pqueries(i).l;
          const size_t qr = pqueries(i).r;
 
          while (cur_r < qr) toggle(++cur_r);
          while (cur_l > ql) toggle(--cur_l);
          while (cur_r > qr) toggle(cur_r--);
          while (cur_l < ql) toggle(cur_l++);
 
          if (pqueries(i).lca_id != NONE)
            {
              pol_.add(node_values_, pqueries(i).lca_id);
              answers(pqueries(i).id) = pol_.answer();
              pol_.remove(node_values_, pqueries(i).lca_id);
            }
          else
            answers(pqueries(i).id) = pol_.answer();
        }
 
      return answers;
    }
 
    Array<answer_type> path_solve(
        std::initializer_list<std::pair<Node*, Node*>> il) const
    {
      const Array<std::pair<Node*, Node*>> pairs(il);
      return path_solve(pairs);
    }
  };
 
  // ================================================================
  // Convenient typedefs
  // ================================================================
 
  template <class GT>
  using Distinct_Count_Mo_On_Trees =
    Gen_Mo_On_Trees<GT,
                    Distinct_Count_Policy<typename GT::Node::Node_Type>>;
 
  template <class GT>
  using Powerful_Array_Mo_On_Trees =
    Gen_Mo_On_Trees<GT,
                    Powerful_Array_Policy<typename GT::Node::Node_Type>>;
 
  template <class GT>
  using Range_Mode_Mo_On_Trees =
    Gen_Mo_On_Trees<GT,
                    Range_Mode_Policy<typename GT::Node::Node_Type>>;
 
  template <typename T>
  using Distinct_Count_Mo_On_Tree_Node =
    Gen_Mo_On_Tree_Node<T, Distinct_Count_Policy<T>>;
 
  template <typename T>
  using Powerful_Array_Mo_On_Tree_Node =
    Gen_Mo_On_Tree_Node<T, Powerful_Array_Policy<T>>;
 
  template <typename T>
  using Range_Mode_Mo_On_Tree_Node =
    Gen_Mo_On_Tree_Node<T, Range_Mode_Policy<T>>;
 
} // namespace Aleph
 
# endif // TPL_MO_ON_TREES_H