/src/htslib/htslib/ksort.h
Line | Count | Source |
1 | | /* The MIT License |
2 | | |
3 | | Copyright (c) 2008, 2012-2013, 2017-2019 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 the C++ standard library's sort() |
47 | | on random input. |
48 | | |
49 | | * Added combsort and in introsort, switch to combsort if the |
50 | | recursion is too deep. |
51 | | |
52 | | 2008-09-13 (0.1.1): |
53 | | |
54 | | * Added k-small algorithm |
55 | | |
56 | | 2008-09-05 (0.1.0): |
57 | | |
58 | | * Initial version |
59 | | |
60 | | */ |
61 | | |
62 | | #ifndef AC_KSORT_H |
63 | | #define AC_KSORT_H |
64 | | |
65 | | #include <stdlib.h> |
66 | | #include <string.h> |
67 | | #include "hts_defs.h" |
68 | | |
69 | | #ifndef klib_unused |
70 | | #if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) |
71 | | #define klib_unused __attribute__ ((__unused__)) |
72 | | #else |
73 | | #define klib_unused |
74 | | #endif |
75 | | #endif /* klib_unused */ |
76 | | |
77 | | #ifdef __cplusplus |
78 | | extern "C" { |
79 | | #endif |
80 | | |
81 | | // Use our own drand48() symbol (used by ks_shuffle) to avoid portability |
82 | | // problems on Windows. Don't include htslib/hts_os.h for this as it |
83 | | // may not get on with older attempts to fix this in code that includes |
84 | | // this file. |
85 | | HTSLIB_EXPORT |
86 | | extern double hts_drand48(void); |
87 | | |
88 | | typedef struct { |
89 | | void *left, *right; |
90 | | int depth; |
91 | | } ks_isort_stack_t; |
92 | | |
93 | | #define KSORT_SWAP(type_t, a, b) { type_t t=(a); (a)=(b); (b)=t; } |
94 | | |
95 | | #define KSORT_INIT(name, type_t, __sort_lt) KSORT_INIT_(_ ## name, , type_t, __sort_lt) |
96 | | #define KSORT_INIT_STATIC(name, type_t, __sort_lt) KSORT_INIT_(_ ## name, static klib_unused, type_t, __sort_lt) |
97 | | #define KSORT_INIT2(name, SCOPE, type_t, __sort_lt) KSORT_INIT_(_ ## name, SCOPE, type_t, __sort_lt) |
98 | | |
99 | | #define KSORT_INIT_(name, SCOPE, type_t, __sort_lt) \ |
100 | | SCOPE int ks_mergesort##name(size_t n, type_t array[], type_t temp[]) \ |
101 | 0 | { \ |
102 | 0 | type_t *a2[2], *a, *b; \ |
103 | 0 | int curr, shift; \ |
104 | 0 | \ |
105 | 0 | a2[0] = array; \ |
106 | 0 | a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n); \ |
107 | 0 | for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \ |
108 | 0 | a = a2[curr]; b = a2[1-curr]; \ |
109 | 0 | if (shift == 0) { \ |
110 | 0 | type_t *p = b, *i, *eb = a + n; \ |
111 | 0 | for (i = a; i < eb; i += 2) { \ |
112 | 0 | if (i == eb - 1) *p++ = *i; \ |
113 | 0 | else { \ |
114 | 0 | if (__sort_lt(*(i+1), *i)) { \ |
115 | 0 | *p++ = *(i+1); *p++ = *i; \ |
116 | 0 | } else { \ |
117 | 0 | *p++ = *i; *p++ = *(i+1); \ |
118 | 0 | } \ |
119 | 0 | } \ |
120 | 0 | } \ |
121 | 0 | } else { \ |
122 | 0 | size_t i, step = 1ul<<shift; \ |
123 | 0 | for (i = 0; i < n; i += step<<1) { \ |
124 | 0 | type_t *p, *j, *k, *ea, *eb; \ |
125 | 0 | if (n < i + step) { \ |
126 | 0 | ea = a + n; eb = a; \ |
127 | 0 | } else { \ |
128 | 0 | ea = a + i + step; \ |
129 | 0 | eb = a + (n < i + (step<<1)? n : i + (step<<1)); \ |
130 | 0 | } \ |
131 | 0 | j = a + i; k = a + i + step; p = b + i; \ |
132 | 0 | while (j < ea && k < eb) { \ |
133 | 0 | if (__sort_lt(*k, *j)) *p++ = *k++; \ |
134 | 0 | else *p++ = *j++; \ |
135 | 0 | } \ |
136 | 0 | while (j < ea) *p++ = *j++; \ |
137 | 0 | while (k < eb) *p++ = *k++; \ |
138 | 0 | } \ |
139 | 0 | } \ |
140 | 0 | curr = 1 - curr; \ |
141 | 0 | } \ |
142 | 0 | if (curr == 1) { \ |
143 | 0 | type_t *p = a2[0], *i = a2[1], *eb = array + n; \ |
144 | 0 | for (; p < eb; ++i) *p++ = *i; \ |
145 | 0 | } \ |
146 | 0 | if (temp == 0) free(a2[1]); \ |
147 | 0 | return 0; \ |
148 | 0 | } \ Unexecuted instantiation: hts.c:ks_mergesort__off Unexecuted instantiation: hts.c:ks_mergesort__off_max |
149 | | SCOPE void ks_heapadjust##name(size_t i, size_t n, type_t l[]) \ |
150 | 0 | { \ |
151 | 0 | size_t k = i; \ |
152 | 0 | type_t tmp = l[i]; \ |
153 | 0 | while ((k = (k << 1) + 1) < n) { \ |
154 | 0 | if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \ |
155 | 0 | if (__sort_lt(l[k], tmp)) break; \ |
156 | 0 | l[i] = l[k]; i = k; \ |
157 | 0 | } \ |
158 | 0 | l[i] = tmp; \ |
159 | 0 | } \ Unexecuted instantiation: hts.c:ks_heapadjust__off Unexecuted instantiation: hts.c:ks_heapadjust__off_max |
160 | | SCOPE void ks_heapmake##name(size_t lsize, type_t l[]) \ |
161 | 0 | { \ |
162 | 0 | size_t i; \ |
163 | 0 | for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \ |
164 | 0 | ks_heapadjust##name(i, lsize, l); \ |
165 | 0 | } \ Unexecuted instantiation: hts.c:ks_heapmake__off Unexecuted instantiation: hts.c:ks_heapmake__off_max |
166 | | SCOPE void ks_heapsort##name(size_t lsize, type_t l[]) \ |
167 | 0 | { \ |
168 | 0 | size_t i; \ |
169 | 0 | for (i = lsize - 1; i > 0; --i) { \ |
170 | 0 | type_t tmp; \ |
171 | 0 | tmp = *l; *l = l[i]; l[i] = tmp; ks_heapadjust##name(0, i, l); \ |
172 | 0 | } \ |
173 | 0 | } \ Unexecuted instantiation: hts.c:ks_heapsort__off Unexecuted instantiation: hts.c:ks_heapsort__off_max |
174 | | static inline void __ks_insertsort##name(type_t *s, type_t *t) \ |
175 | 0 | { \ |
176 | 0 | type_t *i, *j, swap_tmp; \ |
177 | 0 | for (i = s + 1; i < t; ++i) \ |
178 | 0 | for (j = i; j > s && __sort_lt(*j, *(j-1)); --j) { \ |
179 | 0 | swap_tmp = *j; *j = *(j-1); *(j-1) = swap_tmp; \ |
180 | 0 | } \ |
181 | 0 | } \ Unexecuted instantiation: hts.c:__ks_insertsort__off Unexecuted instantiation: hts.c:__ks_insertsort__off_max |
182 | | SCOPE void ks_combsort##name(size_t n, type_t a[]) \ |
183 | 0 | { \ |
184 | 0 | const double shrink_factor = 1.2473309501039786540366528676643; \ |
185 | 0 | int do_swap; \ |
186 | 0 | size_t gap = n; \ |
187 | 0 | type_t tmp, *i, *j; \ |
188 | 0 | do { \ |
189 | 0 | if (gap > 2) { \ |
190 | 0 | gap = (size_t)(gap / shrink_factor); \ |
191 | 0 | if (gap == 9 || gap == 10) gap = 11; \ |
192 | 0 | } \ |
193 | 0 | do_swap = 0; \ |
194 | 0 | for (i = a; i < a + n - gap; ++i) { \ |
195 | 0 | j = i + gap; \ |
196 | 0 | if (__sort_lt(*j, *i)) { \ |
197 | 0 | tmp = *i; *i = *j; *j = tmp; \ |
198 | 0 | do_swap = 1; \ |
199 | 0 | } \ |
200 | 0 | } \ |
201 | 0 | } while (do_swap || gap > 2); \ |
202 | 0 | if (gap != 1) __ks_insertsort##name(a, a + n); \ |
203 | 0 | } \ Unexecuted instantiation: hts.c:ks_combsort__off Unexecuted instantiation: hts.c:ks_combsort__off_max |
204 | | SCOPE int ks_introsort##name(size_t n, type_t a[]) \ |
205 | 0 | { \ |
206 | 0 | int d; \ |
207 | 0 | ks_isort_stack_t *top, *stack; \ |
208 | 0 | type_t rp, swap_tmp; \ |
209 | 0 | type_t *s, *t, *i, *j, *k; \ |
210 | 0 | \ |
211 | 0 | if (n < 1) return 0; \ |
212 | 0 | else if (n == 2) { \ |
213 | 0 | if (__sort_lt(a[1], a[0])) { swap_tmp = a[0]; a[0] = a[1]; a[1] = swap_tmp; } \ |
214 | 0 | return 0; \ |
215 | 0 | } \ |
216 | 0 | for (d = 2; 1ul<<d < n; ++d); \ |
217 | 0 | stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \ |
218 | 0 | top = stack; s = a; t = a + (n-1); d <<= 1; \ |
219 | 0 | while (1) { \ |
220 | 0 | if (s < t) { \ |
221 | 0 | if (--d == 0) { \ |
222 | 0 | ks_combsort##name(t - s + 1, s); \ |
223 | 0 | t = s; \ |
224 | 0 | continue; \ |
225 | 0 | } \ |
226 | 0 | i = s; j = t; k = i + ((j-i)>>1) + 1; \ |
227 | 0 | if (__sort_lt(*k, *i)) { \ |
228 | 0 | if (__sort_lt(*k, *j)) k = j; \ |
229 | 0 | } else k = __sort_lt(*j, *i)? i : j; \ |
230 | 0 | rp = *k; \ |
231 | 0 | if (k != t) { swap_tmp = *k; *k = *t; *t = swap_tmp; } \ |
232 | 0 | for (;;) { \ |
233 | 0 | do ++i; while (__sort_lt(*i, rp)); \ |
234 | 0 | do --j; while (i <= j && __sort_lt(rp, *j)); \ |
235 | 0 | if (j <= i) break; \ |
236 | 0 | swap_tmp = *i; *i = *j; *j = swap_tmp; \ |
237 | 0 | } \ |
238 | 0 | swap_tmp = *i; *i = *t; *t = swap_tmp; \ |
239 | 0 | if (i-s > t-i) { \ |
240 | 0 | if (i-s > 16) { top->left = s; top->right = i-1; top->depth = d; ++top; } \ |
241 | 0 | s = t-i > 16? i+1 : t; \ |
242 | 0 | } else { \ |
243 | 0 | if (t-i > 16) { top->left = i+1; top->right = t; top->depth = d; ++top; } \ |
244 | 0 | t = i-s > 16? i-1 : s; \ |
245 | 0 | } \ |
246 | 0 | } else { \ |
247 | 0 | if (top == stack) { \ |
248 | 0 | free(stack); \ |
249 | 0 | __ks_insertsort##name(a, a+n); \ |
250 | 0 | return 0; \ |
251 | 0 | } else { --top; s = (type_t*)top->left; t = (type_t*)top->right; d = top->depth; } \ |
252 | 0 | } \ |
253 | 0 | } \ |
254 | 0 | return 0; \ |
255 | 0 | } \ Unexecuted instantiation: hts.c:ks_introsort__off Unexecuted instantiation: hts.c:ks_introsort__off_max |
256 | | /* This function is adapted from: http://ndevilla.free.fr/median/ */ \ |
257 | | /* 0 <= kk < n */ \ |
258 | | SCOPE type_t ks_ksmall##name(size_t n, type_t arr[], size_t kk) \ |
259 | 0 | { \ |
260 | 0 | type_t *low, *high, *k, *ll, *hh, *mid; \ |
261 | 0 | low = arr; high = arr + n - 1; k = arr + kk; \ |
262 | 0 | for (;;) { \ |
263 | 0 | if (high <= low) return *k; \ |
264 | 0 | if (high == low + 1) { \ |
265 | 0 | if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \ |
266 | 0 | return *k; \ |
267 | 0 | } \ |
268 | 0 | mid = low + (high - low) / 2; \ |
269 | 0 | if (__sort_lt(*high, *mid)) KSORT_SWAP(type_t, *mid, *high); \ |
270 | 0 | if (__sort_lt(*high, *low)) KSORT_SWAP(type_t, *low, *high); \ |
271 | 0 | if (__sort_lt(*low, *mid)) KSORT_SWAP(type_t, *mid, *low); \ |
272 | 0 | KSORT_SWAP(type_t, *mid, *(low+1)); \ |
273 | 0 | ll = low + 1; hh = high; \ |
274 | 0 | for (;;) { \ |
275 | 0 | do ++ll; while (__sort_lt(*ll, *low)); \ |
276 | 0 | do --hh; while (__sort_lt(*low, *hh)); \ |
277 | 0 | if (hh < ll) break; \ |
278 | 0 | KSORT_SWAP(type_t, *ll, *hh); \ |
279 | 0 | } \ |
280 | 0 | KSORT_SWAP(type_t, *low, *hh); \ |
281 | 0 | if (hh <= k) low = ll; \ |
282 | 0 | if (hh >= k) high = hh - 1; \ |
283 | 0 | } \ |
284 | 0 | } \ Unexecuted instantiation: hts.c:ks_ksmall__off Unexecuted instantiation: hts.c:ks_ksmall__off_max |
285 | | SCOPE void ks_shuffle##name(size_t n, type_t a[]) \ |
286 | 0 | { \ |
287 | 0 | int i, j; \ |
288 | 0 | for (i = n; i > 1; --i) { \ |
289 | 0 | type_t tmp; \ |
290 | 0 | j = (int)(hts_drand48() * i); \ |
291 | 0 | tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; \ |
292 | 0 | } \ |
293 | 0 | } Unexecuted instantiation: hts.c:ks_shuffle__off Unexecuted instantiation: hts.c:ks_shuffle__off_max |
294 | | |
295 | | #define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t) |
296 | 0 | #define ks_introsort(name, n, a) ks_introsort_##name(n, a) |
297 | | #define ks_combsort(name, n, a) ks_combsort_##name(n, a) |
298 | | #define ks_heapsort(name, n, a) ks_heapsort_##name(n, a) |
299 | | #define ks_heapmake(name, n, a) ks_heapmake_##name(n, a) |
300 | | #define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a) |
301 | | #define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k) |
302 | | #define ks_shuffle(name, n, a) ks_shuffle_##name(n, a) |
303 | | |
304 | | #define ks_lt_generic(a, b) ((a) < (b)) |
305 | | #define ks_lt_str(a, b) (strcmp((a), (b)) < 0) |
306 | | |
307 | | typedef const char *ksstr_t; |
308 | | |
309 | | #define KSORT_INIT_GENERIC(type_t) KSORT_INIT_(_ ## type_t, , type_t, ks_lt_generic) |
310 | | #define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str) |
311 | | |
312 | | #define KSORT_INIT_STATIC_GENERIC(type_t) KSORT_INIT_(_ ## type_t, static klib_unused, type_t, ks_lt_generic) |
313 | | #define KSORT_INIT_STATIC_STR KSORT_INIT_STATIC(str, ksstr_t, ks_lt_str) |
314 | | |
315 | | #define KSORT_INIT2_GENERIC(type_t, SCOPE) KSORT_INIT_(_ ## type_t, SCOPE, type_t, ks_lt_generic) |
316 | | #define KSORT_INIT2_STR KSORT_INIT2(str, SCOPE, ksstr_t, ks_lt_str) |
317 | | |
318 | | #ifdef __cplusplus |
319 | | } |
320 | | #endif |
321 | | |
322 | | #endif |