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   2008-11-16 (0.1.4):
  30
  31     * Fixed a bug in introsort() that happens in rare cases.
  32
  33   2008-11-05 (0.1.3):
  34
  35     * Fixed a bug in introsort() for complex comparisons.
  36
  37         * Fixed a bug in mergesort(). The previous version is not stable.
  38
  39   2008-09-15 (0.1.2):
  40
  41         * Accelerated introsort. On my Mac (not on another Linux machine),
  42           my implementation is as fast as std::sort on random input.
  43
  44         * Added combsort and in introsort, switch to combsort if the
  45           recursion is too deep.
  46
  47   2008-09-13 (0.1.1):
  48
  49         * Added k-small algorithm
  50
  51   2008-09-05 (0.1.0):
  52
  53         * Initial version
  54
  55 */
  56
  57 #ifndef AC_KSORT_H
  58 #define AC_KSORT_H
  59
  60 #include <stdlib.h>
  61 #include <string.h>
  62
  63 typedef struct {
  64         void *left, *right;
  65         int depth;
  66 } ks_isort_stack_t;
  67
  68 #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
  69
  70 #define KSORT_INIT(name, type_t, __sort_lt)                                                             \
  71         void ks_mergesort_##name(size_t n, type_t array[], type_t temp[])       \
  72         {                                                                                                                                       \
  73                 type_t *a2[2], *a, *b;                                                                                  \
  74                 int curr, shift;                                                                                                \
  75                                                                                                                                                 \
  76                 a2[0] = array;                                                                                                  \
  77                 a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n);               \
  78                 for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) {                  \
  79                         a = a2[curr]; b = a2[1-curr];                                                           \
  80                         if (shift == 0) {                                                                                       \
  81                                 type_t *p = b, *i, *eb = a + n;                                                 \
  82                                 for (i = a; i < eb; i += 2) {                                                   \
  83                                         if (i == eb - 1) *p++ = *i;                                                     \
  84                                         else {                                                                                          \
  85                                                 if (__sort_lt(*(i+1), *i)) {                                    \
  86                                                         *p++ = *(i+1); *p++ = *i;                                       \
  87                                                 } else {                                                                                \
  88                                                         *p++ = *i; *p++ = *(i+1);                                       \
  89                                                 }                                                                                               \
  90                                         }                                                                                                       \
  91                                 }                                                                                                               \
  92                         } else {                                                                                                        \
  93                                 size_t i, step = 1ul<<shift;                                                    \
  94                                 for (i = 0; i < n; i += step<<1) {                                              \
  95                                         type_t *p, *j, *k, *ea, *eb;                                            \
  96                                         if (n < i + step) {                                                                     \
  97                                                 ea = a + n; eb = a;                                                             \
  98                                         } else {                                                                                        \
  99                                                 ea = a + i + step;                                                              \
 100                                                 eb = a + (n < i + (step<<1)? n : i + (step<<1)); \
 101                                         }                                                                                                       \
 102                                         j = a + i; k = a + i + step; p = b + i;                         \
 103                                         while (j < ea && k < eb) {                                                      \
 104                                                 if (__sort_lt(*k, *j)) *p++ = *k++;                             \
 105                                                 else *p++ = *j++;                                                               \
 106                                         }                                                                                                       \
 107                                         while (j < ea) *p++ = *j++;                                                     \
 108                                         while (k < eb) *p++ = *k++;                                                     \
 109                                 }                                                                                                               \
 110                         }                                                                                                                       \
 111                         curr = 1 - curr;                                                                                        \
 112                 }                                                                                                                               \
 113                 if (curr == 1) {                                                                                                \
 114                         type_t *p = a2[0], *i = a2[1], *eb = array + n;                         \
 115                         for (; p < eb; ++i) *p++ = *i;                                                          \
 116                 }                                                                                                                               \
 117                 if (temp == 0) free(a2[1]);                                                                             \
 118         }                                                                                                                                       \
 119         void ks_heapadjust_##name(size_t i, size_t n, type_t l[])                       \
 120         {                                                                                                                                       \
 121                 size_t k = i;                                                                                                   \
 122                 type_t tmp = l[i];                                                                                              \
 123                 while ((k = (k << 1) + 1) < n) {                                                                \
 124                         if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k;                         \
 125                         if (__sort_lt(l[k], tmp)) break;                                                        \
 126                         l[i] = l[k]; i = k;                                                                                     \
 127                 }                                                                                                                               \
 128                 l[i] = tmp;                                                                                                             \
 129         }                                                                                                                                       \
 130         void ks_heapmake_##name(size_t lsize, type_t l[])                                       \
 131         {                                                                                                                                       \
 132                 size_t i;                                                                                                               \
 133                 for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i)                              \
 134                         ks_heapadjust_##name(i, lsize, l);                                                      \
 135         }                                                                                                                                       \
 136         void ks_heapsort_##name(size_t lsize, type_t l[])                                       \
 137         {                                                                                                                                       \
 138                 size_t i;                                                                                                               \
 139                 for (i = lsize - 1; i > 0; --i) {                                                               \
 140                         type_t tmp;                                                                                                     \
 141                         tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
 142                 }                                                                                                                               \
 143         }                                                                                                                                       \
 144         inline void __ks_insertsort_##name(type_t *s, type_t *t)                        \
 145         {                                                                                                                                       \
 146                 type_t *i, *j, swap_tmp;                                                                                \
 147                 for (i = s + 1; i < t; ++i)                                                                             \
 148                         for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) {                      \
 149                                 swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp;                  \
 150                         }                                                                                                                       \
 151         }                                                                                                                                       \
 152         void ks_combsort_##name(size_t n, type_t a[])                                           \
 153         {                                                                                                                                       \
 154                 const double shrink_factor = 1.2473309501039786540366528676643; \
 155                 int do_swap;                                                                                                    \
 156                 size_t gap = n;                                                                                                 \
 157                 type_t tmp, *i, *j;                                                                                             \
 158                 do {                                                                                                                    \
 159                         if (gap > 2) {                                                                                          \
 160                                 gap = (size_t)(gap / shrink_factor);                                    \
 161                                 if (gap == 9 || gap == 10) gap = 11;                                    \
 162                         }                                                                                                                       \
 163                         do_swap = 0;                                                                                            \
 164                         for (i = a; i < a + n - gap; ++i) {                                                     \
 165                                 j = i + gap;                                                                                    \
 166                                 if (__sort_lt(*j, *i)) {                                                                \
 167                                         tmp = *i; *i = *j; *j = tmp;                                            \
 168                                         do_swap = 1;                                                                            \
 169                                 }                                                                                                               \
 170                         }                                                                                                                       \
 171                 } while (do_swap || gap > 2);                                                                   \
 172                 if (gap != 1) __ks_insertsort_##name(a, a + n);                                 \
 173         }                                                                                                                                       \
 174         void ks_introsort_##name(size_t n, type_t a[])                                          \
 175         {                                                                                                                                       \
 176                 int d;                                                                                                                  \
 177                 ks_isort_stack_t *top, *stack;                                                                  \
 178                 type_t rp, swap_tmp;                                                                                    \
 179                 type_t *s, *t, *i, *j, *k;                                                                              \
 180                                                                                                                                                 \
 181                 if (n < 1) return;                                                                                              \
 182                 else if (n == 2) {                                                                                              \
 183                         if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
 184                         return;                                                                                                         \
 185                 }                                                                                                                               \
 186                 for (d = 2; 1ul<<d < n; ++d);                                                                   \
 187                 stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
 188                 top = stack; s = a; t = a + (n-1); d <<= 1;                                             \
 189                 while (1) {                                                                                                             \
 190                         if (s < t) {                                                                                            \
 191                                 if (--d == 0) {                                                                                 \
 192                                         ks_combsort_##name(t - s + 1, s);                                       \
 193                                         t = s;                                                                                          \
 194                                         continue;                                                                                       \
 195                                 }                                                                                                               \
 196                                 i = s; j = t; k = i + ((j-i)>>1) + 1;                                   \
 197                                 if (__sort_lt(*k, *i)) {                                                                \
 198                                         if (__sort_lt(*k, *j)) k = j;                                           \
 199                                 } else k = __sort_lt(*j, *i)? i : j;                                    \
 200                                 rp = *k;                                                                                                \
 201                                 if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; }  \
 202                                 for (;;) {                                                                                              \
 203                                         do ++i; while (__sort_lt(*i, rp));                                      \
 204                                         do --j; while (i <= j && __sort_lt(rp, *j));            \
 205                                         if (j <= i) break;                                                                      \
 206                                         swap_tmp = *i; *i = *j; *j = swap_tmp;                          \
 207                                 }                                                                                                               \
 208                                 swap_tmp = *i; *i = *t; *t = swap_tmp;                                  \
 209                                 if (i-s > t-i) {                                                                                \
 210                                         if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
 211                                         s = t-i > 16? i+1 : t;                                                          \
 212                                 } else {                                                                                                \
 213                                         if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
 214                                         t = i-s > 16? i-1 : s;                                                          \
 215                                 }                                                                                                               \
 216                         } else {                                                                                                        \
 217                                 if (top == stack) {                                                                             \
 218                                         free(stack);                                                                            \
 219                                         __ks_insertsort_##name(a, a+n);                                         \
 220                                         return;                                                                                         \
 221                                 } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
 222                         }                                                                                                                       \
 223                 }                                                                                                                               \
 224         }                                                                                                                                       \
 225         /* This function is adapted from: http://ndevilla.free.fr/median/ */ \
 226         /* 0 <= kk < n */                                                                                                       \
 227         type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk)                      \
 228         {                                                                                                                                       \
 229                 type_t *low, *high, *k, *ll, *hh, *mid;                                                 \
 230                 low = arr; high = arr + n - 1; k = arr + kk;                                    \
 231                 for (;;) {                                                                                                              \
 232                         if (high <= low) return *k;                                                                     \
 233                         if (high == low + 1) {                                                                          \
 234                                 if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
 235                                 return *k;                                                                                              \
 236                         }                                                                                                                       \
 237                         mid = low + (high - low) / 2;                                                           \
 238                         if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \
 239                         if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
 240                         if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low);      \
 241                         KSORT_SWAP(type_t, *mid, *(low+1));                                                     \
 242                         ll = low + 1; hh = high;                                                                        \
 243                         for (;;) {                                                                                                      \
 244                                 do ++ll; while (__sort_lt(*ll, *low));                                  \
 245                                 do --hh; while (__sort_lt(*low, *hh));                                  \
 246                                 if (hh < ll) break;                                                                             \
 247                                 KSORT_SWAP(type_t, *ll, *hh);                                                   \
 248                         }                                                                                                                       \
 249                         KSORT_SWAP(type_t, *low, *hh);                                                          \
 250                         if (hh <= k) low = ll;                                                                          \
 251                         if (hh >= k) high = hh - 1;                                                                     \
 252                 }                                                                                                                               \
 253         }                                                                                                                                       \
 254         void ks_shuffle_##name(size_t n, type_t a[])                                            \
 255         {                                                                                                                                       \
 256                 int i, j;                                                                                                               \
 257                 for (i = n; i > 1; --i) {                                                                               \
 258                         type_t tmp;                                                                                                     \
 259                         j = (int)(drand48() * i);                                                                       \
 260                         tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp;                                        \
 261                 }                                                                                                                               \
 262         }
 263
 264 #define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t)
 265 #define ks_introsort(name, n, a) ks_introsort_##name(n, a)
 266 #define ks_combsort(name, n, a) ks_combsort_##name(n, a)
 267 #define ks_heapsort(name, n, a) ks_heapsort_##name(n, a)
 268 #define ks_heapmake(name, n, a) ks_heapmake_##name(n, a)
 269 #define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a)
 270 #define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k)
 271 #define ks_shuffle(name, n, a) ks_shuffle_##name(n, a)
 272
 273 #define ks_lt_generic(a, b) ((a) < (b))
 274 #define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
 275
 276 typedef const char *ksstr_t;
 277
 278 #define KSORT_INIT_GENERIC(type_t) KSORT_INIT(type_t, type_t, ks_lt_generic)
 279 #define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str)
 280
 281 #endif