RVec.hxx

// Author: Enrico Guiraud, Enric Tejedor, Danilo Piparo CERN  04/2021
// Implementation adapted from from llvm::SmallVector.
// See /math/vecops/ARCHITECTURE.md for more information.

/*************************************************************************
 * Copyright (C) 1995-2021, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT_RVEC
#define ROOT_RVEC

#if __cplusplus > 201402L
#define R__RVEC_NODISCARD [[nodiscard]]
#else
#define R__RVEC_NODISCARD
#endif

#ifdef _WIN32
   #ifndef M_PI
      #ifndef _USE_MATH_DEFINES
         #define _USE_MATH_DEFINES
      #endif
      #include <math.h>
      #undef _USE_MATH_DEFINES
   #endif
   #define _VECOPS_USE_EXTERN_TEMPLATES false
#else
   #define _VECOPS_USE_EXTERN_TEMPLATES true
#endif

#include <Rtypes.h> // R__CLING_PTRCHECK
#include <TError.h> // R__ASSERT

#include <algorithm>
#include <cmath>
#include <cstring>
#include <limits> // for numeric_limits
#include <memory> // uninitialized_value_construct
#include <new>
#include <numeric> // for inner_product
#include <sstream>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>

#ifdef R__HAS_VDT
#include <vdt/vdtMath.h>
#endif


namespace ROOT {

namespace VecOps {
template<typename T>
class RVec;
}

namespace Internal {
namespace VecOps {

template<typename T>
using RVec = ROOT::VecOps::RVec<T>;

// clang-format off
template <typename>
struct IsRVec : std::false_type {};

template <typename T>
struct IsRVec<ROOT::VecOps::RVec<T>> : std::true_type {};
// clang-format on

constexpr bool All(const bool *vals, std::size_t size)
{
   for (auto i = 0u; i < size; ++i)
      if (!vals[i])
         return false;
   return true;
}

template <typename... T>
std::size_t GetVectorsSize(const std::string &id, const RVec<T> &... vs)
{
   constexpr const auto nArgs = sizeof...(T);
   const std::size_t sizes[] = {vs.size()...};
   if (nArgs > 1) {
      for (auto i = 1UL; i < nArgs; i++) {
         if (sizes[0] == sizes[i])
            continue;
         std::string msg(id);
         msg += ": input RVec instances have different lengths!";
         throw std::runtime_error(msg);
      }
   }
   return sizes[0];
}

template <typename F, typename... RVecs>
auto MapImpl(F &&f, RVecs &&... vs) -> RVec<decltype(f(vs[0]...))>
{
   const auto size = GetVectorsSize("Map", vs...);
   RVec<decltype(f(vs[0]...))> ret(size);

   for (auto i = 0UL; i < size; i++)
      ret[i] = f(vs[i]...);

   return ret;
}

template <typename Tuple_t, std::size_t... Is>
auto MapFromTuple(Tuple_t &&t, std::index_sequence<Is...>)
   -> decltype(MapImpl(std::get<std::tuple_size<Tuple_t>::value - 1>(t), std::get<Is>(t)...))
{
   constexpr const auto tupleSizeM1 = std::tuple_size<Tuple_t>::value - 1;
   return MapImpl(std::get<tupleSizeM1>(t), std::get<Is>(t)...);
}

/// Return the next power of two (in 64-bits) that is strictly greater than A.
/// Return zero on overflow.
inline uint64_t NextPowerOf2(uint64_t A)
{
   A |= (A >> 1);
   A |= (A >> 2);
   A |= (A >> 4);
   A |= (A >> 8);
   A |= (A >> 16);
   A |= (A >> 32);
   return A + 1;
}

/// This is all the stuff common to all SmallVectors.
class R__CLING_PTRCHECK(off) SmallVectorBase {
public:
   // This limits the maximum size of an RVec<char> to ~4GB but we don't expect this to ever be a problem,
   // and we prefer the smaller Size_T to reduce the size of each RVec object.
   using Size_T = int32_t;

protected:
   void *fBeginX;
   /// Always >= 0.
   // Type is signed only for consistency with fCapacity.
   Size_T fSize = 0;
   /// Always >= -1. fCapacity == -1 indicates the RVec is in "memory adoption" mode.
   Size_T fCapacity;

   /// The maximum value of the Size_T used.
   static constexpr size_t SizeTypeMax() { return std::numeric_limits<Size_T>::max(); }

   SmallVectorBase() = delete;
   SmallVectorBase(void *FirstEl, size_t TotalCapacity) : fBeginX(FirstEl), fCapacity(TotalCapacity) {}

   /// This is an implementation of the grow() method which only works
   /// on POD-like data types and is out of line to reduce code duplication.
   /// This function will report a fatal error if it cannot increase capacity.
   void grow_pod(void *FirstEl, size_t MinSize, size_t TSize);

   /// Report that MinSize doesn't fit into this vector's size type. Throws
   /// std::length_error or calls report_fatal_error.
   static void report_size_overflow(size_t MinSize);
   /// Report that this vector is already at maximum capacity. Throws
   /// std::length_error or calls report_fatal_error.
   static void report_at_maximum_capacity();

   /// If false, the RVec is in "memory adoption" mode, i.e. it is acting as a view on a memory buffer it does not own.
   bool Owns() const { return fCapacity != -1; }

public:
   size_t size() const { return fSize; }
   size_t capacity() const noexcept { return Owns() ? fCapacity : fSize; }

   R__RVEC_NODISCARD bool empty() const { return !fSize; }

   /// Set the array size to \p N, which the current array must have enough
   /// capacity for.
   ///
   /// This does not construct or destroy any elements in the vector.
   ///
   /// Clients can use this in conjunction with capacity() to write past the end
   /// of the buffer when they know that more elements are available, and only
   /// update the size later. This avoids the cost of value initializing elements
   /// which will only be overwritten.
   void set_size(size_t N)
   {
      if (N > capacity()) {
         throw std::runtime_error("Setting size to a value greater than capacity.");
      }
      fSize = N;
   }
};

/// Used to figure out the offset of the first element of an RVec
template <class T>
struct SmallVectorAlignmentAndSize {
   alignas(SmallVectorBase) char Base[sizeof(SmallVectorBase)];
   alignas(T) char FirstEl[sizeof(T)];
};

/// This is the part of SmallVectorTemplateBase which does not depend on whether the type T is a POD.
template <typename T>
class R__CLING_PTRCHECK(off) SmallVectorTemplateCommon : public SmallVectorBase {
   using Base = SmallVectorBase;

   /// Find the address of the first element.  For this pointer math to be valid
   /// with small-size of 0 for T with lots of alignment, it's important that
   /// SmallVectorStorage is properly-aligned even for small-size of 0.
   void *getFirstEl() const
   {
      return const_cast<void *>(reinterpret_cast<const void *>(reinterpret_cast<const char *>(this) +
                                                               offsetof(SmallVectorAlignmentAndSize<T>, FirstEl)));
   }
   // Space after 'FirstEl' is clobbered, do not add any instance vars after it.

protected:
   SmallVectorTemplateCommon(size_t Size) : Base(getFirstEl(), Size) {}

   void grow_pod(size_t MinSize, size_t TSize) { Base::grow_pod(getFirstEl(), MinSize, TSize); }

   /// Return true if this is a smallvector which has not had dynamic
   /// memory allocated for it.
   bool isSmall() const { return this->fBeginX == getFirstEl(); }

   /// Put this vector in a state of being small.
   void resetToSmall()
   {
      this->fBeginX = getFirstEl();
      // from the original LLVM implementation:
      // FIXME: Setting fCapacity to 0 is suspect.
      this->fSize = this->fCapacity = 0;
   }

public:
   // note that fSize is a _signed_ integer, but we expose it as an unsigned integer for consistency with STL containers
   // as well as backward-compatibility
   using size_type = size_t;
   using difference_type = ptrdiff_t;
   using value_type = T;
   using iterator = T *;
   using const_iterator = const T *;

   using const_reverse_iterator = std::reverse_iterator<const_iterator>;
   using reverse_iterator = std::reverse_iterator<iterator>;

   using reference = T &;
   using const_reference = const T &;
   using pointer = T *;
   using const_pointer = const T *;

   using Base::capacity;
   using Base::empty;
   using Base::size;

   // forward iterator creation methods.
   iterator begin() noexcept { return (iterator)this->fBeginX; }
   const_iterator begin() const noexcept { return (const_iterator)this->fBeginX; }
   const_iterator cbegin() const noexcept { return (const_iterator)this->fBeginX; }
   iterator end() noexcept { return begin() + size(); }
   const_iterator end() const noexcept { return begin() + size(); }
   const_iterator cend() const noexcept { return begin() + size(); }

   // reverse iterator creation methods.
   reverse_iterator rbegin() noexcept { return reverse_iterator(end()); }
   const_reverse_iterator rbegin() const noexcept { return const_reverse_iterator(end()); }
   const_reverse_iterator crbegin() const noexcept { return const_reverse_iterator(end()); }
   reverse_iterator rend() noexcept { return reverse_iterator(begin()); }
   const_reverse_iterator rend() const noexcept { return const_reverse_iterator(begin()); }
   const_reverse_iterator crend() const noexcept { return const_reverse_iterator(begin()); }

   size_type size_in_bytes() const { return size() * sizeof(T); }
   size_type max_size() const noexcept { return std::min(this->SizeTypeMax(), size_type(-1) / sizeof(T)); }

   size_t capacity_in_bytes() const { return capacity() * sizeof(T); }

   /// Return a pointer to the vector's buffer, even if empty().
   pointer data() noexcept { return pointer(begin()); }
   /// Return a pointer to the vector's buffer, even if empty().
   const_pointer data() const noexcept { return const_pointer(begin()); }

   reference front()
   {
      if (empty()) {
         throw std::runtime_error("`front` called on an empty RVec");
      }
      return begin()[0];
   }

   const_reference front() const
   {
      if (empty()) {
         throw std::runtime_error("`front` called on an empty RVec");
      }
      return begin()[0];
   }

   reference back()
   {
      if (empty()) {
         throw std::runtime_error("`back` called on an empty RVec");
      }
      return end()[-1];
   }

   const_reference back() const
   {
      if (empty()) {
         throw std::runtime_error("`back` called on an empty RVec");
      }
      return end()[-1];
   }
};

/// SmallVectorTemplateBase<TriviallyCopyable = false> - This is where we put
/// method implementations that are designed to work with non-trivial T's.
///
/// We approximate is_trivially_copyable with trivial move/copy construction and
/// trivial destruction. While the standard doesn't specify that you're allowed
/// copy these types with memcpy, there is no way for the type to observe this.
/// This catches the important case of std::pair<POD, POD>, which is not
/// trivially assignable.
template <typename T, bool = (std::is_trivially_copy_constructible<T>::value) &&
                             (std::is_trivially_move_constructible<T>::value) &&
                             std::is_trivially_destructible<T>::value>
class R__CLING_PTRCHECK(off) SmallVectorTemplateBase : public SmallVectorTemplateCommon<T> {
protected:
   SmallVectorTemplateBase(size_t Size) : SmallVectorTemplateCommon<T>(Size) {}

   static void destroy_range(T *S, T *E)
   {
      while (S != E) {
         --E;
         E->~T();
      }
   }

   /// Move the range [I, E) into the uninitialized memory starting with "Dest",
   /// constructing elements as needed.
   template <typename It1, typename It2>
   static void uninitialized_move(It1 I, It1 E, It2 Dest)
   {
      std::uninitialized_copy(std::make_move_iterator(I), std::make_move_iterator(E), Dest);
   }

   /// Copy the range [I, E) onto the uninitialized memory starting with "Dest",
   /// constructing elements as needed.
   template <typename It1, typename It2>
   static void uninitialized_copy(It1 I, It1 E, It2 Dest)
   {
      std::uninitialized_copy(I, E, Dest);
   }

   /// Grow the allocated memory (without initializing new elements), doubling
   /// the size of the allocated memory. Guarantees space for at least one more
   /// element, or MinSize more elements if specified.
   void grow(size_t MinSize = 0);

public:
   void push_back(const T &Elt)
   {
      if (R__unlikely(this->size() >= this->capacity()))
         this->grow();
      ::new ((void *)this->end()) T(Elt);
      this->set_size(this->size() + 1);
   }

   void push_back(T &&Elt)
   {
      if (R__unlikely(this->size() >= this->capacity()))
         this->grow();
      ::new ((void *)this->end()) T(::std::move(Elt));
      this->set_size(this->size() + 1);
   }

   void pop_back()
   {
      this->set_size(this->size() - 1);
      this->end()->~T();
   }
};

// Define this out-of-line to dissuade the C++ compiler from inlining it.
template <typename T, bool TriviallyCopyable>
void R__CLING_PTRCHECK(off) SmallVectorTemplateBase<T, TriviallyCopyable>::grow(size_t MinSize)
{
   // Ensure we can fit the new capacity.
   // This is only going to be applicable when the capacity is 32 bit.
   if (MinSize > this->SizeTypeMax())
      this->report_size_overflow(MinSize);

   // Ensure we can meet the guarantee of space for at least one more element.
   // The above check alone will not catch the case where grow is called with a
   // default MinSize of 0, but the current capacity cannot be increased.
   // This is only going to be applicable when the capacity is 32 bit.
   if (this->capacity() == this->SizeTypeMax())
      this->report_at_maximum_capacity();

   // Always grow, even from zero.
   size_t NewCapacity = size_t(NextPowerOf2(this->capacity() + 2));
   NewCapacity = std::min(std::max(NewCapacity, MinSize), this->SizeTypeMax());
   T *NewElts = static_cast<T *>(malloc(NewCapacity * sizeof(T)));
   R__ASSERT(NewElts != nullptr);

   // Move the elements over.
   this->uninitialized_move(this->begin(), this->end(), NewElts);

   if (this->Owns()) {
      // Destroy the original elements.
      destroy_range(this->begin(), this->end());

      // If this wasn't grown from the inline copy, deallocate the old space.
      if (!this->isSmall())
         free(this->begin());
   }

   this->fBeginX = NewElts;
   this->fCapacity = NewCapacity;
}

/// SmallVectorTemplateBase<TriviallyCopyable = true> - This is where we put
/// method implementations that are designed to work with trivially copyable
/// T's. This allows using memcpy in place of copy/move construction and
/// skipping destruction.
template <typename T>
class R__CLING_PTRCHECK(off) SmallVectorTemplateBase<T, true> : public SmallVectorTemplateCommon<T> {
   using SuperClass = SmallVectorTemplateCommon<T>;

protected:
   SmallVectorTemplateBase(size_t Size) : SmallVectorTemplateCommon<T>(Size) {}

   // No need to do a destroy loop for POD's.
   static void destroy_range(T *, T *) {}

   /// Move the range [I, E) onto the uninitialized memory
   /// starting with "Dest", constructing elements into it as needed.
   template <typename It1, typename It2>
   static void uninitialized_move(It1 I, It1 E, It2 Dest)
   {
      // Just do a copy.
      uninitialized_copy(I, E, Dest);
   }

   /// Copy the range [I, E) onto the uninitialized memory
   /// starting with "Dest", constructing elements into it as needed.
   template <typename It1, typename It2>
   static void uninitialized_copy(It1 I, It1 E, It2 Dest)
   {
      // Arbitrary iterator types; just use the basic implementation.
      std::uninitialized_copy(I, E, Dest);
   }

   /// Copy the range [I, E) onto the uninitialized memory
   /// starting with "Dest", constructing elements into it as needed.
   template <typename T1, typename T2>
   static void uninitialized_copy(
      T1 *I, T1 *E, T2 *Dest,
      typename std::enable_if<std::is_same<typename std::remove_const<T1>::type, T2>::value>::type * = nullptr)
   {
      // Use memcpy for PODs iterated by pointers (which includes SmallVector
      // iterators): std::uninitialized_copy optimizes to memmove, but we can
      // use memcpy here. Note that I and E are iterators and thus might be
      // invalid for memcpy if they are equal.
      if (I != E)
         memcpy(reinterpret_cast<void *>(Dest), I, (E - I) * sizeof(T));
   }

   /// Double the size of the allocated memory, guaranteeing space for at
   /// least one more element or MinSize if specified.
   void grow(size_t MinSize = 0)
   {
      this->grow_pod(MinSize, sizeof(T));
   }

public:
   using iterator = typename SuperClass::iterator;
   using const_iterator = typename SuperClass::const_iterator;
   using reference = typename SuperClass::reference;
   using size_type = typename SuperClass::size_type;

   void push_back(const T &Elt)
   {
      if (R__unlikely(this->size() >= this->capacity()))
         this->grow();
      memcpy(reinterpret_cast<void *>(this->end()), &Elt, sizeof(T));
      this->set_size(this->size() + 1);
   }

   void pop_back() { this->set_size(this->size() - 1); }
};

/// Storage for the SmallVector elements.  This is specialized for the N=0 case
/// to avoid allocating unnecessary storage.
template <typename T, unsigned N>
struct SmallVectorStorage {
   alignas(T) char InlineElts[N * sizeof(T)]{};
};

/// We need the storage to be properly aligned even for small-size of 0 so that
/// the pointer math in \a SmallVectorTemplateCommon::getFirstEl() is
/// well-defined.
template <typename T>
struct alignas(T) SmallVectorStorage<T, 0> {
};

/// The size of the inline storage of an RVec.
/// Our policy is to allocate at least 8 elements (or more if they all fit into one cacheline)
/// unless the size of the buffer with 8 elements would be over a certain maximum size.
template <typename T>
struct RVecInlineStorageSize {
private:
#ifdef R__HAS_HARDWARE_INTERFERENCE_SIZE
   static constexpr std::size_t cacheLineSize = std::hardware_destructive_interference_size;
#else
   // safe bet: assume the typical 64 bytes
   static constexpr std::size_t cacheLineSize = 64;
#endif
   static constexpr unsigned elementsPerCacheLine = (cacheLineSize - sizeof(SmallVectorBase)) / sizeof(T);
   static constexpr unsigned maxInlineByteSize = 1024;

public:
   static constexpr unsigned value =
      elementsPerCacheLine >= 8 ? elementsPerCacheLine : (sizeof(T) * 8 > maxInlineByteSize ? 0 : 8);
};

// A C++14-compatible implementation of std::uninitialized_value_construct
template <typename ForwardIt>
void UninitializedValueConstruct(ForwardIt first, ForwardIt last)
{
#if __cplusplus < 201703L
   for (; first != last; ++first)
      new (static_cast<void *>(std::addressof(*first))) typename std::iterator_traits<ForwardIt>::value_type();
#else
   std::uninitialized_value_construct(first, last);
#endif
}

/// An unsafe function to reset the buffer for which this RVec is acting as a view.
///
/// \note This is a low-level method that _must_ be called on RVecs that are already non-owning:
/// - it does not put the RVec in "non-owning mode" (fCapacity == -1)
/// - it does not free any owned buffer
template <typename T>
void ResetView(RVec<T> &v, T* addr, std::size_t sz)
{
   v.fBeginX = addr;
   v.fSize = sz;
}

} // namespace VecOps
} // namespace Internal

namespace Detail {
namespace VecOps {

/// This class consists of common code factored out of the SmallVector class to
/// reduce code duplication based on the SmallVector 'N' template parameter.
template <typename T>
class R__CLING_PTRCHECK(off) RVecImpl : public Internal::VecOps::SmallVectorTemplateBase<T> {
   using SuperClass = Internal::VecOps::SmallVectorTemplateBase<T>;

public:
   using iterator = typename SuperClass::iterator;
   using const_iterator = typename SuperClass::const_iterator;
   using reference = typename SuperClass::reference;
   using size_type = typename SuperClass::size_type;

protected:
   // Default ctor - Initialize to empty.
   explicit RVecImpl(unsigned N) : ROOT::Internal::VecOps::SmallVectorTemplateBase<T>(N) {}

public:
   RVecImpl(const RVecImpl &) = delete;

   ~RVecImpl()
   {
      // Subclass has already destructed this vector's elements.
      // If this wasn't grown from the inline copy, deallocate the old space.
      if (!this->isSmall() && this->Owns())
         free(this->begin());
   }

   // also give up adopted memory if applicable
   void clear()
   {
      if (this->Owns()) {
         this->destroy_range(this->begin(), this->end());
         this->fSize = 0;
      } else {
         this->resetToSmall();
      }
   }

   void resize(size_type N)
   {
      if (N < this->size()) {
         if (this->Owns())
            this->destroy_range(this->begin() + N, this->end());
         this->set_size(N);
      } else if (N > this->size()) {
         if (this->capacity() < N)
            this->grow(N);
         for (auto I = this->end(), E = this->begin() + N; I != E; ++I)
            new (&*I) T();
         this->set_size(N);
      }
   }

   void resize(size_type N, const T &NV)
   {
      if (N < this->size()) {
         if (this->Owns())
            this->destroy_range(this->begin() + N, this->end());
         this->set_size(N);
      } else if (N > this->size()) {
         if (this->capacity() < N)
            this->grow(N);
         std::uninitialized_fill(this->end(), this->begin() + N, NV);
         this->set_size(N);
      }
   }

   void reserve(size_type N)
   {
      if (this->capacity() < N)
         this->grow(N);
   }

   void pop_back_n(size_type NumItems)
   {
      if (this->size() < NumItems) {
         throw std::runtime_error("Popping back more elements than those available.");
      }
      if (this->Owns())
         this->destroy_range(this->end() - NumItems, this->end());
      this->set_size(this->size() - NumItems);
   }

   R__RVEC_NODISCARD T pop_back_val()
   {
      T Result = ::std::move(this->back());
      this->pop_back();
      return Result;
   }

   void swap(RVecImpl &RHS);

   /// Add the specified range to the end of the SmallVector.
   template <typename in_iter,
             typename = typename std::enable_if<std::is_convertible<
                typename std::iterator_traits<in_iter>::iterator_category, std::input_iterator_tag>::value>::type>
   void append(in_iter in_start, in_iter in_end)
   {
      size_type NumInputs = std::distance(in_start, in_end);
      if (NumInputs > this->capacity() - this->size())
         this->grow(this->size() + NumInputs);

      this->uninitialized_copy(in_start, in_end, this->end());
      this->set_size(this->size() + NumInputs);
   }

   /// Append \p NumInputs copies of \p Elt to the end.
   void append(size_type NumInputs, const T &Elt)
   {
      if (NumInputs > this->capacity() - this->size())
         this->grow(this->size() + NumInputs);

      std::uninitialized_fill_n(this->end(), NumInputs, Elt);
      this->set_size(this->size() + NumInputs);
   }

   void append(std::initializer_list<T> IL) { append(IL.begin(), IL.end()); }

   // from the original LLVM implementation:
   // FIXME: Consider assigning over existing elements, rather than clearing &
   // re-initializing them - for all assign(...) variants.

   void assign(size_type NumElts, const T &Elt)
   {
      clear();
      if (this->capacity() < NumElts)
         this->grow(NumElts);
      this->set_size(NumElts);
      std::uninitialized_fill(this->begin(), this->end(), Elt);
   }

   template <typename in_iter,
             typename = typename std::enable_if<std::is_convertible<
                typename std::iterator_traits<in_iter>::iterator_category, std::input_iterator_tag>::value>::type>
   void assign(in_iter in_start, in_iter in_end)
   {
      clear();
      append(in_start, in_end);
   }

   void assign(std::initializer_list<T> IL)
   {
      clear();
      append(IL);
   }

   iterator erase(const_iterator CI)
   {
      // Just cast away constness because this is a non-const member function.
      iterator I = const_cast<iterator>(CI);

      if (I < this->begin() || I >= this->end()) {
         throw std::runtime_error("The iterator passed to `erase` is out of bounds.");
      }

      iterator N = I;
      // Shift all elts down one.
      std::move(I + 1, this->end(), I);
      // Drop the last elt.
      this->pop_back();
      return (N);
   }

   iterator erase(const_iterator CS, const_iterator CE)
   {
      // Just cast away constness because this is a non-const member function.
      iterator S = const_cast<iterator>(CS);
      iterator E = const_cast<iterator>(CE);

      if (S < this->begin() || E > this->end() || S > E) {
         throw std::runtime_error("Invalid start/end pair passed to `erase` (out of bounds or start > end).");
      }

      iterator N = S;
      // Shift all elts down.
      iterator I = std::move(E, this->end(), S);
      // Drop the last elts.
      if (this->Owns())
         this->destroy_range(I, this->end());
      this->set_size(I - this->begin());
      return (N);
   }

   iterator insert(iterator I, T &&Elt)
   {
      if (I == this->end()) { // Important special case for empty vector.
         this->push_back(::std::move(Elt));
         return this->end() - 1;
      }

      if (I < this->begin() || I > this->end()) {
         throw std::runtime_error("The iterator passed to `insert` is out of bounds.");
      }

      if (this->size() >= this->capacity()) {
         size_t EltNo = I - this->begin();
         this->grow();
         I = this->begin() + EltNo;
      }

      ::new ((void *)this->end()) T(::std::move(this->back()));
      // Push everything else over.
      std::move_backward(I, this->end() - 1, this->end());
      this->set_size(this->size() + 1);

      // If we just moved the element we're inserting, be sure to update
      // the reference.
      T *EltPtr = &Elt;
      if (I <= EltPtr && EltPtr < this->end())
         ++EltPtr;

      *I = ::std::move(*EltPtr);
      return I;
   }

   iterator insert(iterator I, const T &Elt)
   {
      if (I == this->end()) { // Important special case for empty vector.
         this->push_back(Elt);
         return this->end() - 1;
      }

      if (I < this->begin() || I > this->end()) {
         throw std::runtime_error("The iterator passed to `insert` is out of bounds.");
      }

      if (this->size() >= this->capacity()) {
         size_t EltNo = I - this->begin();
         this->grow();
         I = this->begin() + EltNo;
      }
      ::new ((void *)this->end()) T(std::move(this->back()));
      // Push everything else over.
      std::move_backward(I, this->end() - 1, this->end());
      this->set_size(this->size() + 1);

      // If we just moved the element we're inserting, be sure to update
      // the reference.
      const T *EltPtr = &Elt;
      if (I <= EltPtr && EltPtr < this->end())
         ++EltPtr;

      *I = *EltPtr;
      return I;
   }

   iterator insert(iterator I, size_type NumToInsert, const T &Elt)
   {
      // Convert iterator to elt# to avoid invalidating iterator when we reserve()
      size_t InsertElt = I - this->begin();

      if (I == this->end()) { // Important special case for empty vector.
         append(NumToInsert, Elt);
         return this->begin() + InsertElt;
      }

      if (I < this->begin() || I > this->end()) {
         throw std::runtime_error("The iterator passed to `insert` is out of bounds.");
      }

      // Ensure there is enough space.
      reserve(this->size() + NumToInsert);

      // Uninvalidate the iterator.
      I = this->begin() + InsertElt;

      // If there are more elements between the insertion point and the end of the
      // range than there are being inserted, we can use a simple approach to
      // insertion.  Since we already reserved space, we know that this won't
      // reallocate the vector.
      if (size_t(this->end() - I) >= NumToInsert) {
         T *OldEnd = this->end();
         append(std::move_iterator<iterator>(this->end() - NumToInsert), std::move_iterator<iterator>(this->end()));

         // Copy the existing elements that get replaced.
         std::move_backward(I, OldEnd - NumToInsert, OldEnd);

         std::fill_n(I, NumToInsert, Elt);
         return I;
      }

      // Otherwise, we're inserting more elements than exist already, and we're
      // not inserting at the end.

      // Move over the elements that we're about to overwrite.
      T *OldEnd = this->end();
      this->set_size(this->size() + NumToInsert);
      size_t NumOverwritten = OldEnd - I;
      this->uninitialized_move(I, OldEnd, this->end() - NumOverwritten);

      // Replace the overwritten part.
      std::fill_n(I, NumOverwritten, Elt);

      // Insert the non-overwritten middle part.
      std::uninitialized_fill_n(OldEnd, NumToInsert - NumOverwritten, Elt);
      return I;
   }

   template <typename ItTy,
             typename = typename std::enable_if<std::is_convertible<
                typename std::iterator_traits<ItTy>::iterator_category, std::input_iterator_tag>::value>::type>
   iterator insert(iterator I, ItTy From, ItTy To)
   {
      // Convert iterator to elt# to avoid invalidating iterator when we reserve()
      size_t InsertElt = I - this->begin();

      if (I == this->end()) { // Important special case for empty vector.
         append(From, To);
         return this->begin() + InsertElt;
      }

      if (I < this->begin() || I > this->end()) {
         throw std::runtime_error("The iterator passed to `insert` is out of bounds.");
      }

      size_t NumToInsert = std::distance(From, To);

      // Ensure there is enough space.
      reserve(this->size() + NumToInsert);

      // Uninvalidate the iterator.
      I = this->begin() + InsertElt;

      // If there are more elements between the insertion point and the end of the
      // range than there are being inserted, we can use a simple approach to
      // insertion.  Since we already reserved space, we know that this won't
      // reallocate the vector.
      if (size_t(this->end() - I) >= NumToInsert) {
         T *OldEnd = this->end();
         append(std::move_iterator<iterator>(this->end() - NumToInsert), std::move_iterator<iterator>(this->end()));

         // Copy the existing elements that get replaced.
         std::move_backward(I, OldEnd - NumToInsert, OldEnd);

         std::copy(From, To, I);
         return I;
      }

      // Otherwise, we're inserting more elements than exist already, and we're
      // not inserting at the end.

      // Move over the elements that we're about to overwrite.
      T *OldEnd = this->end();
      this->set_size(this->size() + NumToInsert);
      size_t NumOverwritten = OldEnd - I;
      this->uninitialized_move(I, OldEnd, this->end() - NumOverwritten);

      // Replace the overwritten part.
      for (T *J = I; NumOverwritten > 0; --NumOverwritten) {
         *J = *From;
         ++J;
         ++From;
      }

      // Insert the non-overwritten middle part.
      this->uninitialized_copy(From, To, OldEnd);
      return I;
   }

   void insert(iterator I, std::initializer_list<T> IL) { insert(I, IL.begin(), IL.end()); }

   template <typename... ArgTypes>
   reference emplace_back(ArgTypes &&...Args)
   {
      if (R__unlikely(this->size() >= this->capacity()))
         this->grow();
      ::new ((void *)this->end()) T(std::forward<ArgTypes>(Args)...);
      this->set_size(this->size() + 1);
      return this->back();
   }

   RVecImpl &operator=(const RVecImpl &RHS);

   RVecImpl &operator=(RVecImpl &&RHS);
};

template <typename T>
void RVecImpl<T>::swap(RVecImpl<T> &RHS)
{
   if (this == &RHS)
      return;

   // We can only avoid copying elements if neither vector is small.
   if (!this->isSmall() && !RHS.isSmall()) {
      std::swap(this->fBeginX, RHS.fBeginX);
      std::swap(this->fSize, RHS.fSize);
      std::swap(this->fCapacity, RHS.fCapacity);
      return;
   }

   // This block handles the swap of a small and a non-owning vector
   // It is more efficient to first move the non-owning vector, hence the 2 cases
   if (this->isSmall() && !RHS.Owns()) { // the right vector is non-owning
      RVecImpl<T> temp(0);
      temp = std::move(RHS);
      RHS = std::move(*this);
      *this = std::move(temp);
      return;
   } else if (RHS.isSmall() && !this->Owns()) { // the left vector is non-owning
      RVecImpl<T> temp(0);
      temp = std::move(*this);
      *this = std::move(RHS);
      RHS = std::move(temp);
      return;
   }

   if (RHS.size() > this->capacity())
      this->grow(RHS.size());
   if (this->size() > RHS.capacity())
      RHS.grow(this->size());

   // Swap the shared elements.
   size_t NumShared = this->size();
   if (NumShared > RHS.size())
      NumShared = RHS.size();
   for (size_type i = 0; i != NumShared; ++i)
      std::iter_swap(this->begin() + i, RHS.begin() + i);

   // Copy over the extra elts.
   if (this->size() > RHS.size()) {
      size_t EltDiff = this->size() - RHS.size();
      this->uninitialized_copy(this->begin() + NumShared, this->end(), RHS.end());
      RHS.set_size(RHS.size() + EltDiff);
      if (this->Owns())
         this->destroy_range(this->begin() + NumShared, this->end());
      this->set_size(NumShared);
   } else if (RHS.size() > this->size()) {
      size_t EltDiff = RHS.size() - this->size();
      this->uninitialized_copy(RHS.begin() + NumShared, RHS.end(), this->end());
      this->set_size(this->size() + EltDiff);
      if (RHS.Owns())
         this->destroy_range(RHS.begin() + NumShared, RHS.end());
      RHS.set_size(NumShared);
   }
}

template <typename T>
RVecImpl<T> &RVecImpl<T>::operator=(const RVecImpl<T> &RHS)
{
   // Avoid self-assignment.
   if (this == &RHS)
      return *this;

   // If we already have sufficient space, assign the common elements, then
   // destroy any excess.
   size_t RHSSize = RHS.size();
   size_t CurSize = this->size();
   if (CurSize >= RHSSize) {
      // Assign common elements.
      iterator NewEnd;
      if (RHSSize)
         NewEnd = std::copy(RHS.begin(), RHS.begin() + RHSSize, this->begin());
      else
         NewEnd = this->begin();

      // Destroy excess elements.
      if (this->Owns())
         this->destroy_range(NewEnd, this->end());

      // Trim.
      this->set_size(RHSSize);
      return *this;
   }

   // If we have to grow to have enough elements, destroy the current elements.
   // This allows us to avoid copying them during the grow.
   // From the original LLVM implementation:
   // FIXME: don't do this if they're efficiently moveable.
   if (this->capacity() < RHSSize) {
      if (this->Owns()) {
         // Destroy current elements.
         this->destroy_range(this->begin(), this->end());
      }
      this->set_size(0);
      CurSize = 0;
      this->grow(RHSSize);
   } else if (CurSize) {
      // Otherwise, use assignment for the already-constructed elements.
      std::copy(RHS.begin(), RHS.begin() + CurSize, this->begin());
   }

   // Copy construct the new elements in place.
   this->uninitialized_copy(RHS.begin() + CurSize, RHS.end(), this->begin() + CurSize);

   // Set end.
   this->set_size(RHSSize);
   return *this;
}

template <typename T>
RVecImpl<T> &RVecImpl<T>::operator=(RVecImpl<T> &&RHS)
{
   // Avoid self-assignment.
   if (this == &RHS)
      return *this;

   // If the RHS isn't small, clear this vector and then steal its buffer.
   if (!RHS.isSmall()) {
      if (this->Owns()) {
         this->destroy_range(this->begin(), this->end());
         if (!this->isSmall())
            free(this->begin());
      }
      this->fBeginX = RHS.fBeginX;
      this->fSize = RHS.fSize;
      this->fCapacity = RHS.fCapacity;
      RHS.resetToSmall();
      return *this;
   }

   // If we already have sufficient space, assign the common elements, then
   // destroy any excess.
   size_t RHSSize = RHS.size();
   size_t CurSize = this->size();
   if (CurSize >= RHSSize) {
      // Assign common elements.
      iterator NewEnd = this->begin();
      if (RHSSize)
         NewEnd = std::move(RHS.begin(), RHS.end(), NewEnd);

      // Destroy excess elements and trim the bounds.
      if (this->Owns())
         this->destroy_range(NewEnd, this->end());
      this->set_size(RHSSize);

      // Clear the RHS.
      RHS.clear();

      return *this;
   }

   // If we have to grow to have enough elements, destroy the current elements.
   // This allows us to avoid copying them during the grow.
   // From the original LLVM implementation:
   // FIXME: this may not actually make any sense if we can efficiently move
   // elements.
   if (this->capacity() < RHSSize) {
      if (this->Owns()) {
         // Destroy current elements.
         this->destroy_range(this->begin(), this->end());
      }
      this->set_size(0);
      CurSize = 0;
      this->grow(RHSSize);
   } else if (CurSize) {
      // Otherwise, use assignment for the already-constructed elements.
      std::move(RHS.begin(), RHS.begin() + CurSize, this->begin());
   }

   // Move-construct the new elements in place.
   this->uninitialized_move(RHS.begin() + CurSize, RHS.end(), this->begin() + CurSize);

   // Set end.
   this->set_size(RHSSize);

   RHS.clear();
   return *this;
}

template <typename T>
bool IsSmall(const ROOT::VecOps::RVec<T> &v)
{
   return v.isSmall();
}

template <typename T>
bool IsAdopting(const ROOT::VecOps::RVec<T> &v)
{
   return !v.Owns();
}

} // namespace VecOps
} // namespace Detail

namespace VecOps {
// Note that we open here with @{ the Doxygen group vecops and it is
// closed again at the end of the C++ namespace VecOps
/**
  * \defgroup vecops VecOps
  * A "std::vector"-like collection of values implementing handy operation to analyse them
  * @{
*/

// From the original SmallVector code:
// This is a 'vector' (really, a variable-sized array), optimized
// for the case when the array is small.  It contains some number of elements
// in-place, which allows it to avoid heap allocation when the actual number of
// elements is below that threshold.  This allows normal "small" cases to be
// fast without losing generality for large inputs.
//
// Note that this does not attempt to be exception safe.

template <typename T, unsigned int N>
class R__CLING_PTRCHECK(off) RVecN : public Detail::VecOps::RVecImpl<T>, Internal::VecOps::SmallVectorStorage<T, N> {
public:
   RVecN() : Detail::VecOps::RVecImpl<T>(N) {}

   ~RVecN()
   {
      if (this->Owns()) {
         // Destroy the constructed elements in the vector.
         this->destroy_range(this->begin(), this->end());
      }
   }

   explicit RVecN(size_t Size, const T &Value) : Detail::VecOps::RVecImpl<T>(N) { this->assign(Size, Value); }

   explicit RVecN(size_t Size) : Detail::VecOps::RVecImpl<T>(N)
   {
      if (Size > N)
         this->grow(Size);
      this->fSize = Size;
      ROOT::Internal::VecOps::UninitializedValueConstruct(this->begin(), this->end());
   }

   template <typename ItTy,
             typename = typename std::enable_if<std::is_convertible<
                typename std::iterator_traits<ItTy>::iterator_category, std::input_iterator_tag>::value>::type>
   RVecN(ItTy S, ItTy E) : Detail::VecOps::RVecImpl<T>(N)
   {
      this->append(S, E);
   }

   RVecN(std::initializer_list<T> IL) : Detail::VecOps::RVecImpl<T>(N) { this->assign(IL); }

   RVecN(const RVecN &RHS) : Detail::VecOps::RVecImpl<T>(N)
   {
      if (!RHS.empty())
         Detail::VecOps::RVecImpl<T>::operator=(RHS);
   }

   RVecN &operator=(const RVecN &RHS)
   {
      Detail::VecOps::RVecImpl<T>::operator=(RHS);
      return *this;
   }

   RVecN(RVecN &&RHS) : Detail::VecOps::RVecImpl<T>(N)
   {
      if (!RHS.empty())
         Detail::VecOps::RVecImpl<T>::operator=(::std::move(RHS));
   }

   RVecN(Detail::VecOps::RVecImpl<T> &&RHS) : Detail::VecOps::RVecImpl<T>(N)
   {
      if (!RHS.empty())
         Detail::VecOps::RVecImpl<T>::operator=(::std::move(RHS));
   }

   RVecN(const std::vector<T> &RHS) : RVecN(RHS.begin(), RHS.end()) {}

   RVecN &operator=(RVecN &&RHS)
   {
      Detail::VecOps::RVecImpl<T>::operator=(::std::move(RHS));
      return *this;
   }

   RVecN(T* p, size_t n) : Detail::VecOps::RVecImpl<T>(N)
   {
      this->fBeginX = p;
      this->fSize = n;
      this->fCapacity = -1;
   }

   RVecN &operator=(Detail::VecOps::RVecImpl<T> &&RHS)
   {
      Detail::VecOps::RVecImpl<T>::operator=(::std::move(RHS));
      return *this;
   }

   RVecN &operator=(std::initializer_list<T> IL)
   {
      this->assign(IL);
      return *this;
   }

   using reference = typename Internal::VecOps::SmallVectorTemplateCommon<T>::reference;
   using const_reference = typename Internal::VecOps::SmallVectorTemplateCommon<T>::const_reference;
   using size_type = typename Internal::VecOps::SmallVectorTemplateCommon<T>::size_type;
   using value_type = typename Internal::VecOps::SmallVectorTemplateCommon<T>::value_type;
   using Internal::VecOps::SmallVectorTemplateCommon<T>::begin;
   using Internal::VecOps::SmallVectorTemplateCommon<T>::size;

   reference operator[](size_type idx)
   {
      return begin()[idx];
   }

   const_reference operator[](size_type idx) const
   {
      return begin()[idx];
   }

   template <typename V, unsigned M, typename = std::enable_if<std::is_convertible<V, bool>::value>>
   RVecN operator[](const RVecN<V, M> &conds) const
   {
      const size_type n = conds.size();

      if (n != this->size()) {
         std::string msg = "Cannot index RVecN of size " + std::to_string(this->size()) +
                           " with condition vector of different size (" + std::to_string(n) + ").";
         throw std::runtime_error(msg);
      }

      size_type n_true = 0ull;
      for (auto c : conds)
         n_true += c; // relies on bool -> int conversion, faster than branching

      RVecN ret;
      ret.reserve(n_true);
      size_type j = 0u;
      for (size_type i = 0u; i < n; ++i) {
         if (conds[i]) {
            ret.push_back(this->operator[](i));
            ++j;
         }
      }
      return ret;
   }

   // conversion
   template <typename U, unsigned M, typename = std::enable_if<std::is_convertible<T, U>::value>>
   operator RVecN<U, M>() const
   {
      return RVecN<U, M>(this->begin(), this->end());
   }

   reference at(size_type pos)
   {
      if (pos >= size_type(this->fSize)) {
         std::string msg = "RVecN::at: size is " + std::to_string(this->fSize) + " but out-of-bounds index " +
                           std::to_string(pos) + " was requested.";
         throw std::out_of_range(msg);
      }
      return this->operator[](pos);
   }

   const_reference at(size_type pos) const
   {
      if (pos >= size_type(this->fSize)) {
         std::string msg = "RVecN::at: size is " + std::to_string(this->fSize) + " but out-of-bounds index " +
                           std::to_string(pos) + " was requested.";
         throw std::out_of_range(msg);
      }
      return this->operator[](pos);
   }

   /// No exception thrown. The user specifies the desired value in case the RVecN is shorter than `pos`.
   value_type at(size_type pos, value_type fallback)
   {
      if (pos >= size_type(this->fSize))
         return fallback;
      return this->operator[](pos);
   }

   /// No exception thrown. The user specifies the desired value in case the RVecN is shorter than `pos`.
   value_type at(size_type pos, value_type fallback) const
   {
      if (pos >= size_type(this->fSize))
         return fallback;
      return this->operator[](pos);
   }
};

// clang-format off
/**
\class ROOT::VecOps::RVec
\brief A "std::vector"-like collection of values implementing handy operation to analyse them
\tparam T The type of the contained objects

A RVec is a container designed to make analysis of values' collections fast and easy.
Its storage is contiguous in memory and its interface is designed such to resemble to the one
of the stl vector. In addition the interface features methods and
[external functions](https://root.cern/doc/master/namespaceROOT_1_1VecOps.html) to ease the manipulation and analysis
of the data in the RVec.

\note ROOT::VecOps::RVec can also be spelled simply ROOT::RVec. Shorthand aliases such as ROOT::RVecI or ROOT::RVecD
are also available as template instantiations of RVec of fundamental types. The full list of available aliases:
- RVecB (`bool`)
- RVecC (`char`)
- RVecD (`double`)
- RVecF (`float`)
- RVecI (`int`)
- RVecL (`long`)
- RVecLL (`long long`)
- RVecU (`unsigned`)
- RVecUL (`unsigned long`)
- RVecULL (`unsigned long long`)

\note RVec does not attempt to be exception safe. Exceptions thrown by element constructors during insertions, swaps or
other operations will be propagated potentially leaving the RVec object in an invalid state.

\note RVec methods (e.g. `at` or `size`) follow the STL naming convention instead of the ROOT naming convention in order
to make RVec a drop-in replacement for `std::vector`.

\htmlonly
<a href="https://doi.org/10.5281/zenodo.1253756"><img src="https://zenodo.org/badge/DOI/10.5281/zenodo.1253756.svg" alt="DOI"></a>
\endhtmlonly

## Table of Contents
- [Example](\ref example)
- [Arithmetic operations, logical operations and mathematical functions](\ref operationsandfunctions)
- [Owning and adopting memory](\ref owningandadoptingmemory)
- [Sorting and manipulation of indices](\ref sorting)
- [Usage in combination with RDataFrame](\ref usagetdataframe)
- [Reference for the RVec class](\ref RVecdoxyref)
- [Reference for RVec helper functions](https://root.cern/doc/master/namespaceROOT_1_1VecOps.html)

\anchor example
## Example
Suppose to have an event featuring a collection of muons with a certain pseudorapidity,
momentum and charge, e.g.:
~~~{.cpp}
std::vector<short> mu_charge {1, 1, -1, -1, -1, 1, 1, -1};
std::vector<float> mu_pt {56, 45, 32, 24, 12, 8, 7, 6.2};
std::vector<float> mu_eta {3.1, -.2, -1.1, 1, 4.1, 1.6, 2.4, -.5};
~~~
Suppose you want to extract the transverse momenta of the muons satisfying certain
criteria, for example consider only negatively charged muons with a pseudorapidity
smaller or equal to 2 and with a transverse momentum greater than 10 GeV.
Such a selection would require, among the other things, the management of an explicit
loop, for example:
~~~{.cpp}
std::vector<float> goodMuons_pt;
const auto size = mu_charge.size();
for (size_t i=0; i < size; ++i) {
   if (mu_pt[i] > 10 && abs(mu_eta[i]) <= 2. &&  mu_charge[i] == -1) {
      goodMuons_pt.emplace_back(mu_pt[i]);
   }
}
~~~
These operations become straightforward with RVec - we just need to *write what
we mean*:
~~~{.cpp}
auto goodMuons_pt = mu_pt[ (mu_pt > 10.f && abs(mu_eta) <= 2.f && mu_charge == -1) ]
~~~
Now the clean collection of transverse momenta can be used within the rest of the data analysis, for
example to fill a histogram.

\anchor operationsandfunctions
## Arithmetic operations, logical operations and mathematical functions
Arithmetic operations on RVec instances can be performed: for example, they can be added, subtracted, multiplied.
~~~{.cpp}
RVec<double> v1 {1.,2.,3.,4.};
RVec<float> v2 {5.f,6.f,7.f,8.f};
auto v3 = v1+v2;
auto v4 = 3 * v1;
~~~
The supported operators are 
 - +, -, *, /
 - +=, -=, *=, /=
 - <, >, ==, !=, <=, >=, &&, ||
 - ~, !
 - &, |, ^
 - &=, |=, ^=
 - <<=, >>=

The most common mathematical functions are supported. It is possible to invoke them passing 
RVecs as arguments.
 - abs, fdim, fmod, remainder
 - floor, ceil, trunc, round, lround, llround
 - exp, exp2, expm1
 - log, log10, log2, log1p
 - pow
 - sqrt, cbrt
 - sin, cos, tan, asin, acos, atan, atan2, hypot
 - sinh, cosh, tanh, asinh, acosh
 - erf, erfc
 - lgamma, tgamma

If the VDT library is available, the following functions can be invoked. Internally the calculations
are vectorized:
 - fast_expf, fast_logf, fast_sinf, fast_cosf, fast_tanf, fast_asinf, fast_acosf, fast_atanf
 - fast_exp, fast_log, fast_sin, fast_cos, fast_tan, fast_asin, fast_acos, fast_atan

\anchor owningandadoptingmemory
## Owning and adopting memory
RVec has contiguous memory associated to it. It can own it or simply adopt it. In the latter case,
it can be constructed with the address of the memory associated to it and its length. For example:
~~~{.cpp}
std::vector<int> myStlVec {1,2,3};
RVec<int> myRVec(myStlVec.data(), myStlVec.size());
~~~
In this case, the memory associated to myStlVec and myRVec is the same, myRVec simply "adopted it".
If any method which implies a re-allocation is called, e.g. *emplace_back* or *resize*, the adopted
memory is released and new one is allocated. The previous content is copied in the new memory and
preserved.

\anchor sorting
## Sorting and manipulation of indices

### Sorting
RVec complies to the STL interfaces when it comes to iterations. As a result, standard algorithms
can be used, for example sorting:
~~~{.cpp}
RVec<double> v{6., 4., 5.};
std::sort(v.begin(), v.end());
~~~

For convenience, helpers are provided too:
~~~{.cpp}
auto sorted_v = Sort(v);
auto reversed_v = Reverse(v);
~~~

### Manipulation of indices

It is also possible to manipulated the RVecs acting on their indices. For example,
the following syntax
~~~{.cpp}
RVecD v0 {9., 7., 8.};
auto v1 = Take(v0, {1, 2, 0});
~~~
will yield a new RVec<double> the content of which is the first, second and zeroth element of
v0, i.e. `{7., 8., 9.}`.

The `Argsort` and `StableArgsort` helper extracts the indices which order the content of a `RVec`.
For example, this snippet accomplishes in a more expressive way what we just achieved:
~~~{.cpp}
auto v1_indices = Argsort(v0); // The content of v1_indices is {1, 2, 0}.
v1 = Take(v0, v1_indices);
~~~

The `Take` utility allows to extract portions of the `RVec`. The content to be *taken*
can be specified with an `RVec` of indices or an integer. If the integer is negative,
elements will be picked starting from the end of the container:
~~~{.cpp}
RVecF vf {1.f, 2.f, 3.f, 4.f};
auto vf_1 = Take(vf, {1, 3}); // The content is {2.f, 4.f}
auto vf_2 = Take(vf, 2); // The content is {1.f, 2.f}
auto vf_3 = Take(vf, -3); // The content is {2.f, 3.f, 4.f}
~~~

\anchor usagetdataframe
## Usage in combination with RDataFrame
RDataFrame leverages internally RVecs. Suppose to have a dataset stored in a
TTree which holds these columns (here we choose C arrays to represent the
collections, they could be as well std::vector instances):
~~~{.bash}
  nPart            "nPart/I"            An integer representing the number of particles
  px               "px[nPart]/D"        The C array of the particles' x component of the momentum
  py               "py[nPart]/D"        The C array of the particles' y component of the momentum
  E                "E[nPart]/D"         The C array of the particles' Energy
~~~
Suppose you'd like to plot in a histogram the transverse momenta of all particles
for which the energy is greater than 200 MeV.
The code required would just be:
~~~{.cpp}
RDataFrame d("mytree", "myfile.root");
auto cutPt = [](RVecD &pxs, RVecD &pys, RVecD &Es) {
   auto all_pts = sqrt(pxs * pxs + pys * pys);
   auto good_pts = all_pts[Es > 200.];
   return good_pts;
   };

auto hpt = d.Define("pt", cutPt, {"px", "py", "E"})
            .Histo1D("pt");
hpt->Draw();
~~~
And if you'd like to express your selection as a string:
~~~{.cpp}
RDataFrame d("mytree", "myfile.root");
auto hpt = d.Define("pt", "sqrt(pxs * pxs + pys * pys)[E>200]")
            .Histo1D("pt");
hpt->Draw();
~~~
\anchor RVecdoxyref
**/
// clang-format on

template <typename T>
class R__CLING_PTRCHECK(off) RVec : public RVecN<T, Internal::VecOps::RVecInlineStorageSize<T>::value> {
   using SuperClass = RVecN<T, Internal::VecOps::RVecInlineStorageSize<T>::value>;

   friend void Internal::VecOps::ResetView<>(RVec<T> &v, T *addr, std::size_t sz);

public:
   using reference = typename SuperClass::reference;
   using const_reference = typename SuperClass::const_reference;
   using size_type = typename SuperClass::size_type;
   using value_type = typename SuperClass::value_type;
   using SuperClass::begin;
   using SuperClass::size;

   RVec() {}

   explicit RVec(size_t Size, const T &Value) : SuperClass(Size, Value) {}

   explicit RVec(size_t Size) : SuperClass(Size) {}

   template <typename ItTy,
             typename = typename std::enable_if<std::is_convertible<
                typename std::iterator_traits<ItTy>::iterator_category, std::input_iterator_tag>::value>::type>
   RVec(ItTy S, ItTy E) : SuperClass(S, E)
   {
   }

   RVec(std::initializer_list<T> IL) : SuperClass(IL) {}

   RVec(const RVec &RHS) : SuperClass(RHS) {}

   RVec &operator=(const RVec &RHS)
   {
      SuperClass::operator=(RHS);
      return *this;
   }

   RVec(RVec &&RHS) : SuperClass(std::move(RHS)) {}

   RVec &operator=(RVec &&RHS)
   {
      SuperClass::operator=(std::move(RHS));
      return *this;
   }

   RVec(Detail::VecOps::RVecImpl<T> &&RHS) : SuperClass(std::move(RHS)) {}

   template <unsigned N>
   RVec(RVecN<T, N> &&RHS) : SuperClass(std::move(RHS)) {}

   template <unsigned N>
   RVec(const RVecN<T, N> &RHS) : SuperClass(RHS) {}

   RVec(const std::vector<T> &RHS) : SuperClass(RHS) {}

   RVec(T* p, size_t n) : SuperClass(p, n) {}

   // conversion
   template <typename U, typename = std::enable_if<std::is_convertible<T, U>::value>>
   operator RVec<U>() const
   {
      return RVec<U>(this->begin(), this->end());
   }

   using SuperClass::operator[];

   template <typename V, typename = std::enable_if<std::is_convertible<V, bool>::value>>
   RVec operator[](const RVec<V> &conds) const
   {
      return RVec(SuperClass::operator[](conds));
   }

   using SuperClass::at;

   friend bool ROOT::Detail::VecOps::IsSmall<T>(const RVec<T> &v);

   friend bool ROOT::Detail::VecOps::IsAdopting<T>(const RVec<T> &v);
};

template <typename T, unsigned N>
inline size_t CapacityInBytes(const RVecN<T, N> &X)
{
   return X.capacity_in_bytes();
}

///@name RVec Unary Arithmetic Operators
///@{

#define RVEC_UNARY_OPERATOR(OP)                                                \
template <typename T>                                                          \
RVec<T> operator OP(const RVec<T> &v)                                          \
{                                                                              \
   RVec<T> ret(v);                                                             \
   for (auto &x : ret)                                                         \
      x = OP x;                                                                \
return ret;                                                                    \
}                                                                              \

RVEC_UNARY_OPERATOR(+)
RVEC_UNARY_OPERATOR(-)
RVEC_UNARY_OPERATOR(~)
RVEC_UNARY_OPERATOR(!)
#undef RVEC_UNARY_OPERATOR

///@}
///@name RVec Binary Arithmetic Operators
///@{

#define ERROR_MESSAGE(OP) \
 "Cannot call operator " #OP " on vectors of different sizes."

#define RVEC_BINARY_OPERATOR(OP)                                               \
template <typename T0, typename T1>                                            \
auto operator OP(const RVec<T0> &v, const T1 &y)                               \
  -> RVec<decltype(v[0] OP y)>                                                 \
{                                                                              \
   RVec<decltype(v[0] OP y)> ret(v.size());                                    \
   auto op = [&y](const T0 &x) { return x OP y; };                             \
   std::transform(v.begin(), v.end(), ret.begin(), op);                        \
   return ret;                                                                 \
}                                                                              \
                                                                               \
template <typename T0, typename T1>                                            \
auto operator OP(const T0 &x, const RVec<T1> &v)                               \
  -> RVec<decltype(x OP v[0])>                                                 \
{                                                                              \
   RVec<decltype(x OP v[0])> ret(v.size());                                    \
   auto op = [&x](const T1 &y) { return x OP y; };                             \
   std::transform(v.begin(), v.end(), ret.begin(), op);                        \
   return ret;                                                                 \
}                                                                              \
                                                                               \
template <typename T0, typename T1>                                            \
auto operator OP(const RVec<T0> &v0, const RVec<T1> &v1)                       \
  -> RVec<decltype(v0[0] OP v1[0])>                                            \
{                                                                              \
   if (v0.size() != v1.size())                                                 \
      throw std::runtime_error(ERROR_MESSAGE(OP));                             \
                                                                               \
   RVec<decltype(v0[0] OP v1[0])> ret(v0.size());                              \
   auto op = [](const T0 &x, const T1 &y) { return x OP y; };                  \
   std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), op);          \
   return ret;                                                                 \
}                                                                              \

RVEC_BINARY_OPERATOR(+)
RVEC_BINARY_OPERATOR(-)
RVEC_BINARY_OPERATOR(*)
RVEC_BINARY_OPERATOR(/)
RVEC_BINARY_OPERATOR(%)
RVEC_BINARY_OPERATOR(^)
RVEC_BINARY_OPERATOR(|)
RVEC_BINARY_OPERATOR(&)
#undef RVEC_BINARY_OPERATOR

///@}
///@name RVec Assignment Arithmetic Operators
///@{

#define RVEC_ASSIGNMENT_OPERATOR(OP)                                           \
template <typename T0, typename T1>                                            \
RVec<T0>& operator OP(RVec<T0> &v, const T1 &y)                                \
{                                                                              \
   auto op = [&y](T0 &x) { return x OP y; };                                   \
   std::transform(v.begin(), v.end(), v.begin(), op);                          \
   return v;                                                                   \
}                                                                              \
                                                                               \
template <typename T0, typename T1>                                            \
RVec<T0>& operator OP(RVec<T0> &v0, const RVec<T1> &v1)                        \
{                                                                              \
   if (v0.size() != v1.size())                                                 \
      throw std::runtime_error(ERROR_MESSAGE(OP));                             \
                                                                               \
   auto op = [](T0 &x, const T1 &y) { return x OP y; };                        \
   std::transform(v0.begin(), v0.end(), v1.begin(), v0.begin(), op);           \
   return v0;                                                                  \
}                                                                              \

RVEC_ASSIGNMENT_OPERATOR(+=)
RVEC_ASSIGNMENT_OPERATOR(-=)
RVEC_ASSIGNMENT_OPERATOR(*=)
RVEC_ASSIGNMENT_OPERATOR(/=)
RVEC_ASSIGNMENT_OPERATOR(%=)
RVEC_ASSIGNMENT_OPERATOR(^=)
RVEC_ASSIGNMENT_OPERATOR(|=)
RVEC_ASSIGNMENT_OPERATOR(&=)
RVEC_ASSIGNMENT_OPERATOR(>>=)
RVEC_ASSIGNMENT_OPERATOR(<<=)
#undef RVEC_ASSIGNMENT_OPERATOR

///@}
///@name RVec Comparison and Logical Operators
///@{

#define RVEC_LOGICAL_OPERATOR(OP)                                              \
template <typename T0, typename T1>                                            \
auto operator OP(const RVec<T0> &v, const T1 &y)                               \
  -> RVec<int> /* avoid std::vector<bool> */                                   \
{                                                                              \
   RVec<int> ret(v.size());                                                    \
   auto op = [y](const T0 &x) -> int { return x OP y; };                       \
   std::transform(v.begin(), v.end(), ret.begin(), op);                        \
   return ret;                                                                 \
}                                                                              \
                                                                               \
template <typename T0, typename T1>                                            \
auto operator OP(const T0 &x, const RVec<T1> &v)                               \
  -> RVec<int> /* avoid std::vector<bool> */                                   \
{                                                                              \
   RVec<int> ret(v.size());                                                    \
   auto op = [x](const T1 &y) -> int { return x OP y; };                       \
   std::transform(v.begin(), v.end(), ret.begin(), op);                        \
   return ret;                                                                 \
}                                                                              \
                                                                               \
template <typename T0, typename T1>                                            \
auto operator OP(const RVec<T0> &v0, const RVec<T1> &v1)                       \
  -> RVec<int> /* avoid std::vector<bool> */                                   \
{                                                                              \
   if (v0.size() != v1.size())                                                 \
      throw std::runtime_error(ERROR_MESSAGE(OP));                             \
                                                                               \
   RVec<int> ret(v0.size());                                                   \
   auto op = [](const T0 &x, const T1 &y) -> int { return x OP y; };           \
   std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), op);          \
   return ret;                                                                 \
}                                                                              \

RVEC_LOGICAL_OPERATOR(<)
RVEC_LOGICAL_OPERATOR(>)
RVEC_LOGICAL_OPERATOR(==)
RVEC_LOGICAL_OPERATOR(!=)
RVEC_LOGICAL_OPERATOR(<=)
RVEC_LOGICAL_OPERATOR(>=)
RVEC_LOGICAL_OPERATOR(&&)
RVEC_LOGICAL_OPERATOR(||)
#undef RVEC_LOGICAL_OPERATOR

///@}
///@name RVec Standard Mathematical Functions
///@{

/// \cond
template <typename T> struct PromoteTypeImpl;

template <> struct PromoteTypeImpl<float>       { using Type = float;       };
template <> struct PromoteTypeImpl<double>      { using Type = double;      };
template <> struct PromoteTypeImpl<long double> { using Type = long double; };

template <typename T> struct PromoteTypeImpl { using Type = double; };

template <typename T>
using PromoteType = typename PromoteTypeImpl<T>::Type;

template <typename U, typename V>
using PromoteTypes = decltype(PromoteType<U>() + PromoteType<V>());

/// \endcond

#define RVEC_UNARY_FUNCTION(NAME, FUNC)                                        \
   template <typename T>                                                       \
   RVec<PromoteType<T>> NAME(const RVec<T> &v)                                 \
   {                                                                           \
      RVec<PromoteType<T>> ret(v.size());                                      \
      auto f = [](const T &x) { return FUNC(x); };                             \
      std::transform(v.begin(), v.end(), ret.begin(), f);                      \
      return ret;                                                              \
   }

#define RVEC_BINARY_FUNCTION(NAME, FUNC)                                       \
   template <typename T0, typename T1>                                         \
   RVec<PromoteTypes<T0, T1>> NAME(const T0 &x, const RVec<T1> &v)             \
   {                                                                           \
      RVec<PromoteTypes<T0, T1>> ret(v.size());                                \
      auto f = [&x](const T1 &y) { return FUNC(x, y); };                       \
      std::transform(v.begin(), v.end(), ret.begin(), f);                      \
      return ret;                                                              \
   }                                                                           \
                                                                               \
   template <typename T0, typename T1>                                         \
   RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &v, const T1 &y)             \
   {                                                                           \
      RVec<PromoteTypes<T0, T1>> ret(v.size());                                \
      auto f = [&y](const T0 &x) { return FUNC(x, y); };                       \
      std::transform(v.begin(), v.end(), ret.begin(), f);                      \
      return ret;                                                              \
   }                                                                           \
                                                                               \
   template <typename T0, typename T1>                                         \
   RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &v0, const RVec<T1> &v1)     \
   {                                                                           \
      if (v0.size() != v1.size())                                              \
         throw std::runtime_error(ERROR_MESSAGE(NAME));                        \
                                                                               \
      RVec<PromoteTypes<T0, T1>> ret(v0.size());                               \
      auto f = [](const T0 &x, const T1 &y) { return FUNC(x, y); };            \
      std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), f);        \
      return ret;                                                              \
   }                                                                           \

#define RVEC_STD_UNARY_FUNCTION(F) RVEC_UNARY_FUNCTION(F, std::F)
#define RVEC_STD_BINARY_FUNCTION(F) RVEC_BINARY_FUNCTION(F, std::F)

RVEC_STD_UNARY_FUNCTION(abs)
RVEC_STD_BINARY_FUNCTION(fdim)
RVEC_STD_BINARY_FUNCTION(fmod)
RVEC_STD_BINARY_FUNCTION(remainder)

RVEC_STD_UNARY_FUNCTION(exp)
RVEC_STD_UNARY_FUNCTION(exp2)
RVEC_STD_UNARY_FUNCTION(expm1)

RVEC_STD_UNARY_FUNCTION(log)
RVEC_STD_UNARY_FUNCTION(log10)
RVEC_STD_UNARY_FUNCTION(log2)
RVEC_STD_UNARY_FUNCTION(log1p)

RVEC_STD_BINARY_FUNCTION(pow)
RVEC_STD_UNARY_FUNCTION(sqrt)
RVEC_STD_UNARY_FUNCTION(cbrt)
RVEC_STD_BINARY_FUNCTION(hypot)

RVEC_STD_UNARY_FUNCTION(sin)
RVEC_STD_UNARY_FUNCTION(cos)
RVEC_STD_UNARY_FUNCTION(tan)
RVEC_STD_UNARY_FUNCTION(asin)
RVEC_STD_UNARY_FUNCTION(acos)
RVEC_STD_UNARY_FUNCTION(atan)
RVEC_STD_BINARY_FUNCTION(atan2)

RVEC_STD_UNARY_FUNCTION(sinh)
RVEC_STD_UNARY_FUNCTION(cosh)
RVEC_STD_UNARY_FUNCTION(tanh)
RVEC_STD_UNARY_FUNCTION(asinh)
RVEC_STD_UNARY_FUNCTION(acosh)
RVEC_STD_UNARY_FUNCTION(atanh)

RVEC_STD_UNARY_FUNCTION(floor)
RVEC_STD_UNARY_FUNCTION(ceil)
RVEC_STD_UNARY_FUNCTION(trunc)
RVEC_STD_UNARY_FUNCTION(round)
RVEC_STD_UNARY_FUNCTION(lround)
RVEC_STD_UNARY_FUNCTION(llround)

RVEC_STD_UNARY_FUNCTION(erf)
RVEC_STD_UNARY_FUNCTION(erfc)
RVEC_STD_UNARY_FUNCTION(lgamma)
RVEC_STD_UNARY_FUNCTION(tgamma)
#undef RVEC_STD_UNARY_FUNCTION

///@}
///@name RVec Fast Mathematical Functions with Vdt
///@{

#ifdef R__HAS_VDT
#define RVEC_VDT_UNARY_FUNCTION(F) RVEC_UNARY_FUNCTION(F, vdt::F)

RVEC_VDT_UNARY_FUNCTION(fast_expf)
RVEC_VDT_UNARY_FUNCTION(fast_logf)
RVEC_VDT_UNARY_FUNCTION(fast_sinf)
RVEC_VDT_UNARY_FUNCTION(fast_cosf)
RVEC_VDT_UNARY_FUNCTION(fast_tanf)
RVEC_VDT_UNARY_FUNCTION(fast_asinf)
RVEC_VDT_UNARY_FUNCTION(fast_acosf)
RVEC_VDT_UNARY_FUNCTION(fast_atanf)

RVEC_VDT_UNARY_FUNCTION(fast_exp)
RVEC_VDT_UNARY_FUNCTION(fast_log)
RVEC_VDT_UNARY_FUNCTION(fast_sin)
RVEC_VDT_UNARY_FUNCTION(fast_cos)
RVEC_VDT_UNARY_FUNCTION(fast_tan)
RVEC_VDT_UNARY_FUNCTION(fast_asin)
RVEC_VDT_UNARY_FUNCTION(fast_acos)
RVEC_VDT_UNARY_FUNCTION(fast_atan)
#undef RVEC_VDT_UNARY_FUNCTION

#endif // R__HAS_VDT

#undef RVEC_UNARY_FUNCTION

///@}

/// Inner product
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVec<float> v1 {1., 2., 3.};
/// RVec<float> v2 {4., 5., 6.};
/// auto v1_dot_v2 = Dot(v1, v2);
/// v1_dot_v2
/// // (float) 32.0000f
/// ~~~
template <typename T, typename V>
auto Dot(const RVec<T> &v0, const RVec<V> &v1) -> decltype(v0[0] * v1[0])
{
   if (v0.size() != v1.size())
      throw std::runtime_error("Cannot compute inner product of vectors of different sizes");
   return std::inner_product(v0.begin(), v0.end(), v1.begin(), decltype(v0[0] * v1[0])(0));
}

/// Sum elements of an RVec
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 3.f};
/// auto v_sum = Sum(v);
/// v_sum
/// // (float) 6.f
/// auto v_sum_d = Sum(v, 0.);
/// v_sum_d
/// // (double) 6.0000000
/// ~~~
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// const ROOT::Math::PtEtaPhiMVector lv0 {15.5f, .3f, .1f, 105.65f},
///   lv1 {34.32f, 2.2f, 3.02f, 105.65f},
///   lv2 {12.95f, 1.32f, 2.2f, 105.65f};
/// RVec<ROOT::Math::PtEtaPhiMVector> v {lv0, lv1, lv2};
/// auto v_sum_lv = Sum(v, ROOT::Math::PtEtaPhiMVector());
/// v_sum_lv
/// // (ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > &) (30.8489,2.46534,2.58947,361.084)
/// ~~~
template <typename T>
T Sum(const RVec<T> &v, const T zero = T(0))
{
   return std::accumulate(v.begin(), v.end(), zero);
}

inline std::size_t Sum(const RVec<bool> &v, std::size_t zero = 0ul)
{
   return std::accumulate(v.begin(), v.end(), zero);
}

/// Return the product of the elements of the RVec.
template <typename T>
T Product(const RVec<T> &v, const T init = T(1)) // initialize with identity
{
   return std::accumulate(v.begin(), v.end(), init, std::multiplies<T>());
}

/// Get the mean of the elements of an RVec
///
/// The return type is a double precision floating point number.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_mean = Mean(v);
/// v_mean
/// // (double) 2.3333333
/// ~~~
template <typename T>
double Mean(const RVec<T> &v)
{
   if (v.empty()) return 0.;
   return double(Sum(v)) / v.size();
}

/// Get the mean of the elements of an RVec with custom initial value
///
/// The return type will be deduced from the `zero` parameter
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_mean_f = Mean(v, 0.f);
/// v_mean_f
/// // (float) 2.33333f
/// auto v_mean_d = Mean(v, 0.);
/// v_mean_d
/// // (double) 2.3333333
/// ~~~
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// const ROOT::Math::PtEtaPhiMVector lv0 {15.5f, .3f, .1f, 105.65f},
///   lv1 {34.32f, 2.2f, 3.02f, 105.65f},
///   lv2 {12.95f, 1.32f, 2.2f, 105.65f};
/// RVec<ROOT::Math::PtEtaPhiMVector> v {lv0, lv1, lv2};
/// auto v_mean_lv = Mean(v, ROOT::Math::PtEtaPhiMVector());
/// v_mean_lv
/// // (ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > &) (10.283,2.46534,2.58947,120.361)
/// ~~~
template <typename T, typename R = T>
R Mean(const RVec<T> &v, const R zero)
{
   if (v.empty()) return zero;
   return Sum(v, zero) / v.size();
}

/// Get the greatest element of an RVec
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_max = Max(v);
/// v_max
/// (float) 4.00000f
/// ~~~
template <typename T>
T Max(const RVec<T> &v)
{
   return *std::max_element(v.begin(), v.end());
}

/// Get the smallest element of an RVec
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_min = Min(v);
/// v_min
/// (float) 1.00000f
/// ~~~
template <typename T>
T Min(const RVec<T> &v)
{
   return *std::min_element(v.begin(), v.end());
}

/// Get the index of the greatest element of an RVec
/// In case of multiple occurrences of the maximum values,
/// the index corresponding to the first occurrence is returned.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_argmax = ArgMax(v);
/// v_argmax
/// // (unsigned long) 2
/// ~~~
template <typename T>
std::size_t ArgMax(const RVec<T> &v)
{
   return std::distance(v.begin(), std::max_element(v.begin(), v.end()));
}

/// Get the index of the smallest element of an RVec
/// In case of multiple occurrences of the minimum values,
/// the index corresponding to the first occurrence is returned.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_argmin = ArgMin(v);
/// v_argmin
/// // (unsigned long) 0
/// ~~~
template <typename T>
std::size_t ArgMin(const RVec<T> &v)
{
   return std::distance(v.begin(), std::min_element(v.begin(), v.end()));
}

/// Get the variance of the elements of an RVec
///
/// The return type is a double precision floating point number.
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_var = Var(v);
/// v_var
/// // (double) 2.3333333
/// ~~~
template <typename T>
double Var(const RVec<T> &v)
{
   const std::size_t size = v.size();
   if (size < std::size_t(2)) return 0.;
   T sum_squares(0), squared_sum(0);
   auto pred = [&sum_squares, &squared_sum](const T& x) {sum_squares+=x*x; squared_sum+=x;};
   std::for_each(v.begin(), v.end(), pred);
   squared_sum *= squared_sum;
   const auto dsize = (double) size;
   return 1. / (dsize - 1.) * (sum_squares - squared_sum / dsize );
}

/// Get the standard deviation of the elements of an RVec
///
/// The return type is a double precision floating point number.
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_sd = StdDev(v);
/// v_sd
/// // (double) 1.5275252
/// ~~~
template <typename T>
double StdDev(const RVec<T> &v)
{
   return std::sqrt(Var(v));
}

/// Create new collection applying a callable to the elements of the input collection
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v {1.f, 2.f, 4.f};
/// auto v_square = Map(v, [](float f){return f* 2.f;});
/// v_square
/// // (ROOT::VecOps::RVec<float> &) { 2.00000f, 4.00000f, 8.00000f }
///
/// RVecF x({1.f, 2.f, 3.f});
/// RVecF y({4.f, 5.f, 6.f});
/// RVecF z({7.f, 8.f, 9.f});
/// auto mod = [](float x, float y, float z) { return sqrt(x * x + y * y + z * z); };
/// auto v_mod = Map(x, y, z, mod);
/// v_mod
/// // (ROOT::VecOps::RVec<float> &) { 8.12404f, 9.64365f, 11.2250f }
/// ~~~
template <typename... Args>
auto Map(Args &&... args)
{
   /*
   Here the strategy in order to generalise the previous implementation of Map, i.e.
   `RVec Map(RVec, F)`, here we need to move the last parameter of the pack in first
   position in order to be able to invoke the Map function with automatic type deduction.
   This is achieved in two steps:
   1. Forward as tuple the pack to MapFromTuple
   2. Invoke the MapImpl helper which has the signature `template<...T, F> RVec MapImpl(F &&f, RVec<T>...)`
   */

   // check the first N - 1 arguments are RVecs
   constexpr auto nArgs = sizeof...(Args);
   constexpr bool isRVec[]{ROOT::Internal::VecOps::IsRVec<std::remove_cv_t<std::remove_reference_t<Args>>>::value...};
   static_assert(ROOT::Internal::VecOps::All(isRVec, nArgs - 1),
                 "Map: the first N-1 arguments must be RVecs or references to RVecs");

   return ROOT::Internal::VecOps::MapFromTuple(std::forward_as_tuple(args...),
                                               std::make_index_sequence<sizeof...(args) - 1>());
}

/// Create a new collection with the elements passing the filter expressed by the predicate
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecI v {1, 2, 4};
/// auto v_even = Filter(v, [](int i){return 0 == i%2;});
/// v_even
/// // (ROOT::VecOps::RVec<int> &) { 2, 4 }
/// ~~~
template <typename T, typename F>
RVec<T> Filter(const RVec<T> &v, F &&f)
{
   const auto thisSize = v.size();
   RVec<T> w;
   w.reserve(thisSize);
   for (auto &&val : v) {
      if (f(val))
         w.emplace_back(val);
   }
   return w;
}

/// Return true if any of the elements equates to true, return false otherwise.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecI v {0, 1, 0};
/// auto anyTrue = Any(v);
/// anyTrue
/// // (bool) true
/// ~~~
template <typename T>
auto Any(const RVec<T> &v) -> decltype(v[0] == true)
{
   for (auto &&e : v)
      if (static_cast<bool>(e) == true)
         return true;
   return false;
}

/// Return true if all of the elements equate to true, return false otherwise.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecI v {0, 1, 0};
/// auto allTrue = All(v);
/// allTrue
/// // (bool) false
/// ~~~
template <typename T>
auto All(const RVec<T> &v) -> decltype(v[0] == false)
{
   for (auto &&e : v)
      if (static_cast<bool>(e) == false)
         return false;
   return true;
}

template <typename T>
void swap(RVec<T> &lhs, RVec<T> &rhs)
{
   lhs.swap(rhs);
}

/// Return an RVec of indices that sort the input RVec
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto sortIndices = Argsort(v)
/// // (ROOT::VecOps::RVec<unsigned long> &) { 2, 0, 1 }
/// auto values = Take(v, sortIndices)
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<typename RVec<T>::size_type> Argsort(const RVec<T> &v)
{
   using size_type = typename RVec<T>::size_type;
   RVec<size_type> i(v.size());
   std::iota(i.begin(), i.end(), 0);
   std::sort(i.begin(), i.end(), [&v](size_type i1, size_type i2) { return v[i1] < v[i2]; });
   return i;
}

/// Return an RVec of indices that sort the input RVec based on a comparison function.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto sortIndices = Argsort(v, [](double x, double y) {return x > y;})
/// // (ROOT::VecOps::RVec<unsigned long> &) { 1, 0, 2 }
/// auto values = Take(v, sortIndices)
/// // (ROOT::VecOps::RVec<double> &) { 3.0000000, 2.0000000, 1.0000000 }
/// ~~~
template <typename T, typename Compare>
RVec<typename RVec<T>::size_type> Argsort(const RVec<T> &v, Compare &&c)
{
   using size_type = typename RVec<T>::size_type;
   RVec<size_type> i(v.size());
   std::iota(i.begin(), i.end(), 0);
   std::sort(i.begin(), i.end(),
             [&v, &c](size_type i1, size_type i2) { return c(v[i1], v[i2]); });
   return i;
}

/// Return an RVec of indices that sort the input RVec
/// while keeping the order of equal elements.
/// This is the stable variant of `Argsort`.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 2., 1.};
/// auto sortIndices = StableArgsort(v)
/// // (ROOT::VecOps::RVec<unsigned long> &) { 3, 0, 2, 1 }
/// auto values = Take(v, sortIndices)
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 2.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<typename RVec<T>::size_type> StableArgsort(const RVec<T> &v)
{
   using size_type = typename RVec<T>::size_type;
   RVec<size_type> i(v.size());
   std::iota(i.begin(), i.end(), 0);
   std::stable_sort(i.begin(), i.end(), [&v](size_type i1, size_type i2) { return v[i1] < v[i2]; });
   return i;
}

/// Return an RVec of indices that sort the input RVec based on a comparison function
/// while keeping the order of equal elements.
/// This is the stable variant of `Argsort`.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 2., 1.};
/// auto sortIndices = StableArgsort(v, [](double x, double y) {return x > y;})
/// // (ROOT::VecOps::RVec<unsigned long> &) { 1, 0, 2, 3 }
/// auto values = Take(v, sortIndices)
/// // (ROOT::VecOps::RVec<double> &) { 3.0000000, 2.0000000, 2.0000000, 1.0000000 }
/// ~~~
template <typename T, typename Compare>
RVec<typename RVec<T>::size_type> StableArgsort(const RVec<T> &v, Compare &&c)
{
   using size_type = typename RVec<T>::size_type;
   RVec<size_type> i(v.size());
   std::iota(i.begin(), i.end(), 0);
   std::stable_sort(i.begin(), i.end(), [&v, &c](size_type i1, size_type i2) { return c(v[i1], v[i2]); });
   return i;
}

/// Return elements of a vector at given indices
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto vTaken = Take(v, {0,2});
/// vTaken
/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 1.0000000 }
/// ~~~

template <typename T>
RVec<T> Take(const RVec<T> &v, const RVec<typename RVec<T>::size_type> &i)
{
   using size_type = typename RVec<T>::size_type;
   const size_type isize = i.size();
   RVec<T> r(isize);
   for (size_type k = 0; k < isize; k++)
      r[k] = v[i[k]];
   return r;
}

/// Take version that defaults to (user-specified) output value if some index is out of range
template <typename T>
RVec<T> Take(const RVec<T> &v, const RVec<typename RVec<T>::size_type> &i, const T default_val)
{
   using size_type = typename RVec<T>::size_type;
   const size_type isize = i.size();
   RVec<T> r(isize);
   for (size_type k = 0; k < isize; k++)
   {
      if (i[k] < v.size() && i[k]>=0){
         r[k] = v[i[k]];
      }
      else {
         r[k] = default_val;
      }
   }
   return r;
}

/// Return first `n` elements of an RVec if `n > 0` and last `n` elements if `n < 0`.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto firstTwo = Take(v, 2);
/// firstTwo
/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 3.0000000 }
/// auto lastOne = Take(v, -1);
/// lastOne
/// // (ROOT::VecOps::RVec<double>) { 1.0000000 }
/// ~~~
template <typename T>
RVec<T> Take(const RVec<T> &v, const int n)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = v.size();
   const size_type absn = std::abs(n);
   if (absn > size) {
      const auto msg = std::to_string(absn) + " elements requested from Take but input contains only " +
                       std::to_string(size) + " elements.";
      throw std::runtime_error(msg);
   }
   RVec<T> r(absn);
   if (n < 0) {
      for (size_type k = 0; k < absn; k++)
         r[k] = v[size - absn + k];
   } else {
      for (size_type k = 0; k < absn; k++)
         r[k] = v[k];
   }
   return r;
}

/// Return first `n` elements of an RVec if `n > 0` and last `n` elements if `n < 0`.
///
/// This Take version defaults to a user-specified value
/// `default_val` if the absolute value of `n` is
/// greater than the size of the RVec `v`
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using ROOT::VecOps::RVec;
/// RVec<int> x{1,2,3,4};
/// Take(x,-5,1)
/// // (ROOT::VecOps::RVec<int>) { 1, 1, 2, 3, 4 }
/// Take(x,5,20)
/// // (ROOT::VecOps::RVec<int>) { 1, 2, 3, 4, 20 }
/// Take(x,-1,1)
/// // (ROOT::VecOps::RVec<int>) { 4 }
/// Take(x,4,1)
/// // (ROOT::VecOps::RVec<int>) { 1, 2, 3, 4 }
/// ~~~
template <typename T>
RVec<T> Take(const RVec<T> &v, const int n, const T default_val)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = v.size();
   const size_type absn = std::abs(n);
   // Base case, can be handled by another overload of Take
   if (absn <= size) {
      return Take(v, n);
   }
   RVec<T> temp = v;
   // Case when n is positive and n > v.size()
   if (n > 0) {
      temp.resize(n, default_val);
      return temp;
   }
   // Case when n is negative and abs(n) > v.size()
   const auto num_to_fill = absn - size;
   ROOT::VecOps::RVec<T> fill_front(num_to_fill, default_val);
   return Concatenate(fill_front, temp);
}

/// Return a copy of the container without the elements at the specified indices.
///
/// Duplicated and out-of-range indices in idxs are ignored.
template <typename T>
RVec<T> Drop(const RVec<T> &v, RVec<typename RVec<T>::size_type> idxs)
{
   // clean up input indices
   std::sort(idxs.begin(), idxs.end());
   idxs.erase(std::unique(idxs.begin(), idxs.end()), idxs.end());

   RVec<T> r;
   if (v.size() > idxs.size())
      r.reserve(v.size() - idxs.size());

   auto discardIt = idxs.begin();
   using sz_t = typename RVec<T>::size_type;
   for (sz_t i = 0u; i < v.size(); ++i) {
      if (discardIt != idxs.end() && i == *discardIt)
         ++discardIt;
      else
         r.emplace_back(v[i]);
   }

   return r;
}

/// Return copy of reversed vector
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto v_reverse = Reverse(v);
/// v_reverse
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 3.0000000, 2.0000000 }
/// ~~~
template <typename T>
RVec<T> Reverse(const RVec<T> &v)
{
   RVec<T> r(v);
   std::reverse(r.begin(), r.end());
   return r;
}

/// Return copy of RVec with elements sorted in ascending order
///
/// This helper is different from Argsort since it does not return an RVec of indices,
/// but an RVec of values.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto v_sorted = Sort(v);
/// v_sorted
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<T> Sort(const RVec<T> &v)
{
   RVec<T> r(v);
   std::sort(r.begin(), r.end());
   return r;
}

/// Return copy of RVec with elements sorted based on a comparison operator
///
/// The comparison operator has to fulfill the same requirements of the
/// predicate of by std::sort.
///
///
/// This helper is different from Argsort since it does not return an RVec of indices,
/// but an RVec of values.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 1.};
/// auto v_sorted = Sort(v, [](double x, double y) {return 1/x < 1/y;});
/// v_sorted
/// // (ROOT::VecOps::RVec<double> &) { 3.0000000, 2.0000000, 1.0000000 }
/// ~~~
template <typename T, typename Compare>
RVec<T> Sort(const RVec<T> &v, Compare &&c)
{
   RVec<T> r(v);
   std::sort(r.begin(), r.end(), std::forward<Compare>(c));
   return r;
}

/// Return copy of RVec with elements sorted in ascending order
/// while keeping the order of equal elements.
///
/// This is the stable variant of `Sort`.
///
/// This helper is different from StableArgsort since it does not return an RVec of indices,
/// but an RVec of values.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 2, 1.};
/// auto v_sorted = StableSort(v);
/// v_sorted
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 2.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<T> StableSort(const RVec<T> &v)
{
   RVec<T> r(v);
   std::stable_sort(r.begin(), r.end());
   return r;
}

// clang-format off
/// Return copy of RVec with elements sorted based on a comparison operator
/// while keeping the order of equal elements.
///
/// The comparison operator has to fulfill the same requirements of the
/// predicate of std::stable_sort.
///
/// This helper is different from StableArgsort since it does not return an RVec of indices,
/// but an RVec of values.
///
/// This is the stable variant of `Sort`.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 3., 2., 1.};
/// auto v_sorted = StableSort(v, [](double x, double y) {return 1/x < 1/y;});
/// v_sorted
/// // (ROOT::VecOps::RVec<double> &) { 3.0000000, 2.0000000, 2.0000000, 1.0000000 }
/// ~~~
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVec<RVecD> v {{2., 4.}, {3., 1.}, {2, 1.}, {1., 4.}};
/// auto v_sorted = StableSort(StableSort(v, [](const RVecD &x, const RVecD &y) {return x[1] < y[1];}), [](const RVecD &x, const RVecD &y) {return x[0] < y[0];});
/// v_sorted
/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<double> > &) { { 1.0000000, 4.0000000 }, { 2.0000000, 1.0000000 }, { 2.0000000, 4.0000000 }, { 3.0000000, 1.0000000 } }
/// ~~~
// clang-format off
template <typename T, typename Compare>
RVec<T> StableSort(const RVec<T> &v, Compare &&c)
{
   RVec<T> r(v);
   std::stable_sort(r.begin(), r.end(), std::forward<Compare>(c));
   return r;
}

/// Return the indices that represent all combinations of the elements of two
/// RVecs.
///
/// The type of the return value is an RVec of two RVecs containing indices.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// auto comb_idx = Combinations(3, 2);
/// comb_idx
/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0, 0, 1, 1, 2, 2 }, { 0, 1, 0, 1, 0, 1 } }
/// ~~~
inline RVec<RVec<std::size_t>> Combinations(const std::size_t size1, const std::size_t size2)
{
   using size_type = std::size_t;
   RVec<RVec<size_type>> r(2);
   r[0].resize(size1*size2);
   r[1].resize(size1*size2);
   size_type c = 0;
   for(size_type i=0; i<size1; i++) {
      for(size_type j=0; j<size2; j++) {
         r[0][c] = i;
         r[1][c] = j;
         c++;
      }
   }
   return r;
}

/// Return the indices that represent all combinations of the elements of two
/// RVecs.
///
/// The type of the return value is an RVec of two RVecs containing indices.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v1 {1., 2., 3.};
/// RVecD v2 {-4., -5.};
/// auto comb_idx = Combinations(v1, v2);
/// comb_idx
/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0, 0, 1, 1, 2, 2 }, { 0, 1, 0, 1, 0, 1 } }
/// ~~~
template <typename T1, typename T2>
RVec<RVec<typename RVec<T1>::size_type>> Combinations(const RVec<T1> &v1, const RVec<T2> &v2)
{
   return Combinations(v1.size(), v2.size());
}

/// Return the indices that represent all unique combinations of the
/// elements of a given RVec.
///
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {1., 2., 3., 4.};
/// auto v_1 = Combinations(v, 1);
/// v_1
/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0, 1, 2, 3 } }
/// auto v_2 = Combinations(v, 2);
/// v_2
/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0, 0, 0, 1, 1, 2 }, { 1, 2, 3, 2, 3, 3 } }
/// auto v_3 = Combinations(v, 3);
/// v_3
/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0, 0, 0, 1 }, { 1, 1, 2, 2 }, { 2, 3, 3, 3 } }
/// auto v_4 = Combinations(v, 4);
/// v_4
/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<unsigned long> > &) { { 0 }, { 1 }, { 2 }, { 3 } }
/// ~~~
template <typename T>
RVec<RVec<typename RVec<T>::size_type>> Combinations(const RVec<T>& v, const typename RVec<T>::size_type n)
{
   using size_type = typename RVec<T>::size_type;
   const size_type s = v.size();
   if (n > s) {
      throw std::runtime_error("Cannot make unique combinations of size " + std::to_string(n) +
                               " from vector of size " + std::to_string(s) + ".");
   }

   RVec<size_type> indices(s);
   for(size_type k=0; k<s; k++)
      indices[k] = k;

   const auto innersize = [=] {
      size_type inners = s - n + 1;
      for (size_type m = s - n + 2; m <= s; ++m)
         inners *= m;

      size_type factn = 1;
      for (size_type i = 2; i <= n; ++i)
         factn *= i;
      inners /= factn;

      return inners;
   }();

   RVec<RVec<size_type>> c(n, RVec<size_type>(innersize));
   size_type inneridx = 0;
   for (size_type k = 0; k < n; k++)
      c[k][inneridx] = indices[k];
   ++inneridx;

   while (true) {
      bool run_through = true;
      long i = n - 1;
      for (; i>=0; i--) {
         if (indices[i] != i + s - n){
            run_through = false;
            break;
         }
      }
      if (run_through) {
         return c;
      }
      indices[i]++;
      for (long j=i+1; j<(long)n; j++)
         indices[j] = indices[j-1] + 1;
      for (size_type k = 0; k < n; k++)
         c[k][inneridx] = indices[k];
      ++inneridx;
   }
}

/// Return the indices of the elements which are not zero
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v {2., 0., 3., 0., 1.};
/// auto nonzero_idx = Nonzero(v);
/// nonzero_idx
/// // (ROOT::VecOps::RVec<unsigned long> &) { 0, 2, 4 }
/// ~~~
template <typename T>
RVec<typename RVec<T>::size_type> Nonzero(const RVec<T> &v)
{
   using size_type = typename RVec<T>::size_type;
   RVec<size_type> r;
   const auto size = v.size();
   r.reserve(size);
   for(size_type i=0; i<size; i++) {
      if(v[i] != 0) {
         r.emplace_back(i);
      }
   }
   return r;
}

/// Return the intersection of elements of two RVecs.
///
/// Each element of v1 is looked up in v2 and added to the returned vector if
/// found. Following, the order of v1 is preserved. If v2 is already sorted, the
/// optional argument v2_is_sorted can be used to toggle of the internal sorting
/// step, therewith optimising runtime.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v1 {1., 2., 3.};
/// RVecD v2 {-4., -5., 2., 1.};
/// auto v1_intersect_v2 = Intersect(v1, v2);
/// v1_intersect_v2
/// // (ROOT::VecOps::RVec<double> &) { 1.0000000, 2.0000000 }
/// ~~~
template <typename T>
RVec<T> Intersect(const RVec<T>& v1, const RVec<T>& v2, bool v2_is_sorted = false)
{
   RVec<T> v2_sorted;
   if (!v2_is_sorted) v2_sorted = Sort(v2);
   const auto v2_begin = v2_is_sorted ? v2.begin() : v2_sorted.begin();
   const auto v2_end = v2_is_sorted ? v2.end() : v2_sorted.end();
   RVec<T> r;
   const auto size = v1.size();
   r.reserve(size);
   using size_type = typename RVec<T>::size_type;
   for(size_type i=0; i<size; i++) {
      if (std::binary_search(v2_begin, v2_end, v1[i])) {
         r.emplace_back(v1[i]);
      }
   }
   return r;
}

/// Return the elements of v1 if the condition c is true and v2 if the
/// condition c is false.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v1 {1., 2., 3.};
/// RVecD v2 {-1., -2., -3.};
/// auto c = v1 > 1;
/// c
/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
/// auto if_c_v1_else_v2 = Where(c, v1, v2);
/// if_c_v1_else_v2
/// // (ROOT::VecOps::RVec<double> &) { -1.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<T> Where(const RVec<int>& c, const RVec<T>& v1, const RVec<T>& v2)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = c.size();
   RVec<T> r;
   r.reserve(size);
   for (size_type i=0; i<size; i++) {
      r.emplace_back(c[i] != 0 ? v1[i] : v2[i]);
   }
   return r;
}

/// Return the elements of v1 if the condition c is true and sets the value v2
/// if the condition c is false.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecD v1 {1., 2., 3.};
/// double v2 = 4.;
/// auto c = v1 > 1;
/// c
/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
/// auto if_c_v1_else_v2 = Where(c, v1, v2);
/// if_c_v1_else_v2
/// // (ROOT::VecOps::RVec<double>) { 4.0000000, 2.0000000, 3.0000000 }
/// ~~~
template <typename T>
RVec<T> Where(const RVec<int> &c, const RVec<T> &v1, typename RVec<T>::value_type v2)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = c.size();
   RVec<T> r;
   r.reserve(size);
   for (size_type i=0; i<size; i++) {
      r.emplace_back(c[i] != 0 ? v1[i] : v2);
   }
   return r;
}

/// Return the elements of v2 if the condition c is false and sets the value v1
/// if the condition c is true.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// double v1 = 4.;
/// RVecD v2 {1., 2., 3.};
/// auto c = v2 > 1;
/// c
/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
/// auto if_c_v1_else_v2 = Where(c, v1, v2);
/// if_c_v1_else_v2
/// // (ROOT::VecOps::RVec<double>) { 1.0000000, 4.0000000, 4.0000000 }
/// ~~~
template <typename T>
RVec<T> Where(const RVec<int>& c, typename RVec<T>::value_type v1, const RVec<T>& v2)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = c.size();
   RVec<T> r;
   r.reserve(size);
   for (size_type i=0; i<size; i++) {
      r.emplace_back(c[i] != 0 ? v1 : v2[i]);
   }
   return r;
}

/// Return a vector with the value v2 if the condition c is false and sets the
/// value v1 if the condition c is true.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// double v1 = 4.;
/// double v2 = 2.;
/// RVecI c {0, 1, 1};
/// auto if_c_v1_else_v2 = Where(c, v1, v2);
/// if_c_v1_else_v2
/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 4.0000000, 4.0000000 }
/// ~~~
template <typename T>
RVec<T> Where(const RVec<int>& c, T v1, T v2)
{
   using size_type = typename RVec<T>::size_type;
   const size_type size = c.size();
   RVec<T> r;
   r.reserve(size);
   for (size_type i=0; i<size; i++) {
      r.emplace_back(c[i] != 0 ? v1 : v2);
   }
   return r;
}

/// Return the concatenation of two RVecs.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF rvf {0.f, 1.f, 2.f};
/// RVecI rvi {7, 8, 9};
/// Concatenate(rvf, rvi)
/// // (ROOT::VecOps::RVec<float>) { 0.00000f, 1.00000f, 2.00000f, 7.00000f, 8.00000f, 9.00000f }
/// ~~~
template <typename T0, typename T1, typename Common_t = typename std::common_type<T0, T1>::type>
RVec<Common_t> Concatenate(const RVec<T0> &v0, const RVec<T1> &v1)
{
   RVec<Common_t> res;
   res.reserve(v0.size() + v1.size());
   std::copy(v0.begin(), v0.end(), std::back_inserter(res));
   std::copy(v1.begin(), v1.end(), std::back_inserter(res));
   return res;
}

/// Return the angle difference \f$\Delta \phi\f$ of two scalars.
///
/// The function computes the closest angle from v1 to v2 with sign and is
/// therefore in the range \f$[-\pi, \pi]\f$.
/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
/// to degrees \f$c = 180\f$.
template <typename T0, typename T1 = T0, typename Common_t = std::common_type_t<T0, T1>>
Common_t DeltaPhi(T0 v1, T1 v2, const Common_t c = M_PI)
{
   static_assert(std::is_floating_point<T0>::value && std::is_floating_point<T1>::value,
                 "DeltaPhi must be called with floating point values.");
   auto r = std::fmod(v2 - v1, 2.0 * c);
   if (r < -c) {
      r += 2.0 * c;
   }
   else if (r > c) {
      r -= 2.0 * c;
   }
   return r;
}

/// Return the angle difference \f$\Delta \phi\f$ in radians of two vectors.
///
/// The function computes the closest angle from v1 to v2 with sign and is
/// therefore in the range \f$[-\pi, \pi]\f$.
/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
/// to degrees \f$c = 180\f$.
template <typename T0, typename T1 = T0, typename Common_t = typename std::common_type_t<T0, T1>>
RVec<Common_t> DeltaPhi(const RVec<T0>& v1, const RVec<T1>& v2, const Common_t c = M_PI)
{
   using size_type = typename RVec<T0>::size_type;
   const size_type size = v1.size();
   auto r = RVec<Common_t>(size);
   for (size_type i = 0; i < size; i++) {
      r[i] = DeltaPhi(v1[i], v2[i], c);
   }
   return r;
}

/// Return the angle difference \f$\Delta \phi\f$ in radians of a vector and a scalar.
///
/// The function computes the closest angle from v1 to v2 with sign and is
/// therefore in the range \f$[-\pi, \pi]\f$.
/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
/// to degrees \f$c = 180\f$.
template <typename T0, typename T1 = T0, typename Common_t = typename std::common_type_t<T0, T1>>
RVec<Common_t> DeltaPhi(const RVec<T0>& v1, T1 v2, const Common_t c = M_PI)
{
   using size_type = typename RVec<T0>::size_type;
   const size_type size = v1.size();
   auto r = RVec<Common_t>(size);
   for (size_type i = 0; i < size; i++) {
      r[i] = DeltaPhi(v1[i], v2, c);
   }
   return r;
}

/// Return the angle difference \f$\Delta \phi\f$ in radians of a scalar and a vector.
///
/// The function computes the closest angle from v1 to v2 with sign and is
/// therefore in the range \f$[-\pi, \pi]\f$.
/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
/// to degrees \f$c = 180\f$.
template <typename T0, typename T1 = T0, typename Common_t = typename std::common_type_t<T0, T1>>
RVec<Common_t> DeltaPhi(T0 v1, const RVec<T1>& v2, const Common_t c = M_PI)
{
   using size_type = typename RVec<T1>::size_type;
   const size_type size = v2.size();
   auto r = RVec<Common_t>(size);
   for (size_type i = 0; i < size; i++) {
      r[i] = DeltaPhi(v1, v2[i], c);
   }
   return r;
}

/// Return the square of the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
/// the collections eta1, eta2, phi1 and phi2.
///
/// The function computes \f$\Delta R^2 = (\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2\f$
/// of the given collections eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
/// be set to radian or degrees using the optional argument c, see the documentation
/// of the DeltaPhi helper.
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3>>
RVec<Common_t> DeltaR2(const RVec<T0>& eta1, const RVec<T1>& eta2, const RVec<T2>& phi1, const RVec<T3>& phi2, const Common_t c = M_PI)
{
   const auto dphi = DeltaPhi(phi1, phi2, c);
   return (eta1 - eta2) * (eta1 - eta2) + dphi * dphi;
}

/// Return the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
/// the collections eta1, eta2, phi1 and phi2.
///
/// The function computes \f$\Delta R = \sqrt{(\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2}\f$
/// of the given collections eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
/// be set to radian or degrees using the optional argument c, see the documentation
/// of the DeltaPhi helper.
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3>>
RVec<Common_t> DeltaR(const RVec<T0>& eta1, const RVec<T1>& eta2, const RVec<T2>& phi1, const RVec<T3>& phi2, const Common_t c = M_PI)
{
   return sqrt(DeltaR2(eta1, eta2, phi1, phi2, c));
}

/// Return the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
/// the scalars eta1, eta2, phi1 and phi2.
///
/// The function computes \f$\Delta R = \sqrt{(\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2}\f$
/// of the given scalars eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
/// be set to radian or degrees using the optional argument c, see the documentation
/// of the DeltaPhi helper.
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3>>
Common_t DeltaR(T0 eta1, T1 eta2, T2 phi1, T3 phi2, const Common_t c = M_PI)
{
   const auto dphi = DeltaPhi(phi1, phi2, c);
   return std::sqrt((eta1 - eta2) * (eta1 - eta2) + dphi * dphi);
}

/// Return the angle between two three-vectors given the quantities
/// x coordinate (x), y coordinate (y), z coordinate (y).
///
/// The function computes the angle between two three-vectors
/// (x1, y2, z1) and (x2, y2, z2).
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename T4 = T0,
          typename T5 = T0, typename Common_t = std::common_type_t<T0, T1>>
Common_t Angle(T0 x1, T1 y1, T2 z1, T3 x2, T4 y2, T5 z2){
    // cross product
    const auto cx = y1 * z2 - y2 * z1;
    const auto cy = x1 * z2 - x2 * z1;
    const auto cz = x1 * y2 - x2 * y1;
    
    // norm of cross product
    const auto c = std::sqrt(cx * cx + cy * cy + cz * cz);
    
    // dot product
    const auto  d = x1 * x2 + y1 * y2 + z1 * z2;
    
    return std::atan2(c, d);
}

/// Return the invariant mass of two particles given
/// x coordinate (px), y coordinate (py), z coordinate (pz) and mass.
///
/// The function computes the invariant mass of two particles with the four-vectors
/// (x1, y2, z1, mass1) and (x2, py2, pz2, mass2).
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename T4 = T0,
          typename T5 = T0, typename T6 = T0, typename T7 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3, T4, T5, T6, T7>>
Common_t InvariantMasses_PxPyPzM(
   const T0& x1, const T1& y1, const T2& z1, const T3& mass1,
   const T4& x2, const T5& y2, const T6& z2, const T7& mass2)
{

      // Numerically stable computation of Invariant Masses
      const auto p1_sq = x1 * x1 + y1 * y1 + z1 * z1;
      const auto p2_sq = x2 * x2 + y2 * y2 + z2 * z2;

      if (p1_sq <= 0 && p2_sq <= 0)
         return (mass1 + mass2);
      if (p1_sq <= 0) {
         auto mm = mass1 + std::sqrt(mass2*mass2 + p2_sq);
         auto m2 = mm*mm - p2_sq; 
         if (m2 >= 0)
            return std::sqrt( m2 );
         else
            return std::sqrt( -m2 );
      }
      if (p2_sq <= 0) {
         auto mm = mass2 + std::sqrt(mass1*mass1 + p1_sq);
         auto m2 = mm*mm - p1_sq; 
         if (m2 >= 0)
            return std::sqrt( m2 );
         else
            return std::sqrt( -m2 );
      }

      const auto m1_sq =  mass1 * mass1;
      const auto m2_sq =  mass2 * mass2;

      const auto r1 = m1_sq / p1_sq;
      const auto r2 = m2_sq / p2_sq;
      const auto x = r1 + r2 + r1 * r2;
      const auto a = Angle(x1, y1, z1, x2, y2, z2);
      const auto cos_a = std::cos(a);
      auto y = x;
      if ( cos_a >= 0){
         y = (x + std::sin(a) * std::sin(a)) / (std::sqrt(x + 1) + cos_a);
      } else {
         y = std::sqrt(x + 1) - cos_a;
      }

      const auto z = 2 * std::sqrt(p1_sq * p2_sq);

   // Return invariant mass with (+, -, -, -) metric
   return std::sqrt(m1_sq + m2_sq + y * z);
}

/// Return the invariant mass of two particles given the collections of the quantities
/// x coordinate (px), y coordinate (py), z coordinate (pz) and mass.
///
/// The function computes the invariant mass of two particles with the four-vectors
/// (px1, py2, pz1, mass1) and (px2, py2, pz2, mass2).
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename T4 = T0,
          typename T5 = T0, typename T6 = T0, typename T7 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3, T4, T5, T6, T7>>
RVec<Common_t> InvariantMasses_PxPyPzM(
   const RVec<T0>& px1, const RVec<T1>& py1, const RVec<T2>& pz1, const RVec<T3>& mass1,
   const RVec<T4>& px2, const RVec<T5>& py2, const RVec<T6>& pz2, const RVec<T7>& mass2)
{
   std::size_t size = px1.size();

   R__ASSERT(py1.size() == size && pz1.size() == size && mass1.size() == size);
   R__ASSERT(px2.size() == size && py2.size() == size && pz2.size() == size && mass2.size() == size);

   RVec<Common_t> inv_masses(size);

   for (std::size_t i = 0u; i < size; ++i) {
      inv_masses[i] = InvariantMasses_PxPyPzM(px1[i], py1[i], pz1[i], mass1[i], px2[i], py2[i], pz2[i], mass2[i]);
   }

   // Return invariant mass with (+, -, -, -) metric
   return inv_masses;
}

/// Return the invariant mass of two particles given the collections of the quantities
/// transverse momentum (pt), rapidity (eta), azimuth (phi) and mass.
///
/// The function computes the invariant mass of two particles with the four-vectors
/// (pt1, eta2, phi1, mass1) and (pt2, eta2, phi2, mass2).
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename T4 = T0,
          typename T5 = T0, typename T6 = T0, typename T7 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3, T4, T5, T6, T7>>
RVec<Common_t> InvariantMasses(
   const RVec<T0>& pt1, const RVec<T1>& eta1, const RVec<T2>& phi1, const RVec<T3>& mass1,
   const RVec<T4>& pt2, const RVec<T5>& eta2, const RVec<T6>& phi2, const RVec<T7>& mass2)
{
   std::size_t size = pt1.size();

   R__ASSERT(eta1.size() == size && phi1.size() == size && mass1.size() == size);
   R__ASSERT(pt2.size() == size && phi2.size() == size && mass2.size() == size);

   RVec<Common_t> inv_masses(size);

   for (std::size_t i = 0u; i < size; ++i) {
      // Conversion from (pt, eta, phi, mass) to (x, y, z, mass) coordinate system
      const auto x1 = pt1[i] * std::cos(phi1[i]);
      const auto y1 = pt1[i] * std::sin(phi1[i]);
      const auto z1 = pt1[i] * std::sinh(eta1[i]);
      
      const auto x2 = pt2[i] * std::cos(phi2[i]);
      const auto y2 = pt2[i] * std::sin(phi2[i]);
      const auto z2 = pt2[i] * std::sinh(eta2[i]);
      
      // Numerically stable computation of Invariant Masses
      inv_masses[i] = InvariantMasses_PxPyPzM(x1, y1, z1, mass1[i], x2, y2, z2, mass2[i]);
   }

   // Return invariant mass with (+, -, -, -) metric
   return inv_masses;
}

/// Return the invariant mass of multiple particles given the collections of the
/// quantities transverse momentum (pt), rapidity (eta), azimuth (phi) and mass.
///
/// The function computes the invariant mass of multiple particles with the
/// four-vectors (pt, eta, phi, mass).
template <typename T0, typename T1 = T0, typename T2 = T0, typename T3 = T0, typename Common_t = std::common_type_t<T0, T1, T2, T3>>
Common_t InvariantMass(const RVec<T0>& pt, const RVec<T1>& eta, const RVec<T2>& phi, const RVec<T3>& mass)
{
   const std::size_t size = pt.size();

   R__ASSERT(eta.size() == size && phi.size() == size && mass.size() == size);

   Common_t x_sum = 0.;
   Common_t y_sum = 0.;
   Common_t z_sum = 0.;
   Common_t e_sum = 0.;

   for (std::size_t i = 0u; i < size; ++ i) {
      // Convert to (e, x, y, z) coordinate system and update sums
      const auto x = pt[i] * std::cos(phi[i]);
      x_sum += x;
      const auto y = pt[i] * std::sin(phi[i]);
      y_sum += y;
      const auto z = pt[i] * std::sinh(eta[i]);
      z_sum += z;
      const auto e = std::sqrt(x * x + y * y + z * z + mass[i] * mass[i]);
      e_sum += e;
   }

   // Return invariant mass with (+, -, -, -) metric
   return std::sqrt(e_sum * e_sum - x_sum * x_sum - y_sum * y_sum - z_sum * z_sum);
}

////////////////////////////////////////////////////////////////////////////
/// \brief Build an RVec of objects starting from RVecs of input to their constructors.
/// \tparam T Type of the objects contained in the created RVec.
/// \tparam Args_t Pack of types templating the input RVecs.
/// \param[in] args The RVecs containing the values used to initialise the output objects.
/// \return The RVec of objects initialised with the input parameters.
///
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF pts = {15.5, 34.32, 12.95};
/// RVecF etas = {0.3, 2.2, 1.32};
/// RVecF phis = {0.1, 3.02, 2.2};
/// RVecF masses = {105.65, 105.65, 105.65};
/// auto fourVecs = Construct<ROOT::Math::PtEtaPhiMVector>(pts, etas, phis, masses);
/// cout << fourVecs << endl;
/// // { (15.5,0.3,0.1,105.65), (34.32,2.2,3.02,105.65), (12.95,1.32,2.2,105.65) }
/// ~~~
template <typename T, typename... Args_t>
RVec<T> Construct(const RVec<Args_t> &... args)
{
   const auto size = ::ROOT::Internal::VecOps::GetVectorsSize("Construct", args...);
   RVec<T> ret;
   ret.reserve(size);
   for (auto i = 0UL; i < size; ++i) {
      ret.emplace_back(args[i]...);
   }
   return ret;
}

/// For any Rvec v produce another RVec with entries starting from 0, and incrementing by 1 until a N = v.size() is reached.
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// RVecF v = {1., 2., 3.};
/// cout << Enumerate(v1) << "\n";
/// // { 0, 1, 2 }
/// ~~~
template <typename T>
RVec<typename RVec<T>::size_type> Enumerate(const RVec<T> &v)
{
   const auto size = v.size();
   RVec<T> ret;
   ret.reserve(size);
   for (auto i = 0UL; i < size; ++i) {
      ret.emplace_back(i);
   }
   return ret;
}

/// Produce RVec with entries starting from 0, and incrementing by 1 until a user-specified N is reached.
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// cout << Range(3) << "\n";
/// // { 0, 1, 2 }
/// ~~~
inline RVec<std::size_t> Range(std::size_t length)
{
   RVec<std::size_t> ret;
   ret.reserve(length);
   for (auto i = 0UL; i < length; ++i) {
      ret.emplace_back(i);
   }
   return ret;
}

/// Produce RVec with entries equal to begin, begin+1, ..., end-1.
/// An empty RVec is returned if begin >= end.
inline RVec<std::size_t> Range(std::size_t begin, std::size_t end)
{
   RVec<std::size_t> ret;
   ret.reserve(begin < end ? end - begin : 0u);
   for (auto i = begin; i < end; ++i)
      ret.push_back(i);
   return ret;
}

/// Allows for negative begin, end, and/or stride. Produce RVec<int> with entries equal to begin, begin+stride, ... , N,
/// where N is the first integer such that N+stride exceeds or equals N in the positive or negative direction (same as in Python).
/// An empty RVec is returned if begin >= end and stride > 0 or if
/// begin < end and stride < 0. Throws a runtime_error if stride==0
/// Example code, at the ROOT prompt:
/// ~~~{.cpp}
/// using namespace ROOT::VecOps;
/// cout << Range(1, 5, 2) << "\n";
/// // { 1, 3 }
/// cout << Range(-1, -11, -4) << "\n";
/// // { -1, -5, -9 }
/// ~~~
inline RVec<long long int> Range(long long int begin, long long int end, long long int stride)
{
   if (stride==0ll)
   {
      throw std::runtime_error("Range: the stride must not be zero");
   }
   RVec<long long int> ret;
   float ret_cap = std::ceil(static_cast<float>(end-begin) / stride); //the capacity to reserve
   //ret_cap < 0 if either begin > end & stride > 0, or begin < end & stride < 0. In both cases, an empty RVec should be returned
   if (ret_cap < 0)
   {
      return ret;
   }
   ret.reserve(static_cast<size_t>(ret_cap));
   if (stride > 0)
   {
      for (auto i = begin; i < end; i+=stride)
         ret.push_back(i);
   }
   else
   {
      for (auto i = begin; i > end; i+=stride)
         ret.push_back(i);
   }
   return ret;
}

////////////////////////////////////////////////////////////////////////////////
/// Print a RVec at the prompt:
template <class T>
std::ostream &operator<<(std::ostream &os, const RVec<T> &v)
{
   // In order to print properly, convert to 64 bit int if this is a char
   constexpr bool mustConvert = std::is_same<char, T>::value || std::is_same<signed char, T>::value ||
                                std::is_same<unsigned char, T>::value || std::is_same<wchar_t, T>::value ||
                                std::is_same<char16_t, T>::value || std::is_same<char32_t, T>::value;
   using Print_t = typename std::conditional<mustConvert, long long int, T>::type;
   os << "{ ";
   auto size = v.size();
   if (size) {
      for (std::size_t i = 0; i < size - 1; ++i) {
         os << (Print_t)v[i] << ", ";
      }
      os << (Print_t)v[size - 1];
   }
   os << " }";
   return os;
}

#if (_VECOPS_USE_EXTERN_TEMPLATES)

#define RVEC_EXTERN_UNARY_OPERATOR(T, OP) \
   extern template RVec<T> operator OP<T>(const RVec<T> &);

#define RVEC_EXTERN_BINARY_OPERATOR(T, OP)                                     \
   extern template auto operator OP<T, T>(const T &x, const RVec<T> &v)        \
      -> RVec<decltype(x OP v[0])>;                                            \
   extern template auto operator OP<T, T>(const RVec<T> &v, const T &y)        \
      -> RVec<decltype(v[0] OP y)>;                                            \
   extern template auto operator OP<T, T>(const RVec<T> &v0, const RVec<T> &v1)\
      -> RVec<decltype(v0[0] OP v1[0])>;

#define RVEC_EXTERN_ASSIGN_OPERATOR(T, OP)                           \
   extern template RVec<T> &operator OP<T, T>(RVec<T> &, const T &); \
   extern template RVec<T> &operator OP<T, T>(RVec<T> &, const RVec<T> &);

#define RVEC_EXTERN_LOGICAL_OPERATOR(T, OP)                                 \
   extern template RVec<int> operator OP<T, T>(const RVec<T> &, const T &); \
   extern template RVec<int> operator OP<T, T>(const T &, const RVec<T> &); \
   extern template RVec<int> operator OP<T, T>(const RVec<T> &, const RVec<T> &);

#define RVEC_EXTERN_FLOAT_TEMPLATE(T)   \
   extern template class RVec<T>;       \
   RVEC_EXTERN_UNARY_OPERATOR(T, +)     \
   RVEC_EXTERN_UNARY_OPERATOR(T, -)     \
   RVEC_EXTERN_UNARY_OPERATOR(T, !)     \
   RVEC_EXTERN_BINARY_OPERATOR(T, +)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, -)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, *)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, /)    \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, +=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, -=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, *=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, /=)   \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, <)   \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, >)   \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, ==)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, !=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, <=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, >=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, &&)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, ||)

#define RVEC_EXTERN_INTEGER_TEMPLATE(T) \
   extern template class RVec<T>;       \
   RVEC_EXTERN_UNARY_OPERATOR(T, +)     \
   RVEC_EXTERN_UNARY_OPERATOR(T, -)     \
   RVEC_EXTERN_UNARY_OPERATOR(T, ~)     \
   RVEC_EXTERN_UNARY_OPERATOR(T, !)     \
   RVEC_EXTERN_BINARY_OPERATOR(T, +)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, -)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, *)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, /)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, %)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, &)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, |)    \
   RVEC_EXTERN_BINARY_OPERATOR(T, ^)    \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, +=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, -=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, *=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, /=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, %=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, &=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, |=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, ^=)   \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, >>=)  \
   RVEC_EXTERN_ASSIGN_OPERATOR(T, <<=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, <)   \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, >)   \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, ==)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, !=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, <=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, >=)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, &&)  \
   RVEC_EXTERN_LOGICAL_OPERATOR(T, ||)

RVEC_EXTERN_INTEGER_TEMPLATE(char)
RVEC_EXTERN_INTEGER_TEMPLATE(short)
RVEC_EXTERN_INTEGER_TEMPLATE(int)
RVEC_EXTERN_INTEGER_TEMPLATE(long)
//RVEC_EXTERN_INTEGER_TEMPLATE(long long)

RVEC_EXTERN_INTEGER_TEMPLATE(unsigned char)
RVEC_EXTERN_INTEGER_TEMPLATE(unsigned short)
RVEC_EXTERN_INTEGER_TEMPLATE(unsigned int)
RVEC_EXTERN_INTEGER_TEMPLATE(unsigned long)
//RVEC_EXTERN_INTEGER_TEMPLATE(unsigned long long)

RVEC_EXTERN_FLOAT_TEMPLATE(float)
RVEC_EXTERN_FLOAT_TEMPLATE(double)

#undef RVEC_EXTERN_UNARY_OPERATOR
#undef RVEC_EXTERN_BINARY_OPERATOR
#undef RVEC_EXTERN_ASSIGN_OPERATOR
#undef RVEC_EXTERN_LOGICAL_OPERATOR
#undef RVEC_EXTERN_INTEGER_TEMPLATE
#undef RVEC_EXTERN_FLOAT_TEMPLATE

#define RVEC_EXTERN_UNARY_FUNCTION(T, NAME, FUNC) \
   extern template RVec<PromoteType<T>> NAME(const RVec<T> &);

#define RVEC_EXTERN_STD_UNARY_FUNCTION(T, F) RVEC_EXTERN_UNARY_FUNCTION(T, F, std::F)

#define RVEC_EXTERN_BINARY_FUNCTION(T0, T1, NAME, FUNC)                            \
   extern template RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &, const T1 &); \
   extern template RVec<PromoteTypes<T0, T1>> NAME(const T0 &, const RVec<T1> &); \
   extern template RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &, const RVec<T1> &);

#define RVEC_EXTERN_STD_BINARY_FUNCTION(T, F) RVEC_EXTERN_BINARY_FUNCTION(T, T, F, std::F)

#define RVEC_EXTERN_STD_FUNCTIONS(T)             \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, abs)        \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, fdim)      \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, fmod)      \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, remainder) \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, exp)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, exp2)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, expm1)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, log)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, log10)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, log2)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, log1p)      \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, pow)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, sqrt)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, cbrt)       \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, hypot)     \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, sin)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, cos)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, tan)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, asin)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, acos)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, atan)       \
   RVEC_EXTERN_STD_BINARY_FUNCTION(T, atan2)     \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, sinh)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, cosh)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, tanh)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, asinh)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, acosh)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, atanh)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, floor)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, ceil)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, trunc)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, round)      \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, erf)        \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, erfc)       \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, lgamma)     \
   RVEC_EXTERN_STD_UNARY_FUNCTION(T, tgamma)     \

RVEC_EXTERN_STD_FUNCTIONS(float)
RVEC_EXTERN_STD_FUNCTIONS(double)
#undef RVEC_EXTERN_STD_UNARY_FUNCTION
#undef RVEC_EXTERN_STD_BINARY_FUNCTION
#undef RVEC_EXTERN_STD_UNARY_FUNCTIONS

#ifdef R__HAS_VDT

#define RVEC_EXTERN_VDT_UNARY_FUNCTION(T, F) RVEC_EXTERN_UNARY_FUNCTION(T, F, vdt::F)

RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_expf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_logf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_sinf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_cosf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_tanf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_asinf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_acosf)
RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_atanf)

RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_exp)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_log)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_sin)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_cos)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_tan)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_asin)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_acos)
RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_atan)

#endif // R__HAS_VDT

#endif // _VECOPS_USE_EXTERN_TEMPLATES

/** @} */ // end of Doxygen group vecops

} // End of VecOps NS

// Allow to use RVec as ROOT::RVec
using ROOT::VecOps::RVec;

using RVecB = ROOT::VecOps::RVec<bool>;
using RVecC = ROOT::VecOps::RVec<char>;
using RVecD = ROOT::VecOps::RVec<double>;
using RVecF = ROOT::VecOps::RVec<float>;
using RVecI = ROOT::VecOps::RVec<int>;
using RVecL = ROOT::VecOps::RVec<long int>;
using RVecLL = ROOT::VecOps::RVec<long long int>;
using RVecU = ROOT::VecOps::RVec<unsigned int>;
using RVecUL = ROOT::VecOps::RVec<unsigned long int>;
using RVecULL = ROOT::VecOps::RVec<unsigned long long int>;

} // End of ROOT NS

#endif // ROOT_RVEC