Michael Lawrence authored on 11/02/2019 20:33:07 | 1 | 1 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,295 @@ |
| 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 |
+ |
|
| 68 |
+/* Rsamtools */ |
|
| 69 |
+#if defined(_WIN32) |
|
| 70 |
+#define drand48() ((double)rand() / RAND_MAX) |
|
| 71 |
+#else |
|
| 72 |
+#ifndef _SVID_SOURCE |
|
| 73 |
+#define _SVID_SOURCE /* Rsamtools: c99 */ |
|
| 74 |
+#endif |
|
| 75 |
+#endif |
|
| 76 |
+ |
|
| 77 |
+typedef struct {
|
|
| 78 |
+ void *left, *right; |
|
| 79 |
+ int depth; |
|
| 80 |
+} ks_isort_stack_t; |
|
| 81 |
+ |
|
| 82 |
+#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
|
|
| 83 |
+ |
|
| 84 |
+#define KSORT_INIT(name, type_t, __sort_lt) \ |
|
| 85 |
+ void ks_mergesort_##name(size_t n, type_t array[], type_t temp[]) \ |
|
| 86 |
+ { \
|
|
| 87 |
+ type_t *a2[2], *a, *b; \ |
|
| 88 |
+ int curr, shift; \ |
|
| 89 |
+ \ |
|
| 90 |
+ a2[0] = array; \ |
|
| 91 |
+ a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n); \ |
|
| 92 |
+ for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \
|
|
| 93 |
+ a = a2[curr]; b = a2[1-curr]; \ |
|
| 94 |
+ if (shift == 0) { \
|
|
| 95 |
+ type_t *p = b, *i, *eb = a + n; \ |
|
| 96 |
+ for (i = a; i < eb; i += 2) { \
|
|
| 97 |
+ if (i == eb - 1) *p++ = *i; \ |
|
| 98 |
+ else { \
|
|
| 99 |
+ if (__sort_lt(*(i+1), *i)) { \
|
|
| 100 |
+ *p++ = *(i+1); *p++ = *i; \ |
|
| 101 |
+ } else { \
|
|
| 102 |
+ *p++ = *i; *p++ = *(i+1); \ |
|
| 103 |
+ } \ |
|
| 104 |
+ } \ |
|
| 105 |
+ } \ |
|
| 106 |
+ } else { \
|
|
| 107 |
+ size_t i, step = 1ul<<shift; \ |
|
| 108 |
+ for (i = 0; i < n; i += step<<1) { \
|
|
| 109 |
+ type_t *p, *j, *k, *ea, *eb; \ |
|
| 110 |
+ if (n < i + step) { \
|
|
| 111 |
+ ea = a + n; eb = a; \ |
|
| 112 |
+ } else { \
|
|
| 113 |
+ ea = a + i + step; \ |
|
| 114 |
+ eb = a + (n < i + (step<<1)? n : i + (step<<1)); \ |
|
| 115 |
+ } \ |
|
| 116 |
+ j = a + i; k = a + i + step; p = b + i; \ |
|
| 117 |
+ while (j < ea && k < eb) { \
|
|
| 118 |
+ if (__sort_lt(*k, *j)) *p++ = *k++; \ |
|
| 119 |
+ else *p++ = *j++; \ |
|
| 120 |
+ } \ |
|
| 121 |
+ while (j < ea) *p++ = *j++; \ |
|
| 122 |
+ while (k < eb) *p++ = *k++; \ |
|
| 123 |
+ } \ |
|
| 124 |
+ } \ |
|
| 125 |
+ curr = 1 - curr; \ |
|
| 126 |
+ } \ |
|
| 127 |
+ if (curr == 1) { \
|
|
| 128 |
+ type_t *p = a2[0], *i = a2[1], *eb = array + n; \ |
|
| 129 |
+ for (; p < eb; ++i) *p++ = *i; \ |
|
| 130 |
+ } \ |
|
| 131 |
+ if (temp == 0) free(a2[1]); \ |
|
| 132 |
+ } \ |
|
| 133 |
+ void ks_heapadjust_##name(size_t i, size_t n, type_t l[]) \ |
|
| 134 |
+ { \
|
|
| 135 |
+ size_t k = i; \ |
|
| 136 |
+ type_t tmp = l[i]; \ |
|
| 137 |
+ while ((k = (k << 1) + 1) < n) { \
|
|
| 138 |
+ if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \ |
|
| 139 |
+ if (__sort_lt(l[k], tmp)) break; \ |
|
| 140 |
+ l[i] = l[k]; i = k; \ |
|
| 141 |
+ } \ |
|
| 142 |
+ l[i] = tmp; \ |
|
| 143 |
+ } \ |
|
| 144 |
+ void ks_heapmake_##name(size_t lsize, type_t l[]) \ |
|
| 145 |
+ { \
|
|
| 146 |
+ size_t i; \ |
|
| 147 |
+ for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \ |
|
| 148 |
+ ks_heapadjust_##name(i, lsize, l); \ |
|
| 149 |
+ } \ |
|
| 150 |
+ void ks_heapsort_##name(size_t lsize, type_t l[]) \ |
|
| 151 |
+ { \
|
|
| 152 |
+ size_t i; \ |
|
| 153 |
+ for (i = lsize - 1; i > 0; --i) { \
|
|
| 154 |
+ type_t tmp; \ |
|
| 155 |
+ tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \ |
|
| 156 |
+ } \ |
|
| 157 |
+ } \ |
|
| 158 |
+ static inline void __ks_insertsort_##name(type_t *s, type_t *t) \ |
|
| 159 |
+ { \
|
|
| 160 |
+ type_t *i, *j, swap_tmp; \ |
|
| 161 |
+ for (i = s + 1; i < t; ++i) \ |
|
| 162 |
+ for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) { \
|
|
| 163 |
+ swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp; \ |
|
| 164 |
+ } \ |
|
| 165 |
+ } \ |
|
| 166 |
+ void ks_combsort_##name(size_t n, type_t a[]) \ |
|
| 167 |
+ { \
|
|
| 168 |
+ const double shrink_factor = 1.2473309501039786540366528676643; \ |
|
| 169 |
+ int do_swap; \ |
|
| 170 |
+ size_t gap = n; \ |
|
| 171 |
+ type_t tmp, *i, *j; \ |
|
| 172 |
+ do { \
|
|
| 173 |
+ if (gap > 2) { \
|
|
| 174 |
+ gap = (size_t)(gap / shrink_factor); \ |
|
| 175 |
+ if (gap == 9 || gap == 10) gap = 11; \ |
|
| 176 |
+ } \ |
|
| 177 |
+ do_swap = 0; \ |
|
| 178 |
+ for (i = a; i < a + n - gap; ++i) { \
|
|
| 179 |
+ j = i + gap; \ |
|
| 180 |
+ if (__sort_lt(*j, *i)) { \
|
|
| 181 |
+ tmp = *i; *i = *j; *j = tmp; \ |
|
| 182 |
+ do_swap = 1; \ |
|
| 183 |
+ } \ |
|
| 184 |
+ } \ |
|
| 185 |
+ } while (do_swap || gap > 2); \ |
|
| 186 |
+ if (gap != 1) __ks_insertsort_##name(a, a + n); \ |
|
| 187 |
+ } \ |
|
| 188 |
+ void ks_introsort_##name(size_t n, type_t a[]) \ |
|
| 189 |
+ { \
|
|
| 190 |
+ int d; \ |
|
| 191 |
+ ks_isort_stack_t *top, *stack; \ |
|
| 192 |
+ type_t rp, swap_tmp; \ |
|
| 193 |
+ type_t *s, *t, *i, *j, *k; \ |
|
| 194 |
+ \ |
|
| 195 |
+ if (n < 1) return; \ |
|
| 196 |
+ else if (n == 2) { \
|
|
| 197 |
+ if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
|
|
| 198 |
+ return; \ |
|
| 199 |
+ } \ |
|
| 200 |
+ for (d = 2; 1ul<<d < n; ++d); \ |
|
| 201 |
+ stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \ |
|
| 202 |
+ top = stack; s = a; t = a + (n-1); d <<= 1; \ |
|
| 203 |
+ while (1) { \
|
|
| 204 |
+ if (s < t) { \
|
|
| 205 |
+ if (--d == 0) { \
|
|
| 206 |
+ ks_combsort_##name(t - s + 1, s); \ |
|
| 207 |
+ t = s; \ |
|
| 208 |
+ continue; \ |
|
| 209 |
+ } \ |
|
| 210 |
+ i = s; j = t; k = i + ((j-i)>>1) + 1; \ |
|
| 211 |
+ if (__sort_lt(*k, *i)) { \
|
|
| 212 |
+ if (__sort_lt(*k, *j)) k = j; \ |
|
| 213 |
+ } else k = __sort_lt(*j, *i)? i : j; \ |
|
| 214 |
+ rp = *k; \ |
|
| 215 |
+ if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; } \
|
|
| 216 |
+ for (;;) { \
|
|
| 217 |
+ do ++i; while (__sort_lt(*i, rp)); \ |
|
| 218 |
+ do --j; while (i <= j && __sort_lt(rp, *j)); \ |
|
| 219 |
+ if (j <= i) break; \ |
|
| 220 |
+ swap_tmp = *i; *i = *j; *j = swap_tmp; \ |
|
| 221 |
+ } \ |
|
| 222 |
+ swap_tmp = *i; *i = *t; *t = swap_tmp; \ |
|
| 223 |
+ if (i-s > t-i) { \
|
|
| 224 |
+ if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
|
|
| 225 |
+ s = t-i > 16? i+1 : t; \ |
|
| 226 |
+ } else { \
|
|
| 227 |
+ if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
|
|
| 228 |
+ t = i-s > 16? i-1 : s; \ |
|
| 229 |
+ } \ |
|
| 230 |
+ } else { \
|
|
| 231 |
+ if (top == stack) { \
|
|
| 232 |
+ free(stack); \ |
|
| 233 |
+ __ks_insertsort_##name(a, a+n); \ |
|
| 234 |
+ return; \ |
|
| 235 |
+ } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
|
|
| 236 |
+ } \ |
|
| 237 |
+ } \ |
|
| 238 |
+ } \ |
|
| 239 |
+ /* This function is adapted from: http://ndevilla.free.fr/median/ */ \ |
|
| 240 |
+ /* 0 <= kk < n */ \ |
|
| 241 |
+ type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk) \ |
|
| 242 |
+ { \
|
|
| 243 |
+ type_t *low, *high, *k, *ll, *hh, *mid; \ |
|
| 244 |
+ low = arr; high = arr + n - 1; k = arr + kk; \ |
|
| 245 |
+ for (;;) { \
|
|
| 246 |
+ if (high <= low) return *k; \ |
|
| 247 |
+ if (high == low + 1) { \
|
|
| 248 |
+ if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \ |
|
| 249 |
+ return *k; \ |
|
| 250 |
+ } \ |
|
| 251 |
+ mid = low + (high - low) / 2; \ |
|
| 252 |
+ if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \ |
|
| 253 |
+ if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \ |
|
| 254 |
+ if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low); \ |
|
| 255 |
+ KSORT_SWAP(type_t, *mid, *(low+1)); \ |
|
| 256 |
+ ll = low + 1; hh = high; \ |
|
| 257 |
+ for (;;) { \
|
|
| 258 |
+ do ++ll; while (__sort_lt(*ll, *low)); \ |
|
| 259 |
+ do --hh; while (__sort_lt(*low, *hh)); \ |
|
| 260 |
+ if (hh < ll) break; \ |
|
| 261 |
+ KSORT_SWAP(type_t, *ll, *hh); \ |
|
| 262 |
+ } \ |
|
| 263 |
+ KSORT_SWAP(type_t, *low, *hh); \ |
|
| 264 |
+ if (hh <= k) low = ll; \ |
|
| 265 |
+ if (hh >= k) high = hh - 1; \ |
|
| 266 |
+ } \ |
|
| 267 |
+ } \ |
|
| 268 |
+ void ks_shuffle_##name(size_t n, type_t a[]) \ |
|
| 269 |
+ { \
|
|
| 270 |
+ int i, j; \ |
|
| 271 |
+ for (i = n; i > 1; --i) { \
|
|
| 272 |
+ type_t tmp; \ |
|
| 273 |
+ j = (int)(drand48() * i); \ |
|
| 274 |
+ tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; \ |
|
| 275 |
+ } \ |
|
| 276 |
+ } |
|
| 277 |
+ |
|
| 278 |
+#define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t) |
|
| 279 |
+#define ks_introsort(name, n, a) ks_introsort_##name(n, a) |
|
| 280 |
+#define ks_combsort(name, n, a) ks_combsort_##name(n, a) |
|
| 281 |
+#define ks_heapsort(name, n, a) ks_heapsort_##name(n, a) |
|
| 282 |
+#define ks_heapmake(name, n, a) ks_heapmake_##name(n, a) |
|
| 283 |
+#define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a) |
|
| 284 |
+#define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k) |
|
| 285 |
+#define ks_shuffle(name, n, a) ks_shuffle_##name(n, a) |
|
| 286 |
+ |
|
| 287 |
+#define ks_lt_generic(a, b) ((a) < (b)) |
|
| 288 |
+#define ks_lt_str(a, b) (strcmp((a), (b)) < 0) |
|
| 289 |
+ |
|
| 290 |
+typedef const char *ksstr_t; |
|
| 291 |
+ |
|
| 292 |
+#define KSORT_INIT_GENERIC(type_t) KSORT_INIT(type_t, type_t, ks_lt_generic) |
|
| 293 |
+#define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str) |
|
| 294 |
+ |
|
| 295 |
+#endif |