3 Copyright (c) 2008 Genome Research Ltd (GRL).
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:
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
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
26 /* Contact: Heng Li <lh3@sanger.ac.uk> */
31 * Fixed a bug in introsort() that happens in rare cases.
35 * Fixed a bug in introsort() for complex comparisons.
37 * Fixed a bug in mergesort(). The previous version is not stable.
41 * Accelerated introsort. On my Mac (not on another Linux machine),
42 my implementation is as fast as std::sort on random input.
44 * Added combsort and in introsort, switch to combsort if the
45 recursion is too deep.
49 * Added k-small algorithm
68 #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
70 #define KSORT_INIT(name, type_t, __sort_lt) \
71 void ks_mergesort_##name(size_t n, type_t array[], type_t temp[]) \
73 type_t *a2[2], *a, *b; \
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]; \
81 type_t *p = b, *i, *eb = a + n; \
82 for (i = a; i < eb; i += 2) { \
83 if (i == eb - 1) *p++ = *i; \
85 if (__sort_lt(*(i+1), *i)) { \
86 *p++ = *(i+1); *p++ = *i; \
88 *p++ = *i; *p++ = *(i+1); \
93 size_t i, step = 1ul<<shift; \
94 for (i = 0; i < n; i += step<<1) { \
95 type_t *p, *j, *k, *ea, *eb; \
100 eb = a + (n < i + (step<<1)? n : i + (step<<1)); \
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++; \
107 while (j < ea) *p++ = *j++; \
108 while (k < eb) *p++ = *k++; \
114 type_t *p = a2[0], *i = a2[1], *eb = array + n; \
115 for (; p < eb; ++i) *p++ = *i; \
117 if (temp == 0) free(a2[1]); \
119 void ks_heapadjust_##name(size_t i, size_t n, type_t l[]) \
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; \
130 void ks_heapmake_##name(size_t lsize, type_t l[]) \
133 for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \
134 ks_heapadjust_##name(i, lsize, l); \
136 void ks_heapsort_##name(size_t lsize, type_t l[]) \
139 for (i = lsize - 1; i > 0; --i) { \
141 tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust_##name(0, i, l); \
144 inline void __ks_insertsort_##name(type_t *s, type_t *t) \
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; \
152 void ks_combsort_##name(size_t n, type_t a[]) \
154 const double shrink_factor = 1.2473309501039786540366528676643; \
157 type_t tmp, *i, *j; \
160 gap = (size_t)(gap / shrink_factor); \
161 if (gap == 9 || gap == 10) gap = 11; \
164 for (i = a; i < a + n - gap; ++i) { \
166 if (__sort_lt(*j, *i)) { \
167 tmp = *i; *i = *j; *j = tmp; \
171 } while (do_swap || gap > 2); \
172 if (gap != 1) __ks_insertsort_##name(a, a + n); \
174 void ks_introsort_##name(size_t n, type_t a[]) \
177 ks_isort_stack_t *top, *stack; \
178 type_t rp, swap_tmp; \
179 type_t *s, *t, *i, *j, *k; \
183 if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \
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; \
192 ks_combsort_##name(t - s + 1, s); \
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; \
201 if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; } \
203 do ++i; while (__sort_lt(*i, rp)); \
204 do --j; while (i <= j && __sort_lt(rp, *j)); \
206 swap_tmp = *i; *i = *j; *j = swap_tmp; \
208 swap_tmp = *i; *i = *t; *t = swap_tmp; \
210 if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \
211 s = t-i > 16? i+1 : t; \
213 if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \
214 t = i-s > 16? i-1 : s; \
217 if (top == stack) { \
219 __ks_insertsort_##name(a, a+n); \
221 } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \
225 /* This function is adapted from: http://ndevilla.free.fr/median/ */ \
227 type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk) \
229 type_t *low, *high, *k, *ll, *hh, *mid; \
230 low = arr; high = arr + n - 1; k = arr + kk; \
232 if (high <= low) return *k; \
233 if (high == low + 1) { \
234 if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \
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; \
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); \
249 KSORT_SWAP(type_t, *low, *hh); \
250 if (hh <= k) low = ll; \
251 if (hh >= k) high = hh - 1; \
254 void ks_shuffle_##name(size_t n, type_t a[]) \
257 for (i = n; i > 1; --i) { \
259 j = (int)(drand48() * i); \
260 tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; \
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)
273 #define ks_lt_generic(a, b) ((a) < (b))
274 #define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
276 typedef const char *ksstr_t;
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)