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 <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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline void merge(T *a, const long l, const long m, const long r, Array<T> &buf,
                  const Compare &cmp = Compare())
{
  const long s = r - l + 1;
 
  // Grow the buffer's logical size once instead of multiple push_back
  // calls. Array::putn() is additive (adds to the current logical size)
  // and only reallocates past the current capacity, so once `buf` has
  // grown to the largest `s` needed across the whole sort, every further
  // call here is O(1): just bumping the size counter back down would be
  // wrong (contents get fully overwritten below regardless), so the
  // buffer is deliberately never shrunk between merges.
  if (buf.size() < static_cast<size_t>(s))
    buf.putn(static_cast<size_t>(s) - buf.size());
 
  // 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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline void merge(T *a, const long l, const long m, const long r, const Compare &cmp = Compare())
{
  Array<T> buf(static_cast<size_t>(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, Array<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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline void mergesort(T *a, const long l, const long r, const Compare &cmp = Compare())
{
  if (l >= r)
    return;
 
  Array<T> buf(static_cast<size_t>(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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline void mergesort(DynList<T> &list, const Compare &cmp = Compare())
{
  mergesort<DynList<T>, Compare>(list, cmp);
}
 
template <typename T, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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;
 
      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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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;
    }
 
  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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> Compare = Aleph::less<T>>
void quicksort(DynList<T> &list, const Compare &cmp = Compare())
{
  quicksort<T, Compare>(static_cast<HTList &>(list), cmp);
}
 
template <typename T, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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 = std::size(ciura_gaps);
 
  // 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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;
    }
 
  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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline Array<const T *> build_index_ptr(const C<T> &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<const T *> ret = Array<const T *>::create(n);
  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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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;
}
 
 
namespace sort_utils_detail {
template <class C>
decltype(auto) indexed_value(const C &c, const size_t i)
{
  if constexpr (requires { c(i); })
    return c(i);
  else
    return c[i];
}
 
template <class C>
using Indexed_Value_Type =
  std::remove_cvref_t<decltype(indexed_value(std::declval<const C &>(), size_t{}))>;
 
template <class C, class Compare>
bool stable_index_less(const C &c, const size_t i, const size_t j, const Compare &cmp)
{
  auto &&lhs = indexed_value(c, i);
  auto &&rhs = indexed_value(c, j);
 
  if (cmp(lhs, rhs))
    return true;
 
  if (cmp(rhs, lhs))
    return false;
 
  return i < j;
}
}  // namespace sort_utils_detail
 
template <class C, StrictWeakOrder<sort_utils_detail::Indexed_Value_Type<C>> Compare = Aleph::less<sort_utils_detail::Indexed_Value_Type<C>>>
inline Array<size_t> build_index(const C &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<size_t> ret = Array<size_t>::create(n);
  for (size_t i = 0; i < n; ++i)
    ret(i) = i;
 
  quicksort_op(ret, [&a, &cmp](const size_t i, const size_t j)
    {
      return cmp(sort_utils_detail::indexed_value(a, i), sort_utils_detail::indexed_value(a, j));
    });
 
  return ret;
}
 
template <class C, StrictWeakOrder<sort_utils_detail::Indexed_Value_Type<C>> Compare = Aleph::less<sort_utils_detail::Indexed_Value_Type<C>>>
inline Array<size_t> stable_build_index(const C &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<size_t> ret = Array<size_t>::create(n);
  for (size_t i = 0; i < n; ++i)
    ret(i) = i;
 
  quicksort_op(ret, [&a, &cmp](const size_t i, const size_t j)
    {
      return sort_utils_detail::stable_index_less(a, i, j, cmp);
    });
 
  return ret;
}
 
template <template <typename> class C, typename T, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline Array<T *> build_index_ptr(C<T> &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<T *> ret = Array<T *>::create(n);
  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, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline Array<const T *> stable_build_index_ptr(const C<T> &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<size_t> idx = Array<size_t>::create(n);
  for (size_t i = 0; i < n; ++i)
    idx(i) = i;
 
  quicksort_op(idx, [&a, &cmp](const size_t i, const size_t j)
    {
      return sort_utils_detail::stable_index_less(a, i, j, cmp);
    });
 
  Array<const T *> ret = Array<const T *>::create(n);
  for (size_t i = 0; i < n; ++i)
    ret(i) = &a(idx(i));
 
  return ret;
}
 
template <template <typename> class C, typename T, StrictWeakOrder<T> Compare = Aleph::less<T>>
inline Array<T *> stable_build_index_ptr(C<T> &a, const Compare &cmp = Compare())
{
  const size_t n = a.size();
  Array<size_t> idx = Array<size_t>::create(n);
  for (size_t i = 0; i < n; ++i)
    idx(i) = i;
 
  quicksort_op(idx, [&a, &cmp](const size_t i, const size_t j)
    {
      return sort_utils_detail::stable_index_less(a, i, j, cmp);
    });
 
  Array<T *> ret = Array<T *>::create(n);
  for (size_t i = 0; i < n; ++i)
    ret(i) = &a(idx(i));
 
  return ret;
}
 
template <template <typename> class C, typename T, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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 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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 IntegerSortableValue<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 = IntegerSortableValue<T>;
 
template <typename T>
concept RadixSortable = IntegerSortableValue<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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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, StrictWeakOrder<T> 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);
}
}  // namespace Aleph
 
#endif  // TPL_SORT_UTILS_H