ksort.h

   1 /* The MIT License
   2
   3    Copyright (c) 2008 Genome Research Ltd (GRL).
   4
   5    Permission is hereby granted, free of charge, to any person obtaining
   6    a copy of this software and associated documentation files (the
   7    "Software"), to deal in the Software without restriction, including
   8    without limitation the rights to use, copy, modify, merge, publish,
   9    distribute, sublicense, and/or sell copies of the Software, and to
  10    permit persons to whom the Software is furnished to do so, subject to
  11    the following conditions:
  12
  13    The above copyright notice and this permission notice shall be
  14    included in all copies or substantial portions of the Software.
  15
  16    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  18    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  20    BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  21    ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  22    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23    SOFTWARE.
  24 */
  25
  26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
  27
  28 /*
  29   2012-12-11 (0.1.4):
  30
  31     * Defined __ks_insertsort_##name as static to compile with C99.
  32
  33   2008-11-16 (0.1.4):
  34
  35     * Fixed a bug in introsort() that happens in rare cases.
  36
  37   2008-11-05 (0.1.3):
  38
  39     * Fixed a bug in introsort() for complex comparisons.
  40
  41         * Fixed a bug in mergesort(). The previous version is not stable.
  42
  43   2008-09-15 (0.1.2):
  44
  45         * Accelerated introsort. On my Mac (not on another Linux machine),
  46           my implementation is as fast as std::sort on random input.
  47
  48         * Added combsort and in introsort, switch to combsort if the
  49           recursion is too deep.
  50
  51   2008-09-13 (0.1.1):
  52
  53         * Added k-small algorithm
  54
  55   2008-09-05 (0.1.0):
  56
  57         * Initial version
  58
  59 */
  60
  61 #ifndef AC_KSORT_H
  62 #define AC_KSORT_H
  63
  64 #include <stdlib.h>
  65 #include <string.h>
  66
  67 typedef struct {
  68         void *left, *right;
  69         int depth;
  70 } ks_isort_stack_t;
  71
  72 #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
  73
  74 #define KSORT_INIT(name, type_t, __sort_lt)                                                             \
  75         void ks_mergesort_##name(size_t n, type_t array[], type_t temp[])       \
  76         {                                                                                                                                       \
  77                 type_t *a2[2], *a, *b;                                                                                  \
  78                 int curr, shift;                                                                                                \
  79                                                                                                                                                 \
  80                 a2[0] = array;                                                                                                  \
  81                 a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n);               \
  82                 for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) {                  \
  83                         a = a2[curr]; b = a2[1-curr];                                                           \
  84                         if (shift == 0) {                                                                                       \
  85                                 type_t *p = b, *i, *eb = a + n;                                                 \
  86                                 for (i = a; i < eb; i += 2) {                                                   \
  87                                         if (i == eb - 1) *p++ = *i;                                                     \
  88                                         else {                                                                                          \
  89                                                 if (__sort_lt(*(i+1), *i)) {                                    \
  90                                                         *p++ = *(i+1); *p++ = *i;                                       \
  91                                                 } else {                                                                                \
  92                                                         *p++ = *i; *p++ = *(i+1);                                       \
  93                                                 }                                                                                               \
  94                                         }                                                                                                       \
  95                                 }                                                                                                               \
  96                         } else {                                                                                                        \
  97                                 size_t i, step = 1ul<<shift;                                                    \
  98                                 for (i = 0; i < n; i += step<<1) {                                              \
  99                                         type_t *p, *j, *k, *ea, *eb;                                            \
 100                                         if (n < i + step) {                                                                     \
 101                                                 ea = a + n; eb = a;                                                             \
 102                                         } else {                                                                                        \
 103                                                 ea = a + i + step;                                                              \
 104                                                 eb = a + (n < i + (step<<1)? n : i + (step<<1)); \
 105                                         }                                                                                                       \
 106                                         j = a + i; k = a + i + step; p = b + i;                         \
 107                                         while (j < ea && k < eb) {                                                      \
 108                                                 if (__sort_lt(*k, *j)) *p++ = *k++;                             \
 109                                                 else *p++ = *j++;                                                               \
 110                                         }                                                                                                       \
 111                                         while (j < ea) *p++ = *j++;                                                     \
 112                                         while (k < eb) *p++ = *k++;                                                     \
 113                                 }                                                                                                               \
 114                         }                                                                                                                       \
 115                         curr = 1 - curr;                                                                                        \
 116                 }                                                                                                                               \
 117                 if (curr == 1) {                                                                                                \
 118                         type_t *p = a2[0], *i = a2[1], *eb = array + n;                         \
 119                         for (; p < eb; ++i) *p++ = *i;                                                          \
 120                 }                                                                                                                               \
 121                 if (temp == 0) free(a2[1]);                                                                             \
 122         }                                                                                                                                       \
 123         void ks_heapadjust_##name(size_t i, size_t n, type_t l[])                       \
 124         {                                                                                                                                       \
 125                 size_t k = i;                                                                                                   \
 126                 type_t tmp = l[i];                                                                                              \
 127                 while ((k = (k << 1) + 1) < n) {                                                                \
 128                         if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k;                         \
 129                         if (__sort_lt(l[k], tmp)) break;                                                        \
 130                         l[i] = l[k]; i = k;                                                                                     \
 131                 }                                                                                                                               \
 132                 l[i] = tmp;                                                                                                             \
 133         }                                                                                                                                       \
 134         void ks_heapmake_##name(size_t lsize, type_t l[])                                       \
 135         {                                                                                                                                       \
 136                 size_t i;                                                                                                               \
 137                 for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i)                              \
 138                         ks_heapadjust_##name(i, lsize, l);                                                      \
 139         }                                                                                                                                       \
 140         void ks_heapsort_##name(size_t lsize, type_t l[])                                       \
 141         {                                                                                                                                       \
 142                 size_t i;                                                                                                               \
 143                 for (i = lsize - 1; i > 0; --i) {                                                               \
 144                         type_t tmp;                                                                                                     \
 145                         tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
 146                 }                                                                                                                               \
 147         }                                                                                                                                       \
 148         static inline void __ks_insertsort_##name(type_t *s, type_t *t)         \
 149         {                                                                                                                                       \
 150                 type_t *i, *j, swap_tmp;                                                                                \
 151                 for (i = s + 1; i < t; ++i)                                                                             \
 152                         for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) {                      \
 153                                 swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp;                  \
 154                         }                                                                                                                       \
 155         }                                                                                                                                       \
 156         void ks_combsort_##name(size_t n, type_t a[])                                           \
 157         {                                                                                                                                       \
 158                 const double shrink_factor = 1.2473309501039786540366528676643; \
 159                 int do_swap;                                                                                                    \
 160                 size_t gap = n;                                                                                                 \
 161                 type_t tmp, *i, *j;                                                                                             \
 162                 do {                                                                                                                    \
 163                         if (gap > 2) {                                                                                          \
 164                                 gap = (size_t)(gap / shrink_factor);                                    \
 165                                 if (gap == 9 || gap == 10) gap = 11;                                    \
 166                         }                                                                                                                       \
 167                         do_swap = 0;                                                                                            \
 168                         for (i = a; i < a + n - gap; ++i) {                                                     \
 169                                 j = i + gap;                                                                                    \
 170                                 if (__sort_lt(*j, *i)) {                                                                \
 171                                         tmp = *i; *i = *j; *j = tmp;                                            \
 172                                         do_swap = 1;                                                                            \
 173                                 }                                                                                                               \
 174                         }                                                                                                                       \
 175                 } while (do_swap || gap > 2);                                                                   \
 176                 if (gap != 1) __ks_insertsort_##name(a, a + n);                                 \
 177         }                                                                                                                                       \
 178         void ks_introsort_##name(size_t n, type_t a[])                                          \
 179         {                                                                                                                                       \
 180                 int d;                                                                                                                  \
 181                 ks_isort_stack_t *top, *stack;                                                                  \
 182                 type_t rp, swap_tmp;                                                                                    \
 183                 type_t *s, *t, *i, *j, *k;                                                                              \
 184                                                                                                                                                 \
 185                 if (n < 1) return;                                                                                              \
 186                 else if (n == 2) {                                                                                              \
 187                         if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
 188                         return;                                                                                                         \
 189                 }                                                                                                                               \
 190                 for (d = 2; 1ul<<d < n; ++d);                                                                   \
 191                 stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
 192                 top = stack; s = a; t = a + (n-1); d <<= 1;                                             \
 193                 while (1) {                                                                                                             \
 194                         if (s < t) {                                                                                            \
 195                                 if (--d == 0) {                                                                                 \
 196                                         ks_combsort_##name(t - s + 1, s);                                       \
 197                                         t = s;                                                                                          \
 198                                         continue;                                                                                       \
 199                                 }                                                                                                               \
 200                                 i = s; j = t; k = i + ((j-i)>>1) + 1;                                   \
 201                                 if (__sort_lt(*k, *i)) {                                                                \
 202                                         if (__sort_lt(*k, *j)) k = j;                                           \
 203                                 } else k = __sort_lt(*j, *i)? i : j;                                    \
 204                                 rp = *k;                                                                                                \
 205                                 if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; }  \
 206                                 for (;;) {                                                                                              \
 207                                         do ++i; while (__sort_lt(*i, rp));                                      \
 208                                         do --j; while (i <= j && __sort_lt(rp, *j));            \
 209                                         if (j <= i) break;                                                                      \
 210                                         swap_tmp = *i; *i = *j; *j = swap_tmp;                          \
 211                                 }                                                                                                               \
 212                                 swap_tmp = *i; *i = *t; *t = swap_tmp;                                  \
 213                                 if (i-s > t-i) {                                                                                \
 214                                         if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
 215                                         s = t-i > 16? i+1 : t;                                                          \
 216                                 } else {                                                                                                \
 217                                         if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
 218                                         t = i-s > 16? i-1 : s;                                                          \
 219                                 }                                                                                                               \
 220                         } else {                                                                                                        \
 221                                 if (top == stack) {                                                                             \
 222                                         free(stack);                                                                            \
 223                                         __ks_insertsort_##name(a, a+n);                                         \
 224                                         return;                                                                                         \
 225                                 } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
 226                         }                                                                                                                       \
 227                 }                                                                                                                               \
 228         }                                                                                                                                       \
 229         /* This function is adapted from: http://ndevilla.free.fr/median/ */ \
 230         /* 0 <= kk < n */                                                                                                       \
 231         type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk)                      \
 232         {                                                                                                                                       \
 233                 type_t *low, *high, *k, *ll, *hh, *mid;                                                 \
 234                 low = arr; high = arr + n - 1; k = arr + kk;                                    \
 235                 for (;;) {                                                                                                              \
 236                         if (high <= low) return *k;                                                                     \
 237                         if (high == low + 1) {                                                                          \
 238                                 if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
 239                                 return *k;                                                                                              \
 240                         }                                                                                                                       \
 241                         mid = low + (high - low) / 2;                                                           \
 242                         if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \
 243                         if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
 244                         if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low);      \
 245                         KSORT_SWAP(type_t, *mid, *(low+1));                                                     \
 246                         ll = low + 1; hh = high;                                                                        \
 247                         for (;;) {                                                                                                      \
 248                                 do ++ll; while (__sort_lt(*ll, *low));                                  \
 249                                 do --hh; while (__sort_lt(*low, *hh));                                  \
 250                                 if (hh < ll) break;                                                                             \
 251                                 KSORT_SWAP(type_t, *ll, *hh);                                                   \
 252                         }                                                                                                                       \
 253                         KSORT_SWAP(type_t, *low, *hh);                                                          \
 254                         if (hh <= k) low = ll;                                                                          \
 255                         if (hh >= k) high = hh - 1;                                                                     \
 256                 }                                                                                                                               \
 257         }                                                                                                                                       \
 258         void ks_shuffle_##name(size_t n, type_t a[])                                            \
 259         {                                                                                                                                       \
 260                 int i, j;                                                                                                               \
 261                 for (i = n; i > 1; --i) {                                                                               \
 262                         type_t tmp;                                                                                                     \
 263                         j = (int)(drand48() * i);                                                                       \
 264                         tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp;                                        \
 265                 }                                                                                                                               \
 266         }
 267
 268 #define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t)
 269 #define ks_introsort(name, n, a) ks_introsort_##name(n, a)
 270 #define ks_combsort(name, n, a) ks_combsort_##name(n, a)
 271 #define ks_heapsort(name, n, a) ks_heapsort_##name(n, a)
 272 #define ks_heapmake(name, n, a) ks_heapmake_##name(n, a)
 273 #define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a)
 274 #define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k)
 275 #define ks_shuffle(name, n, a) ks_shuffle_##name(n, a)
 276
 277 #define ks_lt_generic(a, b) ((a) < (b))
 278 #define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
 279
 280 typedef const char *ksstr_t;
 281
 282 #define KSORT_INIT_GENERIC(type_t) KSORT_INIT(type_t, type_t, ks_lt_generic)
 283 #define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str)
 284
 285 #endif