VecGeom/base/Algorithms.h

0001 #ifndef VECGEOM_ALGORITHMS_H
0002 #define VECGEOM_ALGORITHMS_H
0003
0004 #include <cmath>
0005 #include <limits>
0006 #include "VecCore/Limits.h"
0007
0008 namespace vecgeom {
0009 namespace algo {
0010
0011 ///< A variant of the quicksort algorithm, avoiding recursiveness by using a limited-size
0012 ///< stack as 'max_levels'. Upon exhausing the stack sorting can fail, so the return value must be
0013 ///< checked. The maximum levels needed for sorting an array of size `len` is approximately:
0014 ///<     max_levels ~ 3 * log10(len)
0015 ///< returns success
0016 template <typename Type, int max_levels = 32>
0017 VECCORE_ATT_HOST_DEVICE
0018 bool quickSort(Type const *arr, size_t elements, size_t *sorted)
0019 {
0020   ///< Implementation credits go to: https://stackoverflow.com/a/55011578
0021
0022   size_t beg[max_levels], end[max_levels], L, R;
0023   int i = 0;
0024
0025   for (size_t j = 0; j < elements; ++j)
0026     sorted[j] = j;
0027
0028   beg[0] = 0;
0029   end[0] = elements;
0030   while (i >= 0) {
0031     L = beg[i];
0032     R = end[i];
0033     if (R - L > 1) {
0034       size_t M   = L + ((R - L) >> 1);
0035       size_t piv = sorted[M];
0036       sorted[M]  = sorted[L];
0037       if (i == max_levels - 1) return false;
0038       R--;
0039       while (L < R) {
0040         while (arr[sorted[R]] >= arr[piv] && L < R)
0041           R--;
0042         if (L < R) sorted[L++] = sorted[R];
0043         while (arr[sorted[L]] <= arr[piv] && L < R)
0044           L++;
0045         if (L < R) sorted[R--] = sorted[L];
0046       }
0047       sorted[L] = piv;
0048       M         = L + 1;
0049       while (L > beg[i] && arr[sorted[L - 1]] == arr[piv])
0050         L--;
0051       while (M < end[i] && arr[sorted[M]] == arr[piv])
0052         M++;
0053       if (L - beg[i] > end[i] - M) {
0054         beg[i + 1] = M;
0055         end[i + 1] = end[i];
0056         end[i++]   = L;
0057       } else {
0058         beg[i + 1] = beg[i];
0059         end[i + 1] = L;
0060         beg[i++]   = M;
0061       }
0062     } else {
0063       i--;
0064     }
0065   }
0066   return true;
0067 }
0068
0069 } // end namespace algo
0070 } // end namespace vecgeom
0071
0072 #endif