//+------------------------------------------------------------------+
//|                                                    RadixSort.mqh |
//|                                        Copyright © 2018, Amr Ali |
//|                             https://www.mql5.com/en/users/amrali |
//+------------------------------------------------------------------+

//+------------------------------------------------------------------+
//| RadixSort                                                        |
//+------------------------------------------------------------------+
//| Sorts the values in the first dimension of a multidimensional    |
//| numeric array in the ascending order.                            |
//|                                                                  |
//| Arguments:                                                       |
//| array[]                                                          |
//|   [in][out] : Numeric array for sorting.                         |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//|                                                                  |
//| Note:                                                            |
//|   An array is always sorted in the ascending order irrespective  |
//|   of the AS_SERIES flag value.                                   |
//|                                                                  |
//|   The function accepts any-dimensional arrays of simple type     |
//|   (char, uchar, short, ushort, int, uint, long, ulong, bool,     |
//|   color, datetime, float, double) as a parameter.  However,      |
//|   sorting is always applied to the first (zero) dimension.       |
//|                                                                  |
//| Performance:                                                     |
//|   This is a highly-optimized implementation of LSD RadixSort     |
//|   in MQL using radix 256 (8-bits). This could be useful for      |
//|   sorting huge numeric arrays with millions of numbers.          |
//|   It is at least 3-10 times faster than built-in ArraySort().    |
//+------------------------------------------------------------------+
template<typename T>
bool RadixSort(T &arr[])
  {
   int n = ArraySize(arr);
//--- check array size
   if(n < 2)
      return(n == 1);
//--- fall back to ArraySort() for a small array
   if(n < 128)
      return(ArraySort(arr));
//--- temporary array to hold the sorted numbers
   T temp[];
   if(ArrayResize(temp, n) != n)
      return(false);
//--- frequncy histogram of the i-th byte
   int counts[sizeof(T)][256];
//--- reset counts array
   ZeroMemory(counts);
//--- compute frequency histograms of the i-th byte
   for(int i = 0; i < n; i++)
     {
      T ai = arr[i];
      for(uchar col = 0; col < sizeof(T); col++)
        {
         uchar key = extract_key(ai, (uchar)(col << 3));
         ++counts[col][key];
        }
     }
//--- checks whether the array is dynamic, in order to utilize ArraySwap() optimization on MT5.
   bool IsDynamic = ArrayIsDynamic(arr);

//--- sort the array elements one byte at a time in LSD order (right-most byte first)
   for(uchar col = 0; col < sizeof(T); col++)
     {
      uchar shift = col << 3;

      //--- determine if any columns can be skipped
      //--- for each digit, check if all the elements are in the same bucket.
      //--- If so, we can skip the whole digit. Instead of checking all the buckets,
      //--- we pick first key and check whether the bucket contains all the elements.
      if(counts[col][extract_key(arr[0], shift)] == n)
        {
         continue;
        }
      //--- transform counts to offset (ranks)
      int offset = 0;
      for(int i = 0; i < 256; i++)
        {
         int old_count = counts[col][i];
         counts[col][i] = offset;
         offset += old_count;
        }
      //--- gather forwards in temp (stable sort, elements ordered by i-th byte)
      for(int i = 0; i < n; i++)
        {
         uchar key = extract_key(arr[i], shift);
         temp[counts[col][key]++] = arr[i];
        }
      //--- copy back from temp to arr
#ifdef __MQL5__
      if(IsDynamic)
         ArraySwap(arr, temp);
      else
         ArrayCopy(arr, temp);
#else
      ArrayCopy(arr, temp);
#endif
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| RadixSort                                                        |
//+------------------------------------------------------------------+
//| Sorts the values in the first dimension of a multidimensional    |
//| numeric array in the ascending order.                            |
//|                                                                  |
//| Function overload for two-dimensional numeric arrays.            |
//|                                                                  |
//| Arguments:                                                       |
//| array[]                                                          |
//|   [in][out] : Numeric array for sorting.                         |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool RadixSort(T &arr[][])
  {
   int n = ArraySize(arr);
   int na = ArrayRange(arr, 0);
   int nb = ArrayRange(arr, 1);
//--- check array size
   if(na < 2)
      return(na == 1);
//--- fall back to ArraySort() for a small array
   if(na < 128)
      return(ArraySort(arr));
//--- prepare temporary array to hold numbers in the first dimension
   T temp[];
   if(ArrayResize(temp, na) != na)
      return(false);
   for(int i = 0; i < na; i++)
      temp[i] = arr[i][0];
//--- calculate order of numbers in the first dimension
   int indices[];
   if(!RadixSortIndices(temp, indices))
      return(false);
//--- flatten the multidimensional numeric array.
   if(ArrayCopy(temp, arr) != n)
      return(false);
//--- sort the flattened multidimensional array using the sorted indices.
   for(int i = 0; i < na; i++)
     {
      int index = indices[i];
      for(int j = 0; j < nb; j++)
         arr[i][j] = temp[nb * index + j];
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| RadixSort                                                        |
//+------------------------------------------------------------------+
//| Sorts the values in the first dimension of a multidimensional    |
//| numeric array in the ascending order.                            |
//|                                                                  |
//| Function overload for three-dimensional numeric arrays.          |
//|                                                                  |
//| Arguments:                                                       |
//| array[]                                                          |
//|   [in][out] : Numeric array for sorting.                         |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool RadixSort(T &arr[][][])
  {
   int n = ArraySize(arr);
   int na = ArrayRange(arr, 0);
   int nb = ArrayRange(arr, 1);
   int nc = ArrayRange(arr, 2);
//--- check array size
   if(na < 2)
      return(na == 1);
//--- fall back to ArraySort() for a small array
   if(na < 128)
      return(ArraySort(arr));
//--- prepare temporary array to hold numbers in the first dimension
   T temp[];
   if(ArrayResize(temp, na) != na)
      return(false);
   for(int i = 0; i < na; i++)
      temp[i] = arr[i][0][0];

//--- calculate order of numbers in the first dimension
   int indices[];
   if(!RadixSortIndices(temp, indices))
      return(false);
//--- flatten the multidimensional numeric array.
   if(ArrayCopy(temp, arr) != n)
      return(false);
//--- sort the flattened multidimensional array using the sorted indices.
   for(int i = 0; i < na; i++)
     {
      int index = indices[i];
      for(int j = 0; j < nb; j++)
         for(int k = 0; k < nc; k++)
            arr[i][j][k] = temp[nc * (nb * index + j) + k];
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| RadixSort                                                        |
//+------------------------------------------------------------------+
//| Sorts the values in the first dimension of a multidimensional    |
//| numeric array in the ascending order.                            |
//|                                                                  |
//| Function overload for four-dimensional numeric arrays.           |
//|                                                                  |
//| Arguments:                                                       |
//| array[]                                                          |
//|   [in][out] : Numeric array for sorting.                         |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool RadixSort(T &arr[][][][])
  {
   int n = ArraySize(arr);
   int na = ArrayRange(arr, 0);
   int nb = ArrayRange(arr, 1);
   int nc = ArrayRange(arr, 2);
   int nd = ArrayRange(arr, 3);
//--- check array size
   if(na < 2)
      return(na == 1);
//--- fall back to ArraySort() for a small array
   if(na < 128)
      return(ArraySort(arr));
//--- prepare temporary array to hold numbers in the first dimension
   T temp[];
   if(ArrayResize(temp, na) != na)
      return(false);
   for(int i = 0; i < na; i++)
      temp[i] = arr[i][0][0][0];
//--- calculate order of numbers in the first dimension
   int indices[];
   if(!RadixSortIndices(temp, indices))
      return(false);
//--- flatten the multidimensional numeric array.
   if(ArrayCopy(temp, arr) != n)
      return(false);
//--- sort the flattened multidimensional array using the sorted indices.
   for(int i = 0; i < na; i++)
     {
      int index = indices[i];
      for(int j = 0; j < nb; j++)
         for(int k = 0; k < nc; k++)
            for(int l = 0; l < nd; l++)
               arr[i][j][k][l] = temp[nd * (nc * (nb * index + j) + k) + l];
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| RadixSortIndices                                                 |
//+------------------------------------------------------------------+
//| Populate an array of indices[] in the order that would sort      |
//| the values of a numeric array[] in the ascending order.          |
//|                                                                  |
//| Arguments:                                                       |
//| array[]     : Array with numeric values to sort                  |
//| indices[]   : Array for sorted indices                           |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T>
bool RadixSortIndices(const T &arr[], int &indices[])
  {
   int n = ArraySize(arr);
//--- prepare array to hold the indices of the numbers
   if(ArrayResize(indices, n) != n)
      return(false);
   for(int i = 0; i < n; i++)
      indices[i] = i;
//--- prepare temporary array to hold the numbers
   T temp[];
   if(ArrayCopy(temp, arr) != n)
      return(false);
//--- calculate order of numbers in the ascending order
   if(!ParallelRadixSort(temp, indices))
      return(false);
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| ParallelRadixSort                                                |
//+------------------------------------------------------------------+
//| The function sorts array[] and items[] simultaneously using      |
//| RadixSort algorithm. After sorting the items[] array is sorted   |
//| by the numeric values of array[] in the ascending order.         |
//|                                                                  |
//| Arguments:                                                       |
//| array[]     : Array with numeric keys to sort                    |
//| items[]     : Array for items to sort based on keys              |
//|                                                                  |
//| Return value: true if successful, otherwise false.               |
//+------------------------------------------------------------------+
template<typename T,typename TItem>
bool ParallelRadixSort(T &arr[], TItem &items[])
  {
   int n = ArraySize(arr);
//--- check array size
   if(n < 2)
      return(n == 1);
   if(ArraySize(items) != n)
      return(false);
//--- temporary array to hold the sorted numbers (keys)
   T temp[];
   if(ArrayResize(temp, n) != n)
      return(false);
//--- temporary array to hold the sorted items
   TItem items2[];
   if(ArrayResize(items2, n) != n)
      return(false);
//--- frequncy histogram of the i-th byte
   int counts[sizeof(T)][256];
//--- reset counts array
   ZeroMemory(counts);
//--- compute frequency histograms of the i-th byte
   for(int i = 0; i < n; i++)
     {
      T ai = arr[i];
      for(uchar col = 0; col < sizeof(T); col++)
        {
         uchar key = extract_key(ai, (uchar)(col << 3));
         ++counts[col][key];
        }
     }
//--- checks whether the arrays are dynamic, in order to utilize ArraySwap() optimization on MT5.
   bool IsDynamic = ArrayIsDynamic(arr) && ArrayIsDynamic(items);

//--- sort the array elements one byte at a time in LSD order (right-most byte first)
   for(uchar col = 0; col < sizeof(T); col++)
     {
      uchar shift = col << 3;

      //--- determine if any columns can be skipped
      //--- for each digit, check if all the elements are in the same bucket.
      //--- If so, we can skip the whole digit. Instead of checking all the buckets,
      //--- we pick first key and check whether the bucket contains all the elements.
      if(counts[col][extract_key(arr[0], shift)] == n)
        {
         continue;
        }
      //--- transform counts to offset (ranks)
      int offset = 0;
      for(int i = 0; i < 256; i++)
        {
         int old_count = counts[col][i];
         counts[col][i] = offset;
         offset += old_count;
        }
      //--- gather forwards in temp (stable sort, elements ordered by i-th byte)
      for(int i = 0; i < n; i++)
        {
         uchar key = extract_key(arr[i], shift);
         offset = counts[col][key]++;
         temp[offset] = arr[i];
         items2[offset] = items[i];
        }
      //--- copy back from temp to arr
#ifdef __MQL5__
      if(IsDynamic)
        {
         ArraySwap(arr, temp);
         ArraySwap(items, items2);
        }
      else
        {
         ArrayCopy(arr, temp);
         ArrayCopy(items, items2);
        }
#else
      ArrayCopy(arr, temp);
      ArrayCopy(items, items2);
#endif
     }
//---
   return(true);
  }
//+------------------------------------------------------------------+
//| Returns an integer key for "signed" primitive types.             |
//+------------------------------------------------------------------+
// https://github.com/eloj/radix-sorting#-key-derivation
// http://stereopsis.com/radix.html
// https://stackoverflow.com/a/42304235/4208440
// https://github.com/skarupke/ska_sort/blob/master/ska_sort.hpp
// https://github.com/JakubValtar/radsort/blob/master/src/scalar.rs
uchar extract_key(const char value, const uchar shift)
  {
   return (uchar)((value ^ 0x80) >> shift);
  }
uchar extract_key(const short value, const uchar shift)
  {
   return (uchar)((value ^ 0x8000) >> shift);
  }
uchar extract_key(const int value, const uchar shift)
  {
   return (uchar)((value ^ 0x80000000) >> shift);
  }
uchar extract_key(const long value, const uchar shift)
  {
   return (uchar)((value ^ 0x8000000000000000) >> shift);
  }
//+------------------------------------------------------------------+
//| Returns an integer key for float.                                |
//|                                                                  |
//| flip a float for sorting                                         |
//|  finds SIGN of fp number.                                        |
//|  if it's 1 (negative float), it flips all bits                   |
//|  if it's 0 (positive float), it flips the sign only              |
//+------------------------------------------------------------------+
uchar extract_key(const float value, const uchar shift)
  {
   union _f { float value; uint bits; } f;
   f.value = value;

//--- extend the sign bit to the whole width with arithmetic
//--- right shift to get a flip mask 0xffffffff or 0x80000000
   return (uchar)((f.bits ^ (((int)f.bits >> 31) | 0x80000000)) >> shift);
  }
//+------------------------------------------------------------------+
//| Returns an integer key for double.                               |
//+------------------------------------------------------------------+
uchar extract_key(const double value, const uchar shift)
  {
   union _d { double value; ulong bits; } d;
   d.value = value;

//--- extend the sign bit to the whole width with arithmetic
//--- right shift to get a flip mask 0xffffffff or 0x80000000
   return (uchar)((d.bits ^ (((long)d.bits >> 63) | 0x8000000000000000)) >> shift);
  }
//+------------------------------------------------------------------+
//| Returns an integer key for string.                               |
//+------------------------------------------------------------------+
uchar extract_key(const string value, const uchar shift)
  {
//--- not implemented.
   return(0);
  }
//+------------------------------------------------------------------+
//| Returns an integer key for other "unsigned" primitive types.     |
//| uchar, ushort, uint, ulong, bool, color, datetime.               |
//+------------------------------------------------------------------+
template<typename T>
uchar extract_key(const T value, const uchar shift)
  {
   return (uchar)(value >> shift);
  }
//+------------------------------------------------------------------+