tpl_sort_utils.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_SORT_UTILS_H
#define TPL_SORT_UTILS_H
 
#include <ahUtils.H>
#include <ahFunctional.H>
#include <ah-concepts.H>
#include <ah-errors.H>
#include <tpl_arrayStack.H>
#include <tpl_array.H>
#include <tpl_dynArray.H>
#include <tpl_dynDlist.H>
#include <tpl_arrayHeap.H>
#include <htlist.H>
#include <limits>
#include <type_traits>
#include <vector>
#include <array>
#include <utility>
#include <cmath>
#include <algorithm>
 
namespace Aleph
{
  extern size_t Insertion_Threshold;
 
  extern size_t Quicksort_Threshold;
 
  extern const int Not_Found;
 
  template <template <typename> class Container, typename T,
            class Compare = Aleph::less<T>>
  [[nodiscard]] bool is_sorted(const Container<T> & cont, const Compare & cmp = Compare())
  {
    if (cont.is_empty())
      return true;
 
    typename Container<T>::Iterator it(cont);
    const T & first = it.get_curr();
    const T *prev = &first;
    it.next_ne();
    for (; it.has_curr(); it.next_ne())
      {
        const T & curr = it.get_curr();
        if (cmp(curr, *prev))
          return false;
        prev = &curr;
      }
    return true;
  }
 
  template <template <typename> class Container, typename T,
            class Compare = Aleph::less<T>>
  [[nodiscard]] std::pair<bool, size_t> search_inversion(const Container<T> & cont,
                                                         const Compare & cmp = Compare())
  {
    if (cont.is_empty())
      return std::make_pair(true, 0);
 
    typename Container<T>::Iterator it(cont);
    const T & first = it.get_curr();
    const T *prev = &first;
    size_t i = 1;
    it.next_ne();
    for (; it.has_curr(); it.next_ne(), ++i)
      {
        const T & curr = it.get_curr();
        if (cmp(curr, *prev))
          return std::make_pair(false, i);
        prev = &curr;
      }
 
    return std::make_pair(true, i);
  }
 
  template <template <typename> class Container, typename T,
            class Compare = Aleph::less<T>>
  [[nodiscard]] bool is_inversely_sorted(const Container<T> & cont, const Compare & cmp = Compare())
  {
    if (cont.is_empty())
      return true;
 
    typename Container<T>::Iterator it(cont);
    const T & first = it.get_curr();
    const T *prev = &first;
    it.next_ne();
    for (; it.has_curr(); it.next_ne())
      {
        const T & curr = it.get_curr();
        if (cmp(*prev, curr))
          return false;
        prev = &curr;
      }
    return true;
  }
 
  template <template <typename> class Container, typename T,
            class Compare = Aleph::less<T>>
  [[nodiscard]] std::pair<bool, size_t>
  test_sorted(const Container<T> & cont, const Compare & cmp = Compare())
  {
    if (cont.is_empty())
      return {true, 0};
 
    typename Container<T>::Iterator it(cont);
    const T & first = it.get_curr();
    const T *prev = &first;
    size_t i = 1;
    it.next_ne();
    for (; it.has_curr(); it.next_ne(), ++i)
      {
        const T & curr = it.get_curr();
        if (cmp(curr, *prev))
          return {false, i};
        prev = &curr;
      }
    return {true, i};
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void selection_sort(T *a, const size_t n, const Compare & cmp = Compare())
    noexcept(noexcept(cmp(a[0], a[0])) && std::is_nothrow_swappable_v<T>)
  {
    if (n < 2)
      return;
 
    for (size_t i = 0, min, j; i < n - 1; ++i)
      {
        for (min = i, j = i + 1; j < n; ++j)
          if (cmp(a[j], a[min]))
            min = j;
 
        if (cmp(a[min], a[i]))
          std::swap(a[min], a[i]);
      }
  }
 
  template <class Link, class Compare>
  inline
  Link * search_extreme(const Link & list, const Compare & cmp)
  {
    typename Link::Iterator it(const_cast<Link &>(list));
    Link *extreme = it.get_curr();
 
    for (it.next(); it.has_curr(); it.next_ne())
      {
        Link *curr = it.get_curr();
        if (cmp(curr, extreme))
          extreme = curr;
      }
 
    return extreme;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  const T * search_extreme(const DynList<T> & list, const Compare & cmp = Compare())
  {
    if (list.is_empty())
      return nullptr;
 
    typename DynList<T>::Iterator it(list);
    const T *extreme = &it.get_curr();
    it.next_ne();
    for (; it.has_curr(); it.next_ne())
      {
        const T & curr = it.get_curr();
        if (cmp(curr, *extreme))
          extreme = &curr;
      }
 
    return extreme;
  }
 
  template <class Compare>
  inline
  Dlink * search_extreme(const Dlink & list, const Compare & cmp = Compare())
  {
    return search_extreme<Dlink, Compare>(list, cmp);
  }
 
  template <class Compare>
  inline
  Slinknc * search_extreme(const Slinknc & list, const Compare & cmp = Compare())
  {
    return search_extreme<Slinknc, Compare>(list, cmp);
  }
 
  template <class Compare>
  inline
  void selection_sort(Dlink & list, Compare cmp)
  {
    Dlink aux;
    while (not list.is_empty())
      {
        Dlink *extreme = search_extreme<Dlink, Compare>(list, cmp);
        extreme->del(); // remove extreme from list
        aux.append(extreme); // insert it ordered into aux
      }
 
    list.swap(&aux);
  }
 
  template <typename Tlink,
            template <class> class Tnode,
            typename T, class Compare>
  class Compare_Tnode
  {
    Compare cmp;
 
  public:
    Compare_Tnode(Compare cmp_fct = Compare()) noexcept(std::is_nothrow_copy_constructible_v<Compare>)
      : cmp(cmp_fct)
    { /* empty */
    }
 
    bool operator ()(Tlink *l1, Tlink *l2) const
      noexcept(noexcept(std::declval<const Compare &>()(std::declval<const T &>(), std::declval<const T &>())))
    {
      auto *n1 = static_cast<Tnode<T> *>(l1);
      auto *n2 = static_cast<Tnode<T> *>(l2);
 
      assert(n1 == l1 and n2 == l2);
 
      return cmp(n1->get_data(), n2->get_data());
    }
 
    bool operator ()(Tlink *l, const T & x) const
      noexcept(noexcept(std::declval<const Compare &>()(std::declval<const T &>(), std::declval<const T &>())))
    {
      auto *n = static_cast<Tnode<T> *>(l);
 
      assert(n == l);
 
      return cmp(n->get_data(), x);
    }
  };
 
  template <typename T, class Compare>
  struct Compare_Dnode : public Compare_Tnode<Dlink, Dnode, T, Compare>
  {
    Compare_Dnode(const Compare & cmp = Compare())
      noexcept(std::is_nothrow_copy_constructible_v<Compare>)
      : Compare_Tnode<Dlink, Dnode, T, Compare>(cmp)
    { /* empty */
    }
  };
 
  template <typename T, class Compare>
  using Compare_Snodenc = Compare_Tnode<Slinknc, Snodenc, T, Compare>;
 
 
  template <typename T, class Compare = Aleph::less<T>>
  inline void selection_sort(Dnode<T> & list, Compare cmp)
  {
    selection_sort<Compare_Dnode<T, Compare>>
        (list, Compare_Dnode<T, Compare>(cmp));
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  [[nodiscard]] inline
  long sequential_search(T *a, const T & x, const long l, const long r,
                         Equal eq = Equal())
  {
    for (long i = l; i <= r; ++i)
      if (eq(a[i], x)) return i;
 
    return Not_Found;
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  long sequential_search(const DynArray<T> & a, const T & x,
                         const long l, const long r, Equal eq = Equal())
  {
    for (long i = l; i <= r; ++i)
      if (a.exist(i))
        if (eq(a(i), x))
          return i;
 
    return Not_Found;
  }
 
  template <class Link, typename T, class Equal>
  Link * sequential_search(const Link & list, const T & x, Equal & eq)
  {
    for (typename Link::Iterator it(const_cast<Link &>(list));
         it.has_curr(); it.next_ne())
      if (Link *curr = it.get_curr(); eq(curr, x))
        return curr;
 
    return nullptr;
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  Dlink * sequential_search(const Dlink & list, const T & x,
                            Equal eq = Equal())
  {
    Compare_Dnode<T, Equal> cmp(eq);
    return sequential_search<Dlink, T, Compare_Dnode<T, Equal>>(list, x, cmp);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  Slinknc * sequential_search(const Slinknc & list, const T & x, Equal eq = Equal())
  {
    Compare_Snodenc<T, Equal> cmp(eq);
    return sequential_search<Slinknc, T, Compare_Snodenc<T, Equal>>(list, x, cmp);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  Dnode<T> * sequential_search(Dnode<T> & list, const T & x, Equal & eq)
  {
    Compare_Dnode<T, Equal> cmp(eq);
    Dlink *ret = sequential_search<Dlink, T, Compare_Dnode<T, Equal>>(list, x, cmp);
 
    return ret == nullptr ? nullptr : static_cast<Dnode<T> *>(ret);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  Dnode<T> * sequential_search(Dnode<T> & list, const T & x, Equal && eq = Equal())
  {
    return sequential_search<T, Equal>(list, x, eq);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  const Dnode<T> * sequential_search(const Dnode<T> & list, const T & x, Equal & eq)
  {
    Compare_Dnode<T, Equal> cmp(eq);
    const auto & base = static_cast<const Dlink &>(list);
    const Dlink *ret = sequential_search<Dlink, T, Compare_Dnode<T, Equal>>(base, x, cmp);
 
    return ret == nullptr ? nullptr : static_cast<const Dnode<T> *>(ret);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  const Dnode<T> * sequential_search(const Dnode<T> & list, const T & x, Equal && eq = Equal())
  {
    return sequential_search<T, Equal>(list, x, eq);
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  T * sequential_search(DynDlist<T> & list, const T & x, Equal eq = Equal())
  {
    Dnode<T> *ret = sequential_search<T, Equal>(static_cast<Dnode<T> &>(list), x, eq);
    return ret != nullptr ? &ret->get_data() : nullptr;
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  const T * sequential_search(const DynDlist<T> & list, const T & x, Equal eq = Equal())
  {
    const Dnode<T> *ret = sequential_search<T, Equal>(static_cast<const Dnode<T> &>(list), x, eq);
    return ret != nullptr ? &ret->get_data() : nullptr;
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  T * sequential_search(DynList<T> & list, const T & x, Equal eq = Equal())
  {
    for (typename DynList<T>::Iterator it(list); it.has_curr(); it.next_ne())
      {
        T & curr = it.get_curr_ne();
        if (eq(curr, x))
          return &curr;
      }
    return nullptr;
  }
 
  template <typename T, class Equal = Aleph::equal_to<T>>
  inline
  const T * sequential_search(const DynList<T> & list, const T & x, Equal eq = Equal())
  {
    for (typename DynList<T>::Iterator it(list); it.has_curr(); it.next_ne())
      {
        const T & curr = it.get_curr_ne();
        if (eq(curr, x))
          return &curr;
      }
    return nullptr;
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long search_extreme(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    long extreme_index = l;
    for (long i = l + 1; i <= r; ++i)
      if (cmp(a[i], a[extreme_index])) // is there a new minimum?
        extreme_index = i; // yes
 
    return extreme_index;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long search_min(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(a, l, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::greater<T>>
  inline
  long search_max(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(a, l, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  Dnode<T> * search_extreme(const Dnode<T> & list, const Compare & cmp = Compare())
  {
    Compare_Dnode<T, Compare> cmp_dnode(cmp);
    Dlink *ret =
        search_extreme<Dlink, Compare_Dnode<T, Compare>>
        (const_cast<Dlink &>(static_cast<const Dlink &>(list)), cmp_dnode);
 
    return static_cast<Dnode<T> *>(ret);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  T * search_extreme(DynDlist<T> & list, const Compare & cmp = Compare())
  {
    const auto & base = static_cast<const Dnode<T> &>(list);
    Dnode<T> *ret = search_extreme<T, Compare>(const_cast<Dnode<T> &>(base), cmp);
 
    return ret != nullptr ? &(ret->get_data()) : nullptr;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  const T * search_extreme(const DynDlist<T> & list, const Compare & cmp = Compare())
  {
    const auto & base = static_cast<const Dnode<T> &>(list);
    Dnode<T> *ret = search_extreme<T, Compare>(const_cast<Dnode<T> &>(base), cmp);
 
    return ret != nullptr ? &(ret->get_data()) : nullptr;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  T * search_extreme(DynList<T> & list, const Compare & cmp = Compare())
  {
    if (list.is_empty())
      return nullptr;
 
    typename DynList<T>::Iterator it(list);
    T *extreme = &it.get_curr_ne();
    it.next_ne();
    for (; it.has_curr(); it.next_ne())
      {
        T & curr = it.get_curr_ne();
        if (cmp(curr, *extreme))
          extreme = &curr;
      }
 
    return extreme;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  T * search_min(DynDlist<T> & list, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  const T * search_min(const DynDlist<T> & list, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::greater<T>>
  inline
  T * search_max(DynDlist<T> & list, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::greater<T>>
  inline
  const T * search_max(const DynDlist<T> & list, const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void insertion_sort(T *a, const long l, const long r,
                      const Compare & cmp = Compare())
    noexcept(noexcept(cmp(std::declval<T&>(), std::declval<T&>())) and
             std::is_nothrow_move_constructible_v<T> and
             std::is_nothrow_move_assignable_v<T>)
  {
    if (l >= r)
      return;
    for (long i = l + 1, j; i <= r; ++i)
      {
        T tmp = std::move(a[i]); // store a[i], since it will be overwritten
        for (j = i; j > l and cmp(tmp, a[j - 1]); --j)
          a[j] = std::move(a[j - 1]); // shift to the right
 
        a[j] = std::move(tmp); // insert tmp into the gap
      }
  }
 
  template <class Compare>
  inline
  void insert_sorted(Dlink & list, Dlink *p, const Compare & cmp)
  {
    Dlink::Iterator it(list);
    while (it.has_curr() and cmp(it.get_curr(), p))
      it.next_ne();
 
    if (it.has_curr())
      it.get_curr()->append(p); // insert before current
    else
      list.append(p);
  }
 
  template <class Compare>
  inline
  void insert_sorted(HTList & list, Slinknc *p, const Compare & cmp)
  {
    // Handle empty list case
    if (list.is_empty())
      {
        list.insert(p);
        return;
      }
 
    if (Slinknc *first = list.get_first(); cmp(p, first) or not cmp(first, p)) // p <= first?
      {
        list.insert(p);
        return;
      }
 
    if (Slinknc *last = list.get_last(); cmp(last, p) or not cmp(p, last)) // p >= last?
      {
        list.append(p);
        return;
      }
 
    Slinknc *prev = list.get_first();
    HTList::Iterator it(list);
    for (it.next(); it.has_curr(); it.next_ne())
      {
        Slinknc *curr = it.get_curr();
        if (cmp(p, curr)) // p < curr
          {
            prev->insert(p);
            return;
          }
        prev = curr;
      }
    ah_logic_error()
        << "insert_sorted(): reached end of list traversal without inserting";
  }
 
  template <class ListType, class Compare>
  inline
  void list_insertion_sort(ListType & list, const Compare & cmp)
  {
    if (list.is_empty())
      return;
 
    ListType aux;
    aux.append(list.remove_first());
    while (not list.is_empty())
      insert_sorted<Compare>(aux, list.remove_first(), cmp);
 
    list.swap(aux);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void insertion_sort(DynList<T> & l, const Compare & cmp = Compare())
  {
    using Cmp = Compare_Snodenc<T, Compare>;
    Cmp c(cmp);
    list_insertion_sort<HTList, Cmp>(l, c);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  DynList<T> insertion_sort(DynList<T> && l, const Compare & cmp = Compare())
  {
    using Cmp = Compare_Snodenc<T, Compare>;
    Cmp c(cmp);
    list_insertion_sort<HTList, Cmp>(l, c);
    return std::move(l);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline void insertion_sort(Dnode<T> & list, const Compare & cmp = Compare())
  {
    using Cmp = Compare_Dnode<T, Compare>;
    Cmp c(cmp);
    list_insertion_sort<Dlink, Cmp>(list, c);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void merge(T *a, const long l, const long m, const long r,
             std::vector<T> & buf, const Compare & cmp = Compare())
  {
    const long s = r - l + 1;
 
    // Resize buffer once instead of multiple push_back calls
    // This avoids repeated capacity checks and potential reallocations
    if (buf.size() < static_cast<size_t>(s))
      buf.resize(static_cast<size_t>(s));
 
    // Copy left partition [l..m] to buffer start
    long buf_idx = 0;
    for (long i = l; i <= m; ++i)
      buf[buf_idx++] = std::move(a[i]);
 
    // Copy right partition [m+1..r] in reverse to buffer end
    buf_idx = s - 1;
    for (long j = m + 1; j <= r; ++j)
      buf[buf_idx--] = std::move(a[j]);
 
    // Merge from both ends toward the middle
    long i = 0;
    long j = s - 1;
 
    for (long k = l; k <= r; ++k)
      if (not cmp(buf[j], buf[i]))
        a[k] = std::move(buf[i++]);
      else
        a[k] = std::move(buf[j--]);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void merge(T *a, const long l, const long m, const long r,
             const Compare & cmp = Compare())
  {
    std::vector<T> buf;
    buf.reserve(r - l + 1);
    merge(a, l, m, r, buf, cmp);
  }
 
  template <typename T, class Compare>
  void mergesort(T *a, const long l, const long r, std::vector<T> & buf, Compare cmp)
  {
    if (l >= r)
      return;
 
    const long m = l + (r - l) / 2;
 
    mergesort<T, Compare>(a, l, m, buf, cmp);
    mergesort<T, Compare>(a, m + 1, r, buf, cmp);
 
    if (not cmp(a[m + 1], a[m]))
      return;
 
    merge<T, Compare>(a, l, m, r, buf, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void mergesort(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    if (l >= r)
      return;
 
    std::vector<T> buf;
    buf.reserve(r - l + 1);
    mergesort(a, l, r, buf, cmp);
  }
 
  template <typename Tlist, class Compare>
  inline
  void merge_lists(Tlist & l1, Tlist & l2, Tlist & result,
                   const Compare & cmp = Compare())
  {
    assert(result.is_empty());
 
    while (not l1.is_empty() and not l2.is_empty())
      if (cmp(l1.get_first_ne(), l2.get_first_ne()))
        result.append(l1.remove_first_ne());
      else
        result.append(l2.remove_first_ne());
 
    if (l1.is_empty())
      result.concat_list(l2);
    else
      result.concat_list(l1);
 
    assert(l1.is_empty() and l2.is_empty());
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void merge_lists(Dnode<T> & l1, Dnode<T> & l2, Dnode<T> & result,
                   const Compare & cmp = Compare())
  {
    merge_lists<Dnode<T>, Compare_Dnode<T, Compare>>(l1, l2, result, cmp);
  }
 
  template <typename Tlist, class Compare>
  inline
  void mergesort(Tlist & list, const Compare & cmp = Compare())
  {
    if (list.is_unitarian_or_empty())
      return;
 
    Tlist l, r;
    list.split_list_ne(l, r); // split into two lists
 
    mergesort<Tlist, Compare>(l, cmp);
    mergesort<Tlist, Compare>(r, cmp);
 
    merge_lists<Tlist, Compare>(l, r, list, cmp); // merge them
  }
 
  template <template <typename> class Tlist, typename T,
            class Compare = Aleph::less<T>>
  inline
  void mergeinsertsort(Tlist<T> & list, const Compare & cmp = Compare(),
                       const size_t lsz = Aleph::Insertion_Threshold)
  {
    if (list.is_unitarian_or_empty())
      return;
 
    if (list.size() <= lsz)
      {
        insertion_sort<T, Compare>(list, cmp);
        return;
      }
 
    Tlist<T> l, r;
    list.split_list(l, r); // split into two lists
 
    mergeinsertsort<Tlist, T, Compare>(l, cmp, lsz);
    mergeinsertsort<Tlist, T, Compare>(r, cmp, lsz);
 
    merge_lists<Tlist<T>, Compare>(l, r, list, cmp); // merge them
  }
 
  template <template <typename> class Tlist, typename T,
            class Compare = Aleph::less<T>>
  inline
  void mergesort(Tlist<T> & list, const Compare & cmp = Compare())
  {
    if (list.is_unitarian_or_empty())
      return;
 
    Tlist<T> l, r;
    list.split_list(l, r); // split into two lists
 
    mergesort<Tlist, T, Compare>(l, cmp);
    mergesort<Tlist, T, Compare>(r, cmp);
 
    merge_lists<Tlist<T>, Compare>(l, r, list, cmp); // merge them
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void mergesort(Dnode<T> & list, const Compare & cmp = Compare())
  {
    mergesort<Dnode<T>, Compare_Dnode<T, Compare>>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void mergesort(DynDlist<T> & list, const Compare & cmp = Compare())
  {
    mergesort<Dnode<T>, Compare_Dnode<T, Compare>>(list, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void mergesort(DynList<T> & list, const Compare & cmp = Compare())
  {
    mergesort<DynList<T>, Compare>(list, cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  DynList<T> mergesort(DynList<T> && list, const Compare & cmp = Compare())
  {
    mergesort<DynList<T>, Compare>(list, cmp);
    return std::move(list);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long select_pivot(T *a, long l, long r, const Compare & cmp = Compare())
    noexcept(noexcept(cmp(a[0], a[0])) && std::is_nothrow_swappable_v<T>);
 
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long partition(T *a, const long l, const long r, const Compare & cmp = Compare())
    noexcept(noexcept(cmp(a[0], a[0])) && std::is_nothrow_swappable_v<T>)
  {
    if (l >= r)
      return l;
 
    const long p = select_pivot<T, Compare>(a, l, r, cmp);
    std::swap(a[p], a[r]);
 
    // Note: pivot stays at a[r] throughout the loop - no copy needed
    // This saves a potentially expensive copy for large T types
    long i = l - 1; // index first element to the left > pivot
    long j = r; // index first element to the right < pivot
    while (true)
      {
        // advance while a[i] < pivot (pivot is at a[r])
        while (cmp(a[++i], a[r])) {}
 
        while (cmp(a[r], a[--j])) // advance while pivot < a[j]
          if (j == l) // Has the left edge been reached?
            break; // yes ==> the iteration must be finished
 
        if (i >= j)
          break;
 
        // At this point there is an inversion a[i] > pivot > a[j]
        std::swap(a[i], a[j]); // Eliminate inversion
      }
 
    std::swap(a[i], a[r]); // Place pivot in its final position
 
    return i;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort_rec(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    if (l >= r)
      return;
 
    const long pivot = partition<T, Compare>(a, l, r, cmp);
 
    quicksort_rec<T, Compare>(a, l, pivot - 1, cmp);
    quicksort_rec<T, Compare>(a, pivot + 1, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort_no_tail(T *a, long l, long r, const Compare & cmp = Compare())
  {
    while (l < r)
      {
        const long pivot = partition<T, Compare>(a, l, r, cmp);
        const long left_size = pivot - l;
        const long right_size = r - pivot;
 
        // Recurse on the smaller partition to ensure O(log n) stack depth.
        if (left_size < right_size)
          {
            quicksort_no_tail<T, Compare>(a, l, pivot - 1, cmp);
            l = pivot + 1; // tail recurse on the right by looping
          }
        else
          {
            quicksort_no_tail<T, Compare>(a, pivot + 1, r, cmp);
            r = pivot - 1; // tail recurse on the left by looping
          }
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort_rec_min(T *a, const long l, const long r,
                         const Compare & cmp = Compare())
  {
    if (r <= l)
      return;
 
    if (const long pivot = partition<T, Compare>(a, l, r, cmp); pivot - l < r - pivot)
    // which partition is smaller?
      { // left partition is smaller
        quicksort_rec_min<T, Compare>(a, l, pivot - 1, cmp);
        quicksort_rec_min<T, Compare>(a, pivot + 1, r, cmp);
      }
    else
      { // right partition is smaller
        quicksort_rec_min<T, Compare>(a, pivot + 1, r, cmp);
        quicksort_rec_min<T, Compare>(a, l, pivot - 1, cmp);
      }
  }
 
  template <typename T, class Compare>
  inline
  long select_pivot(T *a, const long l, const long r, const Compare & cmp)
    noexcept(noexcept(cmp(a[0], a[0])) && std::is_nothrow_swappable_v<T>)
  {
    const long m = l + (r - l) / 2;
 
    if (cmp(a[r], a[l]))
      std::swap(a[l], a[r]);
    if (cmp(a[m], a[l]))
      std::swap(a[m], a[l]);
    if (cmp(a[r], a[m]))
      std::swap(a[r], a[m]);
 
    return m;
  }
 
  inline size_t introsort_depth_limit(size_t n) noexcept;
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    if (r - l < Quicksort_Threshold)
      {
        insertion_sort(a, l, r, cmp);
        return;
      }
 
    const auto n = static_cast<size_t>(r - l + 1);
 
    typedef std::pair<long, long> Partition;
    FixedStack<Partition> stack(introsort_depth_limit(n) + 2);
    stack.push(Partition(l, r));
 
    while (stack.size() > 0)
      {
        auto [fst, snd] = stack.pop();
        const long diff = snd - fst;
 
        if (diff <= 0)
          continue;
 
        if (diff < Quicksort_Threshold)
          {
            insertion_sort(a, fst, snd, cmp);
            continue;
          }
 
        if (const long pivot = partition<T, Compare>(a, fst, snd, cmp); pivot - fst < snd - pivot)
        // which one is smaller?
          { // the left partition is smaller
            stack.push(Partition(pivot + 1, snd));
            stack.push(Partition(fst, pivot - 1));
          }
        else
          { // the right partition is smaller
            stack.push(Partition(fst, pivot - 1));
            stack.push(Partition(pivot + 1, snd));
          }
      }
  }
 
  inline size_t introsort_depth_limit(size_t n) noexcept
  {
    size_t depth = 0;
    while (n > 1)
      {
        n >>= 1;
        ++depth;
      }
    return depth * 2;
  }
 
  template <typename T, class Compare>
  void introsort_loop(T *a, long l, long r, size_t depth_limit,
                      const Compare & cmp)
  {
    while (r - l >= static_cast<long>(Quicksort_Threshold))
      {
        if (depth_limit == 0) [[unlikely]]
          {
            // Depth limit reached: switch to heapsort for guaranteed O(n log n)
            faster_heapsort(a + l, static_cast<size_t>(r - l + 1), cmp);
            return;
          }
 
        --depth_limit;
 
        // Median-of-three pivot selection
        // Recurse on smaller partition, iterate on larger (tail call elimination)
        if (const long pivot = partition<T, Compare>(a, l, r, cmp); pivot - l < r - pivot)
          {
            introsort_loop<T, Compare>(a, l, pivot - 1, depth_limit, cmp);
            l = pivot + 1;
          }
        else
          {
            introsort_loop<T, Compare>(a, pivot + 1, r, depth_limit, cmp);
            r = pivot - 1;
          }
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void introsort(T *a, const long l, const long r, const Compare & cmp = Compare())
  {
    if (r <= l)
      return;
 
    const auto n = static_cast<size_t>(r - l + 1);
 
    if (n < Quicksort_Threshold)
      {
        insertion_sort(a, l, r, cmp);
        return;
      }
 
    // Main introsort loop with depth limit
    introsort_loop<T, Compare>(a, l, r, introsort_depth_limit(n), cmp);
 
    // Final insertion sort pass for small unsorted subarrays
    insertion_sort(a, l, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void introsort(T *a, const size_t n, const Compare & cmp = Compare())
  {
    if (n <= 1)
      return;
    introsort(a, 0L, static_cast<long>(n - 1), cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void introsort(T *begin, T *end, const Compare & cmp = Compare())
  {
    if (begin >= end)
      return;
    const auto n = static_cast<size_t>(end - begin);
    introsort(begin, n, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void introsort(Array<T> & arr, const Compare & cmp = Compare())
  {
    const size_t n = arr.size();
    if (n <= 1)
      return;
    // Array uses contiguous memory via MemArray, so we can use pointer arithmetic
    introsort(&arr(0), n, cmp);
  }
 
  template <class Compare>
  void quicksort(Dlink & list, const Compare & cmp = Compare())
  {
    if (list.is_unitarian_or_empty())
      return;
 
    Dlink *pivot = list.remove_next();
    Dlink smaller, bigger; // lists of smaller and larger than pivot
 
    while (not list.is_empty())
      if (Dlink *p = list.remove_next(); cmp(p, pivot))
        smaller.append(p);
      else
        bigger.append(p);
 
    quicksort<Compare>(bigger, cmp);
    quicksort<Compare>(smaller, cmp);
 
    list.concat_list(&smaller); // restore sorted lists into list
    list.append(pivot);
    list.concat_list(&bigger);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  void quicksort(Dnode<T> & list, const Compare & cmp = Compare())
  {
    quicksort<Compare_Dnode<T, Compare>>(list, cmp);
  }
 
  template <typename T, class Compare>
  void quicksort(HTList & list, const Compare & cmp = Compare())
  {
    if (list.is_unitarian_or_empty())
      return;
 
    auto *pivot = static_cast<Snodenc<T> *>(list.remove_first_ne());
    HTList smaller, bigger; // lists of smaller and larger than pivot
 
    while (not list.is_empty())
      if (auto *p = static_cast<Snodenc<T> *>(list.remove_first_ne()); cmp(p->get_data(), pivot->get_data()))
        smaller.append(p);
      else
        bigger.append(p);
 
    quicksort<T, Compare>(bigger, cmp);
    quicksort<T, Compare>(smaller, cmp);
 
    list.concat_list(smaller); // restore sorted lists into list
    list.append(pivot);
    list.concat_list(bigger);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  void quicksort(DynList<T> & list, const Compare & cmp = Compare())
  {
    quicksort<T, Compare>(static_cast<HTList &>(list), cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort_insertion(T *a, const long l, const long r,
                           const Compare & cmp = Compare())
  {
    if (r <= l)
      return;
 
    const long pivot = partition<T, Compare>(a, l, r, cmp);
 
    const long l_size = pivot - l; // left partition size
    const long r_size = r - pivot; // right partition size
    bool left_done = false; // true if left partition is sorted
    bool right_done = false; // true if der partition is sorted
 
    if (l_size <= Aleph::Insertion_Threshold)
      { // sort left partition by insertion
        insertion_sort<T, Compare>(a, l, pivot - 1, cmp);
        left_done = true;
      }
 
    if (r_size <= Aleph::Insertion_Threshold)
      { // sort der partition by insertion
        insertion_sort<T, Compare>(a, pivot + 1, r, cmp);
        right_done = true;
      }
 
    if (left_done and right_done)
      return; //both partitions sorted by insertion
 
    if (left_done) // left partition sorted by insertion?
      { // Yeah; It only remains to recursively sort the right partition
        quicksort_insertion<T, Compare>(a, pivot + 1, r, cmp);
        return;
      }
 
    if (right_done) // der partition sorted by insertion?
      { // Yeah; It only remains to recursively sort the left partition
        quicksort_insertion<T, Compare>(a, l, pivot - 1, cmp);
        return;
      }
 
    // here, both partitions were not insertion sorted
    if (l_size < r_size) // sort smallest partition first
      { //smaller left partition
        quicksort_insertion<T, Compare>(a, l, pivot - 1, cmp);
        quicksort_insertion<T, Compare>(a, pivot + 1, r, cmp);
      }
    else
      { // smaller right partition
        quicksort_insertion<T, Compare>(a, pivot + 1, r, cmp);
        quicksort_insertion<T, Compare>(a, l, pivot - 1, cmp);
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long random_search(T *a, const T & x, const long l, const long r,
                     const Compare & cmp = Compare())
  {
    if (l > r)
      return Not_Found;
 
    const long pivot = partition<T, Compare>(a, l, r, cmp);
 
    if (cmp(x, a[pivot]))
      return random_search<T, Compare>(a, x, l, pivot - 1, cmp);
    if (cmp(a[pivot], x))
      return random_search<T, Compare>(a, x, pivot + 1, r, cmp);
 
    return pivot; // element found at index x
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long random_search(DynArray<T> & a, const T & x, const long l, const long r,
                     const Compare & cmp = Compare())
  {
    if (l > r)
      return Not_Found;
 
    const long pivot = partition(a, l, r, cmp);
 
    if (cmp(x, a(pivot)))
      return random_search<T, Compare>(a, x, l, pivot - 1, cmp);
    if (cmp(a(pivot), x))
      return random_search<T, Compare>(a, x, pivot + 1, r, cmp);
 
    return pivot; // element found at index x
  }
 
  template <typename T, class Compare>
  inline
  Dnode<T> * dlink_random_search(Dlink & list, const T & x,
                                 const Compare & cmp = Compare())
  {
    if (list.is_empty())
      return nullptr;
 
    Dnode<T> item(x);
    Dnode<T> *item_ptr = &item; // pointer to the cell containing x
 
    Dlink smaller; // list of those smaller than pivot
    Dlink bigger; // list of those greater than pivot
 
    auto *pivot = static_cast<Dnode<T> *>(list.remove_next());
 
    while (not list.is_empty())
      if (Dlink *p = list.remove_next(); cmp(p, pivot))
        smaller.append(p);
      else
        bigger.append(p);
 
    Dnode<T> *ret_val = nullptr;
    if (cmp(item_ptr, pivot))
      ret_val = dlink_random_search<T, Compare>(smaller, x, cmp);
    else if (cmp(pivot, item_ptr))
      ret_val = dlink_random_search<T, Compare>(bigger, x, cmp);
    else
      ret_val = pivot;
 
    assert(list.is_empty());
 
    list.swap(&smaller);
    list.append(pivot);
    list.concat_list(&bigger);
 
    return ret_val;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  Dnode<T> * random_search(Dlink & list, const T & x, const Compare & cmp = Compare())
  {
    return dlink_random_search<T, Compare_Dnode<T, Compare>>(list, x, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  T * random_search(DynDlist<T> & list, const T & x, const Compare & cmp = Compare())
  {
    Dnode<T> *p = dlink_random_search<T, Compare_Dnode<T, Compare>>(list, x, cmp);
 
    return p == nullptr ? nullptr : &(p->get_data());
  }
 
  template <typename T, class Compare>
  static inline
  std::pair<long, long>
  partition_three_way(T *a, long l, long r, const Compare & cmp)
  {
    const long pivot_idx = select_pivot<T, Compare>(a, l, r, cmp);
    std::swap(a[pivot_idx], a[r]);
 
    long lt = l;
    long gt = r - 1;
    long current = l;
 
    while (current <= gt)
      if (cmp(a[current], a[r]))
        {
          std::swap(a[lt], a[current]);
          lt++;
          current++;
        }
      else if (cmp(a[r], a[current]))
        {
          std::swap(a[current], a[gt]);
          gt--;
        }
      else
        current++;
 
    std::swap(a[current], a[r]);
    return {lt, current};
  }
 
  template <class Container, typename T, class Compare>
  static inline
  std::pair<long, long>
  partition_three_way_op(Container & a, long l, long r, const Compare & cmp)
  {
    const long pivot_idx = select_pivot_op<T, Compare>(a, l, r, cmp);
    std::swap(a(pivot_idx), a(r));
 
    long lt = l;
    long gt = r - 1;
    long current = l;
 
    while (current <= gt)
      if (cmp(a(current), a(r)))
        {
          std::swap(a(lt), a(current));
          lt++;
          current++;
        }
      else if (cmp(a(r), a(current)))
        {
          std::swap(a(current), a(gt));
          gt--;
        }
      else
        current++;
 
    std::swap(a(current), a(r));
    return {lt, current};
  }
 
  template <typename T, class Compare>
  static inline
  const T &__random_select(T *a, const long i, long l, long r,
                           const Compare & cmp)
  {
    while (l < r)
      {
        if (r - l < 10)
          {
            insertion_sort<T, Compare>(a, l, r, cmp);
            return a[i];
          }
 
        const auto [lt, current] = partition_three_way<T, Compare>(a, l, r, cmp);
 
        if (i >= lt and i <= current)
          return a[i];
 
        if (i < lt)
          r = lt - 1;
        else
          l = current + 1;
      }
 
    return a[l];
  }
 
  template <typename T, class Compare, typename Container>
  inline
  long select_pivot_op_impl(const Container & a, const long l, const long r,
                            const Compare & cmp)
  {
    assert(l <= r);
 
    if (r - l <= 5)
      return r;
 
    const long m = l + (r - l) / 2; // central element
 
    // Access array entries only once
    const T & la = a(l);
    const T & ra = a(r);
    const T & ma = a(m);
 
    const T *med_ptr = Aleph::median<T, Compare>(la, ma, ra, cmp);
 
    if (med_ptr == &la)
      return l;
 
    if (med_ptr == &ma)
      return m;
 
    assert(med_ptr == &ra);
 
    return r;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long select_pivot_op(const DynArray<T> & a, const long l, const long r,
                       const Compare & cmp = Compare())
  {
    return select_pivot_op_impl<T, Compare>(a, l, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long select_pivot_op(const Array<T> & a, const long l, const long r,
                       const Compare & cmp = Compare())
  {
    return select_pivot_op_impl<T, Compare>(a, l, r, cmp);
  }
 
  template <typename T, class Compare, typename Container>
  inline
  long partition_op_impl(Container & a, const long l, const long r,
                         const Compare & cmp)
  {
    if (l == r)
      return l;
 
    long i = l - 1;
    long j = r;
 
    // Move the selected pivot to position r
    if (const long pivot_idx = select_pivot_op<T, Compare>(a, l, r, cmp); pivot_idx != r)
      std::swap(a(pivot_idx), a(r));
 
    // Now pivot value is at a(r) - use direct access instead of reference
    while (true)
      {
        while (cmp(a(++i), a(r))) {}
 
        while (cmp(a(r), a(--j)))
          if (j == l)
            break;
 
        if (i >= j)
          break;
 
        std::swap(a(i), a(j));
      }
 
    std::swap(a(i), a(r));
 
    return i;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long partition_op(DynArray<T> & a, const long l, const long r,
                    const Compare & cmp = Compare())
  {
    return partition_op_impl<T, Compare>(a, l, r, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long partition_op(Array<T> & a, const long l, const long r,
                    const Compare & cmp = Compare())
  {
    return partition_op_impl<T, Compare>(a, l, r, cmp);
  }
 
  template <typename Container, class Compare>
  static inline
  void insertion_sort_range(Container & a, long l, long r, const Compare & cmp)
  {
    using T = typename Container::Item_Type;
    for (long p = l + 1; p <= r; ++p)
      {
        T key = std::move(a(p));
        long j = p - 1;
        while (j >= l and cmp(key, a(j)))
          {
            a(j + 1) = std::move(a(j));
            j--;
          }
        a(j + 1) = std::move(key);
      }
  }
 
  template <typename T, class Compare>
  static inline
  const T &__random_select(DynArray<T> & a, const long i,
                           long l, long r,
                           const Compare & cmp = Compare())
  {
    assert(i >= l and i <= r);
 
    while (l < r)
      {
        if (r - l < 10)
          {
            insertion_sort_range(a, l, r, cmp);
            return a(i);
          }
 
        const auto [lt, current] = partition_three_way_op<DynArray<T>, T, Compare>(a, l, r, cmp);
 
        if (i >= lt and i <= current)
          return a(i);
 
        if (i < lt)
          r = lt - 1;
        else
          l = current + 1;
      }
 
    return a(l);
  }
 
  template <typename T, class Compare>
  static inline
  const T &__random_select(Array<T> & a, const long i,
                           long l, long r,
                           const Compare & cmp = Compare())
  {
    assert(i >= l and i <= r);
 
    while (l < r)
      {
        if (r - l < 10)
          {
            insertion_sort_range(a, l, r, cmp);
            return a(i);
          }
 
        const auto [lt, current] = partition_three_way_op<Array<T>, T, Compare>(a, l, r, cmp);
 
        if (i >= lt and i <= current)
          return a(i);
 
        if (i < lt)
          r = lt - 1;
        else
          l = current + 1;
      }
 
    return a(l);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  const T &random_select(DynArray<T> & a, const long i, const Compare & cmp = Compare())
  {
    const auto n_sz = a.size();
    ah_out_of_range_error_if(n_sz == 0 or i < 0) << "index out of range";
 
    const long n = static_cast<long>(n_sz) - 1;
    ah_out_of_range_error_if(i > n) << "index out of range";
 
    return __random_select<T, Compare>(a, i, 0, n, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  const T &random_select(Array<T> & a, const long i, const Compare & cmp = Compare())
  {
    const auto n_sz = a.size();
    ah_out_of_range_error_if(n_sz == 0 or i < 0) << "index out of range";
 
    const long n = static_cast<long>(n_sz) - 1;
    ah_out_of_range_error_if(i > n) << "index out of range";
 
    return __random_select<T, Compare>(a, i, 0, n, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  const T &random_select(T *a, const long i, const long n,
                         const Compare & cmp = Compare())
  {
    ah_invalid_argument_if(a == nullptr and n > 0) << "null pointer";
 
    ah_out_of_range_error_if(n <= 0 or i < 0 or i >= n) << "index out of range";
 
    return __random_select<T, Compare>(a, i, 0, n - 1, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  const T &random_select(T *a, const long i, const long n,
                         Compare && cmp = Compare())
  {
    return random_select<T, Compare>(a, i, n, cmp);
  }
 
  template <class Compare>
  Dlink * dlink_random_select(Dlink & list, const size_t i,
                              const Compare & cmp = Compare())
  {
    if (list.is_empty())
      return nullptr;
 
    Dlink smaller; // list of minors than pivot
    Dlink bigger; // list of majors than pivot
 
    size_t smaller_count = 0, // number of elements in smaller
        bigger_count = 0; // number of elements in bigger
 
    Dlink *pivot = list.remove_next();
 
    while (not list.is_empty())
      if (Dlink *p = list.remove_next(); cmp(p, pivot)) // p < pivot?
        {
          smaller.append(p);
          ++smaller_count;
        }
      else
        {
          bigger.append(p);
          ++bigger_count;
        }
 
    if (i >= smaller_count + bigger_count + 1)
      {
        list.concat_list(&smaller);
        list.append(pivot);
        list.concat_list(&bigger);
        ah_out_of_range_error() << "index of selection greater than list's size";
      }
 
    Dlink *ret_val = nullptr;
    if (i == smaller_count)
      ret_val = pivot;
    else if (i < smaller_count)
      ret_val = dlink_random_select<Compare>(smaller, i, cmp);
    else
      ret_val = dlink_random_select<Compare>(bigger, i - (smaller_count + 1), cmp);
 
    list.concat_list(&smaller);
    list.append(pivot);
    list.concat_list(&bigger);
 
    return ret_val;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  Dnode<T> * random_select(Dlink & list, const size_t i, const Compare & cmp = Compare())
  {
    return static_cast<Dnode<T> *>(dlink_random_select<Compare_Dnode<T, Compare>>(list, i, cmp));
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  T * random_select(DynDlist<T> & list, const size_t i, const Compare & cmp = Compare())
  {
    Dlink *link = dlink_random_select<Compare_Dnode<T, Compare>>(list, i, cmp);
 
    auto *p = static_cast<Dnode<T> *>(link);
 
    return p != nullptr ? &(p->get_data()) : nullptr;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void selection_sort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    const long n = static_cast<long>(a.size());
    if (n < 2)
      return;
 
    for (long i = 0; i < n - 1; ++i)
      {
        long min = i;
 
        for (long j = i + 1; j < n; ++j)
          if (cmp(a(j), a(min)))
            min = j;
 
        if (cmp(a(min), a(i)))
          std::swap(a(min), a(i));
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void bubble_sort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    const long n = static_cast<long>(a.size());
    if (n < 2)
      return;
 
    for (long i = 0; i < n - 1; ++i)
      for (long j = n - 1; j > i; --j)
        if (cmp(a(j), a(j - 1)))
          std::swap(a(j - 1), a(j));
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  void insertion_sort(C<T> & a, const long l, const long r, const Compare & cmp = Compare())
  {
    for (long i = l + 1; i <= r; i++)
      {
        T tmp = std::move(a(i));
        long j = i;
        for (/* nothing */; j > l and cmp(tmp, a(j - 1)); --j)
          a(j) = std::move(a(j - 1));
 
        a(j) = std::move(tmp);
      }
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  void insertion_sort(C<T> & a, const Compare & cmp = Compare())
  {
    const auto n = a.size();
    if (n <= 1)
      return;
    insertion_sort(a, 0, static_cast<long>(n) - 1, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void shellsort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    const long n = a.size();
    if (n <= 1)
      return;
 
    // Ciura's gap sequence - empirically optimal for most data distributions
    // Extended with ~2.25x multiplier for larger arrays
    static constexpr long ciura_gaps[] =
        {
          1, 4, 10, 23, 57, 132, 301, 701, 1750, 4024, 9233,
          21223, 48805, 112217, 258100, 593630, 1365513, 3140680
        };
    static constexpr int num_gaps = sizeof(ciura_gaps) / sizeof(ciura_gaps[0]);
 
    // Find the first gap smaller than the array size
    int gap_idx = num_gaps - 1;
    while (gap_idx >= 0 && ciura_gaps[gap_idx] >= n)
      --gap_idx;
 
    // Sort with decreasing gaps
    for (; gap_idx >= 0; --gap_idx)
      {
        const long h = ciura_gaps[gap_idx];
        for (long i = h; i < n; i++)
          {
            T tmp = std::move(a(i));
            long j = i;
 
            while (j >= h and cmp(tmp, a(j - h)))
              {
                a(j) = std::move(a(j - h));
                j -= h;
              }
 
            a(j) = std::move(tmp);
          }
      }
  }
 
 
  inline static long back_index(const long i) noexcept { return i - 1; }
 
  template <typename T, class Compare>
  inline
  void sift_up(DynArray<T> & table, const long n, const Compare & cmp)
  {
    long p;
    for (long i = n; i > 1; i = p)
      {
        p = i >> 1; // c = i/2
 
        if (cmp(table(back_index(p)), table(back_index(i))))
          return;
 
        std::swap(table(back_index(p)), table(back_index(i)));
      }
  }
 
  template <typename T, class Compare>
  inline
  void sift_down(DynArray<T> & table, const long start, const long n, const Compare & cmp)
  {
    long i = start;
 
    while (true)
      {
        long c = i << 1; // c = 2*i (left child)
 
        if (c > n)
          return;
 
        // Select the smaller/larger child depending on comparator
        if (c + 1 <= n)
          if (cmp(table(back_index(c + 1)), table(back_index(c))))
            c++;
 
        if (cmp(table(back_index(i)), table(back_index(c))))
          return;
 
        std::swap(table(back_index(c)), table(back_index(i)));
        i = c;
      }
  }
 
  template <typename T, class Compare>
  inline
  void sift_down(DynArray<T> & table, const long n, const Compare & cmp)
  {
    sift_down<T, Compare>(table, 1, n, cmp);
  }
 
  template <typename T, class Compare>
  class Negate_Compare
  {
    Compare cmp;
 
  public:
    Negate_Compare(Compare __cmp = Compare())
      noexcept(std::is_nothrow_copy_assignable_v<Compare>)
      : cmp(__cmp)
    { /* Empty */
    }
 
    bool operator ()(const T & e1, const T & e2) const
      noexcept(noexcept(std::declval<const Compare &>()(e2, e1)))
    {
      return cmp(e2, e1);
    }
  };
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void heapsort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    const long n = a.size();
    if (n <= 1)
      return;
 
    auto nc = Negate_Compare<T, Compare>(cmp);
 
    // Floyd's heap construction: O(n) instead of O(n log n)
    // Start from the last non-leaf node and sift down each node
    for (long i = n / 2; i >= 1; --i)
      sift_down<T, Negate_Compare<T, Compare>>(a, i, n, nc);
 
    // Extract elements from the heap one by one
    for (long i = n; i >= 2; --i)
      {
        std::swap(a(0), a(i - 1));
        sift_down<T, Negate_Compare<T, Compare>>(a, 1, i - 1, nc);
      }
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long partition(C<T> & a, const long l, const long r, const Compare & cmp = Compare())
  {
    if (l == r)
      return l;
 
    long i = l - 1;
    long j = r;
    const T & pivot = a(r);
 
    while (true)
      {
        while (cmp(a(++i), pivot)) { /* Nothing else */ }
 
        while (cmp(pivot, a(--j)))
          if (j == l)
            break;
 
        if (i >= j)
          break;
 
        std::swap(a(i), a(j));
      }
 
    std::swap(a(i), a(r));
 
    return i;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort_rec(DynArray<T> & a, const long l, const long r,
                     const Compare & cmp = Compare())
  {
    if (r <= l)
      return;
 
    if (const long i = partition(a, l, r, cmp); i - l < r - i)
      {
        quicksort_rec<T, Compare>(a, l, i - 1, cmp);
        quicksort_rec<T, Compare>(a, i + 1, r, cmp);
      }
    else
      {
        quicksort_rec<T, Compare>(a, i + 1, r, cmp);
        quicksort_rec<T, Compare>(a, l, i - 1, cmp);
      }
  }
 
  template <class Stack, class A, class B>
  inline
  void push2(Stack & stack, const A & a, const B & b)
  {
    stack.push(b);
    stack.push(a);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void quicksort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    long l = 0, r = a.size() - 1;
 
    const size_t n = a.size();
    FixedStack<long> stack((introsort_depth_limit(n) + 2) * 2);
 
    push2(stack, l, r);
 
    while (not stack.is_empty())
      {
        l = stack.pop();
        r = stack.pop();
 
        if (r <= l)
          continue;
 
        if (const long i = partition(a, l, r, cmp); i - l > r - i)
          {
            push2(stack, l, i - 1);
            push2(stack, i + 1, r);
          }
        else
          {
            push2(stack, i + 1, r);
            push2(stack, l, i - 1);
          }
      }
  }
 
  template <typename T, class Compare>
  inline
  void sift_down_subrange(DynArray<T> & a, long start, const long n,
                          const long offset, const Compare & cmp)
  {
    while (true)
      {
        long c = start << 1; // c = 2*start (left child)
        if (c > n)
          return;
 
        // Select the child with higher priority
        if (c + 1 <= n)
          if (cmp(a(offset + c), a(offset + c - 1)))
            c++;
 
        if (cmp(a(offset + start - 1), a(offset + c - 1)))
          return;
 
        std::swap(a(offset + start - 1), a(offset + c - 1));
        start = c;
      }
  }
 
  template <typename T, class Compare>
  inline
  void heapsort_subrange(DynArray<T> & a, const long l, const long r,
                         const Compare & cmp)
  {
    const long n = r - l + 1;
    if (n <= 1)
      return;
 
    auto nc = Negate_Compare<T, Compare>(cmp);
 
    // Floyd's heap construction: O(n)
    for (long i = n / 2; i >= 1; --i)
      sift_down_subrange<T, Negate_Compare<T, Compare>>(a, i, n, l, nc);
 
    // Extract elements
    for (long i = n; i >= 2; --i)
      {
        std::swap(a(l), a(l + i - 1));
        sift_down_subrange<T, Negate_Compare<T, Compare>>(a, 1, i - 1, l, nc);
      }
  }
 
  template <typename T, class Compare>
  void introsort_loop(DynArray<T> & a, long l, long r, size_t depth_limit,
                      const Compare & cmp)
  {
    while (r - l >= static_cast<long>(Quicksort_Threshold))
      {
        if (depth_limit == 0) [[unlikely]]
          {
            // Depth limit reached: switch to heapsort for guaranteed O(n log n)
            heapsort_subrange<T, Compare>(a, l, r, cmp);
            return;
          }
 
        --depth_limit;
 
        // Partition and continue
 
        // Recurse on smaller partition, iterate on larger (tail call elimination)
        if (const long pivot = partition(a, l, r, cmp); pivot - l < r - pivot)
          {
            introsort_loop<T, Compare>(a, l, pivot - 1, depth_limit, cmp);
            l = pivot + 1;
          }
        else
          {
            introsort_loop<T, Compare>(a, pivot + 1, r, depth_limit, cmp);
            r = pivot - 1;
          }
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  void introsort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    const long n = a.size();
    if (n <= 1)
      return;
 
    if (n < static_cast<long>(Quicksort_Threshold))
      {
        insertion_sort(a, 0, n - 1, cmp);
        return;
      }
 
    // Main introsort loop with depth limit
    introsort_loop<T, Compare>(a, 0, n - 1, introsort_depth_limit(n), cmp);
 
    // Final insertion sort pass for small unsorted subarrays
    insertion_sort(a, 0, n - 1, cmp);
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long search_extreme(const DynArray<T> & a, const long l, const long r,
                      const Compare & cmp = Compare())
  {
    long extreme_index = l;
 
    for (long i = l + 1; i <= r; i++)
      if (cmp(a(i), a(extreme_index)))
        extreme_index = i;
 
    return extreme_index;
  }
 
  template <typename T, class Compare = Aleph::greater<T>>
  inline
  long search_max(const DynArray<T> & a, const long l, const long r,
                  const Compare & cmp = Compare())
  {
    return search_extreme<T, Compare>(a, l, r, cmp);
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long binary_search(const C<T> & a, const T & x, long l, long r,
                     const Compare & cmp = Compare())
  {
    if (l > r)
      return l;
 
    while (l <= r)
      if (long m = l + (r - l) / 2; cmp(x, a(m)))
        r = m - 1;
      else if (cmp(a(m), x))
        l = m + 1;
      else
        return m; // key found
 
    return l;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long binary_search(const C<T *> & a, const T & x, long l, long r,
                     const Compare & cmp = Compare())
  {
    if (l > r)
      return l;
 
    while (l <= r)
      if (long m = l + (r - l) / 2; cmp(x, *a(m)))
        r = m - 1;
      else if (cmp(*a(m), x))
        l = m + 1;
      else
        return m; // found key
 
    return l;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long binary_search(const C<T *> & a, const T & x, Compare && cmp = Compare())
  {
    const auto n = a.size();
    if (n == 0)
      return 0;
 
    return binary_search(a, x, 0, static_cast<long>(n) - 1, cmp);
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long binary_search(const C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    const auto n = a.size();
    if (n == 0)
      return 0;
 
    return binary_search(a, x, 0, static_cast<long>(n) - 1, cmp);
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<size_t> binary_search_dup(const C<T> & a, const T & x,
                                    const Compare & cmp = Compare())
  {
    DynList<size_t> ret;
    const auto n = a.size();
    if (n == 0)
      return ret;
 
    long idx = binary_search(a, x, 0, static_cast<long>(n) - 1, cmp);
    if (idx < 0)
      return ret;
 
    if (not are_equals(a(idx), x, cmp))
      return ret;
 
    ret.append(idx);
    for (long i = idx - 1; i >= 0; --i)
      {
        if (not are_equals(a(i), x, cmp))
          break;
        ret.insert(i);
      }
 
    for (long i = idx + 1; i < static_cast<long>(n); ++i)
      {
        if (not are_equals(a(i), x, cmp))
          break;
        ret.append(i);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynArray<const T *> build_index_ptr(const C<T> & a, const Compare & cmp = Compare())
  {
    const size_t & n = a.size();
    DynArray<const T *> ret;
    ret.reserve(a.size());
    for (size_t i = 0; i < n; ++i)
      ret(i) = &a(i);
 
    quicksort_op(ret, [&cmp](const T *ptr1, const T *ptr2)
                   {
                     return cmp(*ptr1, *ptr2);
                   });
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<const T *> bsearch_dup(const C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    DynList<const T *> ret;
    long idx = binary_search(a, x, cmp);
    if (idx < 0)
      return ret;
 
    const T *found_ptr = &a(idx);
    if (not are_equals(*found_ptr, x, cmp))
      return ret;
 
    for (long i = idx - 1; i >= 0; --i)
      {
        const T *ptr = &a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.insert(ptr);
      }
 
    ret.append(found_ptr);
 
    for (long i = idx + 1, n = a.size(); i < n; ++i)
      {
        const T *ptr = &a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.append(ptr);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  T * bsearch(C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    long i = binary_search(a, x, cmp);
    if (i < 0)
      return nullptr;
    T *ptr = &a(i);
    return are_equals(*ptr, x, cmp) ? ptr : nullptr;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  const T * bsearch(const C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    long i = binary_search(a, x, cmp);
    if (i < 0)
      return nullptr;
    const T *ptr = &a(i);
    return are_equals(*ptr, x, cmp) ? ptr : nullptr;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  T * bsearch(const C<T *> & a, const T & x, const Compare & cmp = Compare())
  {
    long i = binary_search(a, x, cmp);
    if (i < 0)
      return nullptr;
    T *ptr = a(i);
    return are_equals(*ptr, x, cmp) ? ptr : nullptr;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<T *> bsearch_dup(C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    DynList<T *> ret;
    long idx = binary_search(a, x, cmp);
    if (idx < 0)
      return ret;
 
    T *found_ptr = &a(idx);
    if (not are_equals(*found_ptr, x, cmp))
      return ret;
 
    for (long i = idx - 1; i >= 0; --i)
      {
        T *ptr = &a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.insert(ptr);
      }
 
    ret.append(found_ptr);
 
    for (long i = idx + 1, n = a.size(); i < n; ++i)
      {
        T *ptr = &a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.append(ptr);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<T *> bsearch_dup(const C<T *> & a, const T & x, const Compare & cmp = Compare())
  {
    DynList<T *> ret;
    const auto n = a.size();
    if (n == 0)
      return ret;
 
    long idx = binary_search(a, x, 0, static_cast<long>(n) - 1, cmp);
    if (idx < 0)
      return ret;
 
    T *found_ptr = a(idx);
    if (not are_equals(*found_ptr, x, cmp))
      return ret;
 
    for (long i = idx - 1; i >= 0; --i)
      {
        T *ptr = a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.insert(ptr);
      }
 
    ret.append(found_ptr);
 
    for (long i = idx + 1, n = a.size(); i < n; ++i)
      {
        T *ptr = a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.append(ptr);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  long binindex(const C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    return binary_search(a, x, cmp);
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<long> binindex_dup(const C<T> & a, const T & x, const Compare & cmp = Compare())
  {
    DynList<long> ret;
    long idx = binary_search(a, x, cmp);
    if (idx < 0)
      return ret;
 
    if (not are_equals(a(idx), x, cmp))
      return ret;
 
    const long mid = idx;
 
    for (long i = idx - 1; i >= 0; --i)
      {
        if (not are_equals(a(i), x, cmp))
          break;
        ret.insert(i);
      }
 
    ret.append(mid);
 
    for (long i = idx + 1, n = a.size(); i < n; ++i)
      {
        if (not are_equals(a(i), x, cmp))
          break;
        ret.append(i);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynList<long> binindex_dup(const C<T *> & a, const T & x, const Compare & cmp = Compare())
  {
    DynList<long> ret;
    long idx = binary_search(a, x, cmp);
    if (idx < 0)
      return ret;
 
    T *ptr = a(idx);
    if (not are_equals(*ptr, x, cmp))
      return ret;
 
    for (long i = idx - 1; i >= 0; --i)
      {
        ptr = a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.insert(i);
      }
 
    ret.append(idx);
 
    for (long i = idx + 1, n = a.size(); i < n; ++i)
      {
        T *ptr = a(i);
        if (not are_equals(*ptr, x, cmp))
          break;
        ret.append(i);
      }
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynArray<size_t> build_index(const C<T> & a, const Compare & cmp = Compare())
  {
    const size_t & n = a.size();
    DynArray<size_t> ret;
    ret.reserve(a.size());
    for (size_t i = 0; i < n; ++i)
      ret(i) = i;
 
    quicksort_op(ret,
                 [&a, &cmp](size_t i, size_t j) { return cmp(a(i), a(j)); });
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  DynArray<T *> build_index_ptr(C<T> & a, const Compare & cmp = Compare())
  {
    const size_t & n = a.size();
    DynArray<T *> ret;
    ret.reserve(a.size());
    for (size_t i = 0; i < n; ++i)
      ret(i) = &a(i);
 
    quicksort_op(ret, [&cmp](const T *ptr1, const T *ptr2)
                   {
                     return cmp(*ptr1, *ptr2);
                   });
 
    return ret;
  }
 
  template <template <typename> class C, typename T,
            class Compare = Aleph::less<T>>
  inline
  void quicksort_op(C<T> & a, const Compare & cmp = Compare(),
                    const size_t threshold = Quicksort_Threshold)
  {
    const size_t n = a.size();
    if (n <= 1)
      return;
 
    size_t l = 0, r = n - 1;
 
    FixedStack<size_t> stack((introsort_depth_limit(n) + 2) * 2);
 
    push2(stack, l, r);
 
    while (not stack.is_empty())
      {
        l = stack.pop();
        r = stack.pop();
 
        const size_t partition_size = r - l + 1;
        if (partition_size <= 1)
          continue;
 
        if (partition_size <= threshold)
          {
            insertion_sort(a, l, r, cmp);
            continue;
          }
 
        const auto i =
            static_cast<size_t>(partition_op<T, Compare>(a,
                                                         static_cast<long>(l),
                                                         static_cast<long>(r),
                                                         cmp));
 
        const size_t left_size = i > l ? i - l : 0;
        const size_t right_size = r > i ? r - i : 0;
 
        if (left_size > right_size)
          {
            if (left_size > 0)
              push2(stack, l, i - 1);
            if (right_size > 0)
              push2(stack, i + 1, r);
          }
        else
          {
            if (right_size > 0)
              push2(stack, i + 1, r);
            if (left_size > 0)
              push2(stack, l, i - 1);
          }
      }
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long binary_search_rec(T *a, const T & x, const long l, const long r,
                         const Compare & cmp = Compare())
  {
    if (l > r)
      return l;
 
    const long m = l + (r - l) / 2;
 
    if (cmp(x, a[m]))
      return binary_search_rec<T, Compare>(a, x, l, m - 1, cmp);
    if (cmp(a[m], x))
      return binary_search_rec<T, Compare>(a, x, m + 1, r, cmp);
 
    return m;
  }
 
  template <typename T, class Compare = Aleph::less<T>>
  inline
  long binary_search(T *a, const T & x, long l, long r,
                     const Compare & cmp = Compare())
  {
    while (l <= r)
      {
        long m = l + (r - l) / 2;
        if (cmp(x, a[m]))
          r = m - 1;
        else if (cmp(a[m], x))
          l = m + 1;
        else
          return m;
      }
    return l;
  }
 
 
  template <typename KeyFn>
    requires std::invocable<const KeyFn &, size_t> and
             std::convertible_to<std::invoke_result_t<const KeyFn &, size_t>, int>
  void counting_sort_indices(Array<size_t> & sa,
                             Array<size_t> & tmp,
                             const size_t n,
                             const int min_key,
                             const int max_key,
                             const KeyFn & key_of)
  {
    ah_domain_error_if(max_key < min_key)
      << "counting_sort_indices(): max_key < min_key";
 
    if (n == 0)
      return;
 
    ah_out_of_range_error_if(sa.size() < n or tmp.size() < n)
      << "counting_sort_indices(): sa/tmp size is smaller than n";
 
    const auto range_minus_one =
        static_cast<unsigned long long>(static_cast<long long>(max_key) -
                                        static_cast<long long>(min_key));
    constexpr auto max_count =
        static_cast<unsigned long long>(std::numeric_limits<size_t>::max() - size_t(1));
    ah_runtime_error_unless(range_minus_one <= max_count)
      << "counting_sort_indices(): key range is too large";
 
    const auto k = static_cast<size_t>(range_minus_one) + 1;
    Array<size_t> count = Array<size_t>::create(k);
    for (size_t i = 0; i < k; ++i)
      count[i] = 0;
 
    for (size_t i = 0; i < n; ++i)
      {
        const int key = key_of(sa[i]);
        ah_out_of_range_error_if(key < min_key or key > max_key)
          << "counting_sort_indices(): key out of expected range";
        const auto bucket =
            static_cast<size_t>(static_cast<unsigned long long>(
              static_cast<long long>(key) - static_cast<long long>(min_key)));
        ++count[bucket];
      }
 
    for (size_t j = 1; j < k; ++j)
      count[j] += count[j - 1];
 
    for (size_t i = n; i > 0; --i)
      {
        const int key = key_of(sa[i - 1]);
        const auto idx =
            static_cast<size_t>(static_cast<unsigned long long>(
              static_cast<long long>(key) - static_cast<long long>(min_key)));
        tmp[--count[idx]] = sa[i - 1];
      }
 
    for (size_t i = 0; i < n; ++i)
      sa[i] = tmp[i];
  }
 
 
  namespace sort_utils_detail
  {
    template <typename T>
    concept counting_integer = IntegerSortableValue<T>;
 
 
    template <typename T>
    concept radix_integer = IntegerSortableValue<T>;
 
 
    template <typename IntT>
    using counting_unsigned_t = std::make_unsigned_t<std::remove_cv_t<IntT>>;
 
 
    template <typename IntT>
    using radix_unsigned_t = std::make_unsigned_t<std::remove_cv_t<IntT>>;
 
 
    template <typename IntT>
    [[nodiscard]] size_t
    counting_bucket_count(const IntT min_key, const IntT max_key)
    {
      ah_domain_error_if(max_key < min_key) << "counting_sort(): max_key < min_key";
 
      using U = counting_unsigned_t<IntT>;
      const U range_minus_one =
          static_cast<U>(max_key) - static_cast<U>(min_key);
 
      if constexpr (sizeof(U) >= sizeof(size_t))
        {
          constexpr size_t max_count =
              std::numeric_limits<size_t>::max() - static_cast<size_t>(1);
          ah_runtime_error_unless(range_minus_one <= static_cast<U>(max_count))
            << "counting_sort(): key range is too large";
        }
 
      return static_cast<size_t>(range_minus_one) + 1;
    }
 
 
    template <typename IntT>
    [[nodiscard]] inline size_t
    counting_bucket(const IntT value, const IntT min_key) noexcept
    {
      return static_cast<size_t>(
        static_cast<counting_unsigned_t<IntT>>(value) -
        static_cast<counting_unsigned_t<IntT>>(min_key));
    }
 
 
    template <typename IntT>
    [[nodiscard]] inline IntT
    counting_value_from_bucket(const size_t bucket, const IntT min_key) noexcept
    {
      using U = counting_unsigned_t<IntT>;
      return static_cast<IntT>(static_cast<U>(min_key) + static_cast<U>(bucket));
    }
 
 
    template <template <typename> class C, typename IntT>
      requires counting_integer<IntT>
    void counting_sort_impl(C<IntT> & a)
    {
      const size_t n = a.size();
      if (n < 2)
        return;
 
      IntT min_key = a(0), max_key = a(0);
      for (size_t i = 1; i < n; ++i)
        {
          if (a(i) < min_key)
            min_key = a(i);
          if (a(i) > max_key)
            max_key = a(i);
        }
 
      const size_t k = counting_bucket_count(min_key, max_key);
      Array<size_t> count = Array<size_t>::create(k);
      for (size_t i = 0; i < k; ++i)
        count[i] = 0;
 
      for (size_t i = 0; i < n; ++i)
        ++count[counting_bucket(a(i), min_key)];
 
      for (size_t i = 1; i < k; ++i)
        count[i] += count[i - 1];
 
      Array<IntT> tmp = Array<IntT>::create(n);
      for (size_t i = n; i > 0; --i)
        {
          const auto bucket = counting_bucket(a(i - 1), min_key);
          tmp[--count[bucket]] = a(i - 1);
        }
 
      for (size_t i = 0; i < n; ++i)
        a(i) = tmp[i];
    }
 
 
    template <typename IntT>
      requires counting_integer<IntT>
    void counting_sort_impl(IntT *a, const size_t n)
    {
      if (n < 2)
        return;
 
      IntT min_key = a[0], max_key = a[0];
      for (size_t i = 1; i < n; ++i)
        {
          if (a[i] < min_key)
            min_key = a[i];
          if (a[i] > max_key)
            max_key = a[i];
        }
 
      const size_t k = counting_bucket_count(min_key, max_key);
      Array<size_t> count = Array<size_t>::create(k);
      for (size_t i = 0; i < k; ++i)
        count[i] = 0;
 
      for (size_t i = 0; i < n; ++i)
        ++count[counting_bucket(a[i], min_key)];
 
      for (size_t i = 1; i < k; ++i)
        count[i] += count[i - 1];
 
      Array<IntT> tmp = Array<IntT>::create(n);
      for (size_t i = n; i > 0; --i)
        {
          const auto bucket = counting_bucket(a[i - 1], min_key);
          tmp[--count[bucket]] = a[i - 1];
        }
 
      for (size_t i = 0; i < n; ++i)
        a[i] = tmp[i];
    }
 
 
    template <typename IntT>
      requires counting_integer<IntT>
    void counting_sort_impl(DynList<IntT> & list)
    {
      const size_t n = list.size();
      if (n < 2)
        return;
 
      typename DynList<IntT>::Iterator it(list);
      IntT min_key = it.get_curr_ne();
      IntT max_key = min_key;
      for (it.next_ne(); it.has_curr(); it.next_ne())
        {
          const IntT value = it.get_curr_ne();
          if (value < min_key)
            min_key = value;
          if (value > max_key)
            max_key = value;
        }
 
      const size_t k = counting_bucket_count(min_key, max_key);
      Array<size_t> count = Array<size_t>::create(k);
      for (size_t i = 0; i < k; ++i)
        count[i] = 0;
 
      for (typename DynList<IntT>::Iterator scan(list); scan.has_curr(); scan.next_ne())
        ++count[counting_bucket(scan.get_curr_ne(), min_key)];
 
      typename DynList<IntT>::Iterator out(list);
      for (size_t bucket = 0; bucket < k; ++bucket)
        {
          const size_t times = count[bucket];
          if (times == 0)
            continue;
 
          const IntT value = counting_value_from_bucket(bucket, min_key);
          for (size_t c = 0; c < times; ++c)
            {
              out.get_curr_ne() = value;
              out.next_ne();
            }
        }
    }
 
 
    template <typename IntT>
      requires counting_integer<IntT>
    void counting_sort_impl(DynDlist<IntT> & list)
    {
      const size_t n = list.size();
      if (n < 2)
        return;
 
      typename DynDlist<IntT>::Iterator it(list);
      IntT min_key = it.get_curr_ne();
      IntT max_key = min_key;
      for (it.next_ne(); it.has_curr(); it.next_ne())
        {
          const IntT value = it.get_curr_ne();
          if (value < min_key)
            min_key = value;
          if (value > max_key)
            max_key = value;
        }
 
      const size_t k = counting_bucket_count(min_key, max_key);
      Array<size_t> count = Array<size_t>::create(k);
      for (size_t i = 0; i < k; ++i)
        count[i] = 0;
 
      for (typename DynDlist<IntT>::Iterator scan(list); scan.has_curr(); scan.next_ne())
        ++count[counting_bucket(scan.get_curr_ne(), min_key)];
 
      typename DynDlist<IntT>::Iterator out(list);
      for (size_t bucket = 0; bucket < k; ++bucket)
        {
          const size_t times = count[bucket];
          if (times == 0)
            continue;
 
          const IntT value = counting_value_from_bucket(bucket, min_key);
          for (size_t c = 0; c < times; ++c)
            {
              out.get_curr_ne() = value;
              out.next_ne();
            }
        }
    }
 
 
    template <typename IntT>
    [[nodiscard]] constexpr radix_unsigned_t<IntT>
    radix_key(const IntT value) noexcept
    {
      using U = radix_unsigned_t<IntT>;
      if constexpr (std::is_signed_v<std::remove_cv_t<IntT>>)
        return static_cast<U>(value) ^
               (U{1} << (std::numeric_limits<U>::digits - 1));
      else
        return static_cast<U>(value);
    }
 
 
    template <template <typename> class C, typename IntT>
      requires radix_integer<IntT>
    void radix_sort_impl(C<IntT> & a)
    {
      using U = radix_unsigned_t<IntT>;
 
      const size_t n = a.size();
      if (n < 2)
        return;
 
      Array<IntT> tmp = Array<IntT>::create(n);
      Array<size_t> count = Array<size_t>::create(256);
 
      for (size_t pass = 0; pass < sizeof(U); ++pass)
        {
          for (size_t i = 0; i < 256; ++i)
            count[i] = 0;
 
          const size_t shift = pass * 8;
          for (size_t i = 0; i < n; ++i)
            {
              const U key = radix_key<IntT>(a(i));
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              ++count[bucket];
            }
 
          for (size_t i = 1; i < 256; ++i)
            count[i] += count[i - 1];
 
          for (size_t i = n; i > 0; --i)
            {
              const U key = radix_key<IntT>(a(i - 1));
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              tmp[--count[bucket]] = a(i - 1);
            }
 
          for (size_t i = 0; i < n; ++i)
            a(i) = tmp[i];
        }
    }
 
 
    template <typename IntT>
      requires radix_integer<IntT>
    void radix_sort_impl(IntT *a, const size_t n)
    {
      using U = radix_unsigned_t<IntT>;
 
      if (n < 2)
        return;
 
      Array<IntT> tmp = Array<IntT>::create(n);
      Array<size_t> count = Array<size_t>::create(256);
 
      for (size_t pass = 0; pass < sizeof(U); ++pass)
        {
          for (size_t i = 0; i < 256; ++i)
            count[i] = 0;
 
          const size_t shift = pass * 8;
          for (size_t i = 0; i < n; ++i)
            {
              const U key = radix_key<IntT>(a[i]);
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              ++count[bucket];
            }
 
          for (size_t i = 1; i < 256; ++i)
            count[i] += count[i - 1];
 
          for (size_t i = n; i > 0; --i)
            {
              const U key = radix_key<IntT>(a[i - 1]);
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              tmp[--count[bucket]] = a[i - 1];
            }
 
          for (size_t i = 0; i < n; ++i)
            a[i] = tmp[i];
        }
    }
 
 
    template <typename IntT>
      requires radix_integer<IntT>
    void radix_sort_impl(DynList<IntT> & list)
    {
      using U = radix_unsigned_t<IntT>;
 
      if (const size_t n = list.size(); n < 2)
        return;
 
      Array<HTList> buckets;
      buckets.reserve(256);
      for (size_t i = 0; i < 256; ++i)
        buckets.append(HTList());
 
      for (size_t pass = 0; pass < sizeof(U); ++pass)
        {
          const size_t shift = pass * 8;
 
          while (not list.is_empty())
            {
              Slinknc *link = list.HTList::remove_first_ne();
              auto *node = static_cast<Snodenc<IntT> *>(link);
              const U key = radix_key<IntT>(node->get_data());
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              buckets[bucket].append(link);
            }
 
          for (auto & bucket: buckets)
            list.HTList::append(bucket);
        }
    }
 
 
    template <typename IntT>
      requires radix_integer<IntT>
    void radix_sort_impl(DynDlist<IntT> & list)
    {
      using U = radix_unsigned_t<IntT>;
 
      if (const size_t n = list.size(); n < 2)
        return;
 
      Array<Dnode<IntT>> buckets;
      buckets.reserve(256);
      for (size_t i = 0; i < 256; ++i)
        buckets.append(Dnode<IntT>());
 
      for (size_t pass = 0; pass < sizeof(U); ++pass)
        {
          const size_t shift = pass * 8;
 
          while (not list.template Dnode<IntT>::is_empty())
            {
              Dnode<IntT> *node = list.template Dnode<IntT>::remove_first_ne();
              const U key = radix_key<IntT>(node->get_data());
              const auto bucket = static_cast<size_t>((key >> shift) & U{0xFF});
              buckets[bucket].append(node);
            }
 
          for (auto & bucket: buckets)
            list.template Dnode<IntT>::append_list(&bucket);
        }
    }
  } // namespace sort_utils_detail
 
 
  template <typename T>
  concept CountingSortable = sort_utils_detail::counting_integer<T>;
 
 
  template <typename T>
  concept RadixSortable = sort_utils_detail::radix_integer<T>;
 
 
  template <typename T>
    requires CountingSortable<T>
  void counting_sort(DynArray<T> & a)
  {
    sort_utils_detail::counting_sort_impl(a);
  }
 
 
  template <typename T>
    requires CountingSortable<T>
  void counting_sort(Array<T> & a)
  {
    sort_utils_detail::counting_sort_impl(a);
  }
 
 
  template <typename T>
    requires CountingSortable<T>
  void counting_sort(DynList<T> & list)
  {
    sort_utils_detail::counting_sort_impl(list);
  }
 
 
  template <typename T>
    requires CountingSortable<T>
  void counting_sort(DynDlist<T> & list)
  {
    sort_utils_detail::counting_sort_impl(list);
  }
 
 
  template <typename T>
    requires CountingSortable<T>
  void counting_sort(T *a, const size_t n)
  {
    ah_invalid_argument_if(a == nullptr and n > 0)
      << "counting_sort(): null pointer with non-zero length";
    sort_utils_detail::counting_sort_impl(a, n);
  }
 
 
  template <typename T, size_t N>
    requires CountingSortable<T>
  void counting_sort(T (&a)[N])
  {
    counting_sort(a, N);
  }
 
 
  template <typename T>
    requires RadixSortable<T>
  void radix_sort(DynArray<T> & a)
  {
    sort_utils_detail::radix_sort_impl(a);
  }
 
 
  template <typename T>
    requires RadixSortable<T>
  void radix_sort(Array<T> & a)
  {
    sort_utils_detail::radix_sort_impl(a);
  }
 
 
  template <typename T>
    requires RadixSortable<T>
  void radix_sort(DynList<T> & list)
  {
    sort_utils_detail::radix_sort_impl(list);
  }
 
 
  template <typename T>
    requires RadixSortable<T>
  void radix_sort(DynDlist<T> & list)
  {
    sort_utils_detail::radix_sort_impl(list);
  }
 
 
  template <typename T>
    requires RadixSortable<T>
  void radix_sort(T *a, const size_t n)
  {
    ah_invalid_argument_if(a == nullptr and n > 0)
      << "radix_sort(): null pointer with non-zero length";
    sort_utils_detail::radix_sort_impl(a, n);
  }
 
 
  template <typename T, size_t N>
    requires RadixSortable<T>
  void radix_sort(T (&a)[N])
  {
    radix_sort(a, N);
  }
 
 
  // ================================================================
  //                        Bucket Sort
  // ================================================================
 
  namespace bucket_sort_detail
  {
    template <typename T, class Compare>
    void insertion_sort_range(T *a, const size_t lo, const size_t hi,
                              const Compare & cmp)
      noexcept(std::is_nothrow_move_constructible_v<T> and
               std::is_nothrow_move_assignable_v<T> and
               noexcept(std::declval<const Compare &>()(std::declval<const T &>(),
                                                        std::declval<const T &>())))
    {
      for (size_t i = lo + 1; i < hi; ++i)
        {
          T tmp = std::move(a[i]);
          size_t j = i;
          while (j > lo and cmp(tmp, a[j - 1]))
            {
              a[j] = std::move(a[j - 1]);
              --j;
            }
          a[j] = std::move(tmp);
        }
    }
 
 
    template <typename T, class Compare, class BucketKey>
    void bucket_sort_impl(T *a, const size_t n, const size_t num_buckets,
                          const BucketKey & bucket_key, const Compare & cmp)
    {
      if (n < 2)
        return;
 
      ah_domain_error_if(num_buckets == 0)
        << "bucket_sort: num_buckets must be > 0";
 
      // Count elements per bucket
      Array<size_t> counts = Array<size_t>::create(num_buckets);
      for (size_t i = 0; i < num_buckets; ++i)
        counts[i] = 0;
 
      for (size_t i = 0; i < n; ++i)
        {
          const size_t b = bucket_key(a[i]);
          ah_domain_error_if(b >= num_buckets)
            << "bucket_sort: bucket_key returned " << b
            << " which is >= num_buckets (" << num_buckets << ")";
          ++counts[b];
        }
 
      // Compute bucket start offsets (prefix sum)
      Array<size_t> offsets = Array<size_t>::create(num_buckets);
      offsets[0] = 0;
      for (size_t i = 1; i < num_buckets; ++i)
        offsets[i] = offsets[i - 1] + counts[i - 1];
 
      // Distribute elements into a temporary array (stable)
      Array<T> tmp = Array<T>::create(n);
      // Reset counts for placement
      for (size_t i = 0; i < num_buckets; ++i)
        counts[i] = 0;
 
      for (size_t i = 0; i < n; ++i)
        {
          const size_t b = bucket_key(a[i]);
          tmp[offsets[b] + counts[b]] = std::move(a[i]);
          ++counts[b];
        }
 
      // Sort each bucket with insertion sort
      for (size_t b = 0; b < num_buckets; ++b)
        if (counts[b] > 1)
          insertion_sort_range(&tmp[0], offsets[b],
                               offsets[b] + counts[b], cmp);
 
      // Copy back
      for (size_t i = 0; i < n; ++i)
        a[i] = std::move(tmp[i]);
    }
 
    template <typename T, class Container, class SortFn, class Compare>
    void apply_sort_with_temp(Container & a, SortFn sort_fn, const Compare & cmp)
    {
      const size_t n = a.size();
      if (n < 2)
        return;
 
      Array<T> tmp = Array<T>::create(n);
      size_t i = 0;
      for (auto it = a.get_it(); it.has_curr(); it.next_ne())
        tmp[i++] = std::move(it.get_curr_ne());
 
      sort_fn(&tmp[0], n, cmp);
 
      i = 0;
      for (auto it = a.get_it(); it.has_curr(); it.next_ne())
        it.get_curr_ne() = std::move(tmp[i++]);
    }
 
    template <auto SortFn>
    struct sort_invoker_factory
    {
      template <class Compare>
      auto operator()(const Compare &) const
      {
        return [](auto *p, size_t n, const Compare & c) { SortFn(p, n, c); };
      }
    };
 
    template <class BucketKey>
    struct bucket_sort_invoker_factory
    {
      size_t num_buckets;
      const BucketKey & bucket_key;
 
      template <class Compare>
      auto operator()(const Compare &) const
      {
        return [this](auto *p, size_t n, const Compare & c)
               {
                 bucket_sort(p, n, num_buckets, bucket_key, c);
               };
      }
    };
  } // namespace bucket_sort_detail
 
 
  template <typename T, class Compare = Aleph::less<T>, class BucketKey>
  void bucket_sort(T *a, const size_t n, const size_t num_buckets,
                   const BucketKey & bucket_key,
                   const Compare & cmp = Compare())
  {
    ah_invalid_argument_if(a == nullptr and n > 0)
      << "bucket_sort(): null pointer with non-zero length";
    bucket_sort_detail::bucket_sort_impl(a, n, num_buckets, bucket_key, cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
    requires std::floating_point<T>
  void bucket_sort(T *a, const size_t n, const Compare & cmp = Compare())
  {
    ah_invalid_argument_if(a == nullptr and n > 0)
      << "bucket_sort(): null pointer with non-zero length";
 
    if (n < 2)
      return;
 
    T min_val = a[0], max_val = a[0];
    for (size_t i = 0; i < n; ++i)
      {
        if (not std::isfinite(a[i])) [[unlikely]]
          {
            timsort(a, n, cmp);
            return;
          }
        if (a[i] < min_val) min_val = a[i];
        if (a[i] > max_val) max_val = a[i];
      }
 
    const T range = max_val - min_val;
    if (range == T{0})
      return; // all equal
 
    // Numeric interpolation is only correct for the default arithmetic
    // comparators (Aleph::less / Aleph::greater).  For any other comparator
    // the mapping from value to bucket index is undefined, so fall back to
    // the stable timsort which respects the caller's comparator directly.
    if constexpr (std::is_same_v<Compare, Aleph::less<T>> or
                  std::is_same_v<Compare, Aleph::greater<T>>)
      {
        const bool descending = cmp(max_val, min_val) and not cmp(min_val, max_val);
        const size_t num_buckets = n;
        auto key = [min_val, range, num_buckets, descending](const T & val) -> size_t
          {
            auto b = static_cast<size_t>((val - min_val) / range *
                                         static_cast<T>(num_buckets - 1));
            if (b >= num_buckets)
              b = num_buckets - 1;
 
            if (descending)
              b = (num_buckets - 1) - b;
 
            return b;
          };
 
        bucket_sort_detail::bucket_sort_impl(a, n, num_buckets, key, cmp);
      }
    else
      {
        // Custom comparator: bucket interpolation is not order-preserving;
        // delegate to timsort which respects any strict-weak-ordering comparator.
        timsort(a, n, cmp);
      }
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
    requires std::floating_point<T>
  void bucket_sort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(a, 
      bucket_sort_detail::sort_invoker_factory<static_cast<void(*)(T*, size_t, const Compare&)>(bucket_sort)>{}(cmp),
      cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
    requires std::floating_point<T>
  void bucket_sort(Array<T> & a, const Compare & cmp = Compare())
  {
    if (a.size() < 2)
      return;
    bucket_sort(&a[0], a.size(), cmp);
  }
 
 
  template <typename T, size_t N, class Compare = Aleph::less<T>>
    requires std::floating_point<T>
  void bucket_sort(T (&a)[N], const Compare & cmp = Compare())
  {
    bucket_sort(static_cast<T *>(a), N, cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
    requires std::floating_point<T>
  void bucket_sort(DynList<T> & list, const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(list, 
      bucket_sort_detail::sort_invoker_factory<static_cast<void(*)(T*, size_t, const Compare&)>(bucket_sort)>{}(cmp),
      cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>, class BucketKey>
  void bucket_sort(DynArray<T> & a, const size_t num_buckets,
                   const BucketKey & bucket_key,
                   const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(a, 
      bucket_sort_detail::bucket_sort_invoker_factory<BucketKey>{num_buckets, bucket_key}(cmp),
      cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>, class BucketKey>
  void bucket_sort(Array<T> & a, const size_t num_buckets,
                   const BucketKey & bucket_key,
                   const Compare & cmp = Compare())
  {
    if (a.size() < 2)
      return;
    bucket_sort(&a[0], a.size(), num_buckets, bucket_key, cmp);
  }
 
 
  // ================================================================
  //                          Timsort
  // ================================================================
 
  namespace timsort_detail
  {
    [[nodiscard]] constexpr size_t compute_minrun(size_t n) noexcept
    {
      size_t r = 0;
      while (n >= 64)
        {
          r |= n & 1;
          n >>= 1;
        }
      return n + r;
    }
 
 
    template <typename T>
    void reverse_range(T *a, size_t lo, size_t hi)
      noexcept(std::is_nothrow_swappable_v<T>)
    {
      while (lo + 1 < hi)
        {
          std::swap(a[lo], a[hi - 1]);
          ++lo;
          --hi;
        }
    }
 
 
    template <typename T, class Compare>
    size_t count_run_and_make_ascending(T *a, const size_t lo,
                                        const size_t hi,
                                        const Compare & cmp)
    {
      if (hi - lo <= 1)
        return hi - lo;
 
      size_t run_hi = lo + 1;
      if (cmp(a[run_hi], a[lo])) // descending
        {
          while (run_hi < hi and cmp(a[run_hi], a[run_hi - 1]))
            ++run_hi;
          reverse_range(a, lo, run_hi);
        }
      else // ascending
        while (run_hi < hi and not cmp(a[run_hi], a[run_hi - 1]))
          ++run_hi;
      return run_hi - lo;
    }
 
 
    template <typename T, class Compare>
    void binary_insertion_sort(T *a, const size_t lo, const size_t hi,
                               const size_t start, const Compare & cmp)
      noexcept(std::is_nothrow_move_constructible_v<T> and
               std::is_nothrow_move_assignable_v<T> and
               noexcept(std::declval<const Compare &>()(std::declval<const T &>(),
                                                        std::declval<const T &>())))
 
    {
      size_t i = (start <= lo) ? lo + 1 : start;
      for (; i < hi; ++i)
        {
          T pivot = std::move(a[i]);
 
          // Binary search for insertion point in [lo, i)
          size_t left = lo;
          size_t right = i;
          while (left < right)
            {
              const size_t mid = left + (right - left) / 2;
              if (cmp(pivot, a[mid]))
                right = mid;
              else
                left = mid + 1;
            }
 
          // Shift elements [left, i) right by one
          for (size_t p = i; p > left; --p)
            a[p] = std::move(a[p - 1]);
          a[left] = std::move(pivot);
        }
    }
 
 
    template <typename T, class Compare>
    size_t gallop_left(const T & key, const T *a, const size_t len,
                       const size_t hint, const Compare & cmp)
      noexcept(noexcept(std::declval<const Compare &>()(std::declval<const T &>(),
                                                        std::declval<const T &>())))
    {
      size_t last_ofs = 0;
      size_t ofs = 1;
 
      if (cmp(a[hint], key)) // key > a[hint]: search right
        {
          const size_t max_ofs = len - hint;
          while (ofs < max_ofs and cmp(a[hint + ofs], key))
            {
              last_ofs = ofs;
              ofs = (ofs << 1) + 1;
              if (ofs > max_ofs)
                ofs = max_ofs;
            }
          last_ofs += hint;
          ofs += hint;
        }
      else // key <= a[hint]: search left
        {
          const size_t max_ofs = hint + 1;
          while (ofs < max_ofs and not cmp(a[hint - ofs], key))
            {
              last_ofs = ofs;
              ofs = (ofs << 1) + 1;
              if (ofs > max_ofs)
                ofs = max_ofs;
            }
          const size_t tmp = last_ofs;
          // Note: last_ofs may intentionally underflow if hint < ofs.
          // This is compensated by the subsequent ++last_ofs to effectively
          // start the binary search at index 0. This is safe for unsigned types.
          last_ofs = hint - ofs;
          ofs = hint - tmp;
        }
 
      // Binary search in (last_ofs, ofs]
      ++last_ofs;
      while (last_ofs < ofs)
        {
          const size_t mid = last_ofs + (ofs - last_ofs) / 2;
          if (cmp(a[mid], key))
            last_ofs = mid + 1;
          else
            ofs = mid;
        }
      return ofs;
    }
 
 
    template <typename T, class Compare>
    size_t gallop_right(const T & key, const T *a, const size_t len,
                        const size_t hint, const Compare & cmp)
      noexcept(noexcept(std::declval<const Compare &>()(std::declval<const T &>(),
                                                        std::declval<const T &>())))
    {
      size_t last_ofs = 0;
      size_t ofs = 1;
 
      if (cmp(key, a[hint])) // key < a[hint]: search left
        {
          const size_t max_ofs = hint + 1;
          while (ofs < max_ofs and cmp(key, a[hint - ofs]))
            {
              last_ofs = ofs;
              ofs = (ofs << 1) + 1;
              if (ofs > max_ofs)
                ofs = max_ofs;
            }
          const size_t tmp = last_ofs;
          // Note: last_ofs may intentionally underflow if hint < ofs.
          // This is compensated by the subsequent ++last_ofs to effectively
          // start the binary search at index 0. This is safe for unsigned types.
          last_ofs = hint - ofs;
          ofs = hint - tmp;
        }
      else // key >= a[hint]: search right
        {
          const size_t max_ofs = len - hint;
          while (ofs < max_ofs and not cmp(key, a[hint + ofs]))
            {
              last_ofs = ofs;
              ofs = (ofs << 1) + 1;
              if (ofs > max_ofs)
                ofs = max_ofs;
            }
          last_ofs += hint;
          ofs += hint;
        }
 
      ++last_ofs;
      while (last_ofs < ofs)
        {
          const size_t mid = last_ofs + (ofs - last_ofs) / 2;
          if (cmp(key, a[mid]))
            ofs = mid;
          else
            last_ofs = mid + 1;
        }
      return ofs;
    }
 
 
    template <typename T, class Compare>
    void merge_lo(T *a, const size_t base1, const size_t len1,
                  const size_t base2, const size_t len2,
                  T *tmp, const Compare & cmp)
    {
      // Copy run1 into tmp
      for (size_t i = 0; i < len1; ++i)
        tmp[i] = std::move(a[base1 + i]);
 
      size_t c1 = 0; // cursor into tmp (run1 copy)
      size_t c2 = base2; // cursor into a (run2)
      size_t dest = base1; // write cursor into a
      const size_t end1 = len1;
      const size_t end2 = base2 + len2;
 
      while (c1 < end1 and c2 < end2)
        if (cmp(a[c2], tmp[c1]))
          a[dest++] = std::move(a[c2++]);
        else
          a[dest++] = std::move(tmp[c1++]);
 
      // Copy remaining from tmp (run1)
      while (c1 < end1)
        a[dest++] = std::move(tmp[c1++]);
      // Remaining run2 elements are already in place
    }
 
 
    template <typename T, class Compare>
    void merge_hi(T *a, const size_t base1, const size_t len1,
                  const size_t base2, const size_t len2,
                  T *tmp, const Compare & cmp)
    {
      // Copy run2 into tmp
      for (size_t i = 0; i < len2; ++i)
        tmp[i] = std::move(a[base2 + i]);
 
      // Merge from the right (highest to lowest)
      size_t dest = base2 + len2; // one past end
      size_t c1 = base1 + len1; // one past end of run1 in a
      size_t c2 = len2; // one past end of run2 in tmp
 
      while (c1 > base1 and c2 > 0)
        if (cmp(tmp[c2 - 1], a[c1 - 1]))
          a[--dest] = std::move(a[--c1]);
        else
          a[--dest] = std::move(tmp[--c2]);
 
      // Copy remaining from tmp (run2)
      while (c2 > 0)
        a[--dest] = std::move(tmp[--c2]);
      // Remaining run1 elements are already in place
    }
 
 
    template <typename T, class Compare>
    void merge_at(T *a, size_t *run_base, size_t *run_len,
                  size_t & stack_size, const size_t at,
                  T *tmp, const Compare & cmp)
    {
      size_t base1 = run_base[at];
      size_t len1 = run_len[at];
      size_t base2 = run_base[at + 1];
      size_t len2 = run_len[at + 1];
 
      // Record merged run
      run_len[at] = len1 + len2;
      if (at == stack_size - 3)
        {
          run_base[at + 1] = run_base[at + 2];
          run_len[at + 1] = run_len[at + 2];
        }
      --stack_size;
 
      // Where does the first element of run2 go in run1?
      size_t k = gallop_right(a[base2], a + base1, len1, 0, cmp);
      base1 += k;
      len1 -= k;
      if (len1 == 0)
        return;
 
      // Where does the last element of run1 go in run2?
      len2 = gallop_left(a[base1 + len1 - 1], a + base2, len2,
                         len2 - 1, cmp);
      if (len2 == 0)
        return;
 
      if (len1 <= len2)
        merge_lo(a, base1, len1, base2, len2, tmp, cmp);
      else
        merge_hi(a, base1, len1, base2, len2, tmp, cmp);
    }
 
 
    template <typename T, class Compare>
    void merge_collapse(T *a, size_t *run_base, size_t *run_len,
                        size_t & stack_size, T *tmp,
                        const Compare & cmp)
    {
      while (stack_size > 1)
        {
          long n = static_cast<long>(stack_size) - 2;
          if ((n > 0 and run_len[n - 1] <= run_len[n] + run_len[n + 1]) or
              (n > 1 and run_len[n - 2] <= run_len[n - 1] + run_len[n]))
            {
              if (run_len[n - 1] < run_len[n + 1])
                --n;
              merge_at(a, run_base, run_len, stack_size, n, tmp, cmp);
            }
          else if (run_len[n] <= run_len[n + 1])
            merge_at(a, run_base, run_len, stack_size, n, tmp, cmp);
          else
            break;
        }
    }
 
 
    template <typename T, class Compare>
    void merge_force_collapse(T *a, size_t *run_base, size_t *run_len,
                              size_t & stack_size, T *tmp,
                              const Compare & cmp)
    {
      while (stack_size > 1)
        {
          long n = static_cast<long>(stack_size) - 2;
          if (n > 0 and run_len[n - 1] < run_len[n + 1])
            --n;
          merge_at(a, run_base, run_len, stack_size, n, tmp, cmp);
        }
    }
 
 
    template <typename T, class Compare>
    void timsort_impl(T *a, const size_t n, const Compare & cmp)
    {
      if (n < 2)
        return;
 
      // For very small arrays, just use binary insertion sort
      if (n < 64)
        {
          const size_t run_len = count_run_and_make_ascending(a, 0, n, cmp);
          binary_insertion_sort(a, 0, n, run_len, cmp);
          return;
        }
 
      const size_t minrun = compute_minrun(n);
 
      // Run stack — log2(n) + 1 entries is sufficient
      // (timsort invariant guarantees Fibonacci-like growth)
      constexpr size_t MAX_STACK = 85; // enough for 2^64 elements
      size_t run_base_arr[MAX_STACK];
      size_t run_len_arr[MAX_STACK];
      size_t stack_size = 0;
 
      // Temporary merge buffer — allocated once, sized to n/2
      // (merge_lo/merge_hi copy the smaller run)
      Array<T> tmp_buf = Array<T>::create(n / 2 + 1);
      T *tmp = &tmp_buf[0];
 
      size_t lo = 0;
      size_t remaining = n;
 
      while (remaining > 0)
        {
          // Find the next natural run
          size_t run_len = count_run_and_make_ascending(a, lo, lo + remaining,
                                                        cmp);
 
          // If the run is too short, extend it with insertion sort
          if (run_len < minrun)
            {
              const size_t force = remaining < minrun ? remaining : minrun;
              binary_insertion_sort(a, lo, lo + force, lo + run_len, cmp);
              run_len = force;
            }
 
          // Push this run onto the stack
          ah_runtime_error_unless(stack_size < MAX_STACK)
            << "timsort: run stack overflow (stack_size=" << stack_size
            << " >= MAX_STACK=" << MAX_STACK << "); this should never happen";
          run_base_arr[stack_size] = lo;
          run_len_arr[stack_size] = run_len;
          ++stack_size;
 
          // Maintain the stack invariants
          merge_collapse(a, run_base_arr, run_len_arr, stack_size, tmp, cmp);
 
          lo += run_len;
          remaining -= run_len;
        }
 
      // Force-merge all remaining runs
      merge_force_collapse(a, run_base_arr, run_len_arr, stack_size, tmp, cmp);
    }
  } // namespace timsort_detail
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(T *a, const size_t n, const Compare & cmp = Compare())
  {
    if (n == 0)
      return;
 
    ah_invalid_argument_if(a == nullptr)
      << "timsort(): null pointer with non-zero length";
 
    if (n < 2)
      return;
 
    timsort_detail::timsort_impl(a, n, cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(T *a, const long l, const long r,
               const Compare & cmp = Compare())
  {
    ah_range_error_if(l < 0 or r < 0)
      << "timsort(): negative bounds [" << l << ", " << r << "] are not allowed";
 
    if (l >= r)
      return;
 
    ah_invalid_argument_if(a == nullptr)
      << "timsort(): null pointer contract violation for range ["
      << l << ", " << r << "]; cannot call timsort_detail::timsort_impl()";
 
    timsort_detail::timsort_impl(a + l, static_cast<size_t>(r - l + 1), cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(DynArray<T> & a, const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(a, 
      bucket_sort_detail::sort_invoker_factory<static_cast<void(*)(T*, size_t, const Compare&)>(timsort)>{}(cmp),
      cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(Array<T> & a, const Compare & cmp = Compare())
  {
    if (a.size() < 2)
      return;
    timsort(&a[0], a.size(), cmp);
  }
 
 
  template <typename T, size_t N, class Compare = Aleph::less<T>>
    requires std::is_invocable_r_v<bool, Compare, const T &, const T &>
  void timsort(T (&a)[N], const Compare & cmp = Compare())
  {
    timsort(static_cast<T *>(a), N, cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(DynList<T> & list, const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(list, 
      bucket_sort_detail::sort_invoker_factory<static_cast<void(*)(T*, size_t, const Compare&)>(timsort)>{}(cmp),
      cmp);
  }
 
 
  template <typename T, class Compare = Aleph::less<T>>
  void timsort(DynDlist<T> & list, const Compare & cmp = Compare())
  {
    bucket_sort_detail::apply_sort_with_temp<T>(list, 
      bucket_sort_detail::sort_invoker_factory<static_cast<void(*)(T*, size_t, const Compare&)>(timsort)>{}(cmp),
      cmp);
  }
}
 
#endif // TPL_SORT_UTILS_H