]>
Commit | Line | Data |
---|---|---|
6d7e8eda | 1 | /* Copyright (C) 1991-2023 Free Software Foundation, Inc. |
6d52618b | 2 | This file is part of the GNU C Library. |
28f540f4 | 3 | |
6d52618b | 4 | The GNU C Library is free software; you can redistribute it and/or |
41bdb6e2 AJ |
5 | modify it under the terms of the GNU Lesser General Public |
6 | License as published by the Free Software Foundation; either | |
7 | version 2.1 of the License, or (at your option) any later version. | |
28f540f4 | 8 | |
6d52618b UD |
9 | The GNU C Library is distributed in the hope that it will be useful, |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
41bdb6e2 | 12 | Lesser General Public License for more details. |
28f540f4 | 13 | |
41bdb6e2 | 14 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 | 15 | License along with the GNU C Library; if not, see |
5a82c748 | 16 | <https://www.gnu.org/licenses/>. */ |
28f540f4 | 17 | |
061d137b UD |
18 | /* If you consider tuning this algorithm, you should consult first: |
19 | Engineering a sort function; Jon Bentley and M. Douglas McIlroy; | |
20 | Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993. */ | |
21 | ||
22 | #include <alloca.h> | |
23 | #include <limits.h> | |
21d30c77 | 24 | #include <memswap.h> |
28f540f4 RM |
25 | #include <stdlib.h> |
26 | #include <string.h> | |
21d30c77 | 27 | #include <stdbool.h> |
28f540f4 | 28 | |
21d30c77 AZ |
29 | /* Swap SIZE bytes between addresses A and B. These helpers are provided |
30 | along the generic one as an optimization. */ | |
31 | ||
32 | enum swap_type_t | |
33 | { | |
34 | SWAP_WORDS_64, | |
35 | SWAP_WORDS_32, | |
36 | SWAP_BYTES | |
37 | }; | |
38 | ||
39 | /* If this function returns true, elements can be safely copied using word | |
40 | loads and stores. Otherwise, it might not be safe. BASE (as an integer) | |
41 | must be a multiple of the word alignment. SIZE must be a multiple of | |
42 | WORDSIZE. Since WORDSIZE must be a multiple of the word alignment, and | |
43 | WORDSIZE is a power of two on all supported platforms, this function for | |
44 | speed merely checks that BASE and SIZE are both multiples of the word | |
45 | size. */ | |
46 | static inline bool | |
47 | is_aligned (const void *base, size_t size, size_t wordsize) | |
48 | { | |
49 | return (((uintptr_t) base | size) & (wordsize - 1)) == 0; | |
50 | } | |
51 | ||
52 | static inline void | |
53 | swap_words_64 (void * restrict a, void * restrict b, size_t n) | |
54 | { | |
55 | typedef uint64_t __attribute__ ((__may_alias__)) u64_alias_t; | |
56 | do | |
57 | { | |
58 | n -= 8; | |
59 | u64_alias_t t = *(u64_alias_t *)(a + n); | |
60 | *(u64_alias_t *)(a + n) = *(u64_alias_t *)(b + n); | |
61 | *(u64_alias_t *)(b + n) = t; | |
62 | } while (n); | |
63 | } | |
64 | ||
65 | static inline void | |
66 | swap_words_32 (void * restrict a, void * restrict b, size_t n) | |
67 | { | |
68 | typedef uint32_t __attribute__ ((__may_alias__)) u32_alias_t; | |
69 | do | |
70 | { | |
71 | n -= 4; | |
72 | u32_alias_t t = *(u32_alias_t *)(a + n); | |
73 | *(u32_alias_t *)(a + n) = *(u32_alias_t *)(b + n); | |
74 | *(u32_alias_t *)(b + n) = t; | |
75 | } while (n); | |
76 | } | |
77 | ||
78 | /* Replace the indirect call with a serie of if statements. It should help | |
79 | the branch predictor. */ | |
80 | static void | |
81 | do_swap (void * restrict a, void * restrict b, size_t size, | |
82 | enum swap_type_t swap_type) | |
83 | { | |
84 | if (swap_type == SWAP_WORDS_64) | |
85 | swap_words_64 (a, b, size); | |
86 | else if (swap_type == SWAP_WORDS_32) | |
87 | swap_words_32 (a, b, size); | |
88 | else | |
89 | __memswap (a, b, size); | |
90 | } | |
28f540f4 RM |
91 | |
92 | /* Discontinue quicksort algorithm when partition gets below this size. | |
93 | This particular magic number was chosen to work best on a Sun 4/260. */ | |
94 | #define MAX_THRESH 4 | |
95 | ||
96 | /* Stack node declarations used to store unfulfilled partition obligations. */ | |
6d52618b | 97 | typedef struct |
28f540f4 RM |
98 | { |
99 | char *lo; | |
100 | char *hi; | |
274a46c9 | 101 | size_t depth; |
28f540f4 RM |
102 | } stack_node; |
103 | ||
061d137b UD |
104 | /* The stack needs log (total_elements) entries (we could even subtract |
105 | log(MAX_THRESH)). Since total_elements has type size_t, we get as | |
106 | upper bound for log (total_elements): | |
107 | bits per byte (CHAR_BIT) * sizeof(size_t). */ | |
d097f3c7 AZ |
108 | enum { STACK_SIZE = CHAR_BIT * sizeof (size_t) }; |
109 | ||
110 | static inline stack_node * | |
274a46c9 | 111 | push (stack_node *top, char *lo, char *hi, size_t depth) |
d097f3c7 AZ |
112 | { |
113 | top->lo = lo; | |
114 | top->hi = hi; | |
274a46c9 | 115 | top->depth = depth; |
d097f3c7 AZ |
116 | return ++top; |
117 | } | |
118 | ||
119 | static inline stack_node * | |
274a46c9 | 120 | pop (stack_node *top, char **lo, char **hi, size_t *depth) |
d097f3c7 AZ |
121 | { |
122 | --top; | |
123 | *lo = top->lo; | |
124 | *hi = top->hi; | |
274a46c9 | 125 | *depth = top->depth; |
d097f3c7 AZ |
126 | return top; |
127 | } | |
28f540f4 | 128 | |
274a46c9 AZ |
129 | /* NB: N is inclusive bound for BASE. */ |
130 | static inline void | |
131 | siftdown (void *base, size_t size, size_t k, size_t n, | |
132 | enum swap_type_t swap_type, __compar_d_fn_t cmp, void *arg) | |
133 | { | |
134 | while (k <= n / 2) | |
135 | { | |
136 | size_t j = 2 * k; | |
137 | if (j < n && cmp (base + (j * size), base + ((j + 1) * size), arg) < 0) | |
138 | j++; | |
139 | ||
140 | if (cmp (base + (k * size), base + (j * size), arg) >= 0) | |
141 | break; | |
142 | ||
143 | do_swap (base + (size * j), base + (k * size), size, swap_type); | |
144 | k = j; | |
145 | } | |
146 | } | |
147 | ||
148 | static inline void | |
149 | heapify (void *base, size_t size, size_t n, enum swap_type_t swap_type, | |
150 | __compar_d_fn_t cmp, void *arg) | |
151 | { | |
152 | size_t k = n / 2; | |
153 | while (1) | |
154 | { | |
155 | siftdown (base, size, k, n, swap_type, cmp, arg); | |
156 | if (k-- == 0) | |
157 | break; | |
158 | } | |
159 | } | |
160 | ||
161 | /* A non-recursive heapsort, used on introsort implementation as a fallback | |
162 | routine with worst-case performance of O(nlog n) and worst-case space | |
163 | complexity of O(1). It sorts the array starting at BASE and ending at | |
164 | END, with each element of SIZE bytes. The SWAP_TYPE is the callback | |
165 | function used to swap elements, and CMP is the function used to compare | |
166 | elements. */ | |
167 | static void | |
168 | heapsort_r (void *base, void *end, size_t size, enum swap_type_t swap_type, | |
169 | __compar_d_fn_t cmp, void *arg) | |
170 | { | |
171 | const size_t count = ((uintptr_t) end - (uintptr_t) base) / size; | |
172 | ||
173 | if (count < 2) | |
174 | return; | |
175 | ||
176 | size_t n = count - 1; | |
177 | ||
178 | /* Build the binary heap, largest value at the base[0]. */ | |
179 | heapify (base, size, n, swap_type, cmp, arg); | |
180 | ||
181 | /* On each iteration base[0:n] is the binary heap, while base[n:count] | |
182 | is sorted. */ | |
183 | while (n > 0) | |
184 | { | |
185 | do_swap (base, base + (n * size), size, swap_type); | |
186 | n--; | |
187 | siftdown (base, size, 0, n, swap_type, cmp, arg); | |
188 | } | |
189 | } | |
28f540f4 | 190 | |
a035a985 AZ |
191 | static inline void |
192 | insertion_sort_qsort_partitions (void *const pbase, size_t total_elems, | |
193 | size_t size, enum swap_type_t swap_type, | |
194 | __compar_d_fn_t cmp, void *arg) | |
195 | { | |
196 | char *base_ptr = (char *) pbase; | |
197 | char *const end_ptr = &base_ptr[size * (total_elems - 1)]; | |
198 | char *tmp_ptr = base_ptr; | |
199 | #define min(x, y) ((x) < (y) ? (x) : (y)) | |
200 | const size_t max_thresh = MAX_THRESH * size; | |
201 | char *thresh = min(end_ptr, base_ptr + max_thresh); | |
202 | char *run_ptr; | |
203 | ||
204 | /* Find smallest element in first threshold and place it at the | |
205 | array's beginning. This is the smallest array element, | |
206 | and the operation speeds up insertion sort's inner loop. */ | |
207 | ||
208 | for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size) | |
209 | if (cmp (run_ptr, tmp_ptr, arg) < 0) | |
210 | tmp_ptr = run_ptr; | |
211 | ||
212 | if (tmp_ptr != base_ptr) | |
213 | do_swap (tmp_ptr, base_ptr, size, swap_type); | |
214 | ||
215 | /* Insertion sort, running from left-hand-side up to right-hand-side. */ | |
216 | ||
217 | run_ptr = base_ptr + size; | |
218 | while ((run_ptr += size) <= end_ptr) | |
219 | { | |
220 | tmp_ptr = run_ptr - size; | |
221 | while (cmp (run_ptr, tmp_ptr, arg) < 0) | |
222 | tmp_ptr -= size; | |
223 | ||
224 | tmp_ptr += size; | |
225 | if (tmp_ptr != run_ptr) | |
226 | { | |
227 | char *trav; | |
228 | ||
229 | trav = run_ptr + size; | |
230 | while (--trav >= run_ptr) | |
231 | { | |
232 | char c = *trav; | |
233 | char *hi, *lo; | |
234 | ||
235 | for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo) | |
236 | *hi = *lo; | |
237 | *hi = c; | |
238 | } | |
239 | } | |
240 | } | |
241 | } | |
242 | ||
28f540f4 RM |
243 | /* Order size using quicksort. This implementation incorporates |
244 | four optimizations discussed in Sedgewick: | |
245 | ||
6d52618b UD |
246 | 1. Non-recursive, using an explicit stack of pointer that store the |
247 | next array partition to sort. To save time, this maximum amount | |
061d137b UD |
248 | of space required to store an array of SIZE_MAX is allocated on the |
249 | stack. Assuming a 32-bit (64 bit) integer for size_t, this needs | |
250 | only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes). | |
251 | Pretty cheap, actually. | |
28f540f4 RM |
252 | |
253 | 2. Chose the pivot element using a median-of-three decision tree. | |
6d52618b | 254 | This reduces the probability of selecting a bad pivot value and |
28f540f4 RM |
255 | eliminates certain extraneous comparisons. |
256 | ||
257 | 3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving | |
6d52618b | 258 | insertion sort to order the MAX_THRESH items within each partition. |
28f540f4 | 259 | This is a big win, since insertion sort is faster for small, mostly |
6d52618b | 260 | sorted array segments. |
28f540f4 RM |
261 | |
262 | 4. The larger of the two sub-partitions is always pushed onto the | |
263 | stack first, with the algorithm then concentrating on the | |
061d137b | 264 | smaller partition. This *guarantees* no more than log (total_elems) |
28f540f4 RM |
265 | stack size is needed (actually O(1) in this case)! */ |
266 | ||
267 | void | |
061d137b | 268 | _quicksort (void *const pbase, size_t total_elems, size_t size, |
e458144c | 269 | __compar_d_fn_t cmp, void *arg) |
28f540f4 | 270 | { |
2e09a79a | 271 | char *base_ptr = (char *) pbase; |
28f540f4 | 272 | |
7cc27f44 | 273 | const size_t max_thresh = MAX_THRESH * size; |
28f540f4 | 274 | |
274a46c9 | 275 | if (total_elems <= 1) |
28f540f4 RM |
276 | /* Avoid lossage with unsigned arithmetic below. */ |
277 | return; | |
278 | ||
21d30c77 AZ |
279 | enum swap_type_t swap_type; |
280 | if (is_aligned (pbase, size, 8)) | |
281 | swap_type = SWAP_WORDS_64; | |
282 | else if (is_aligned (pbase, size, 4)) | |
283 | swap_type = SWAP_WORDS_32; | |
284 | else | |
285 | swap_type = SWAP_BYTES; | |
286 | ||
274a46c9 AZ |
287 | /* Maximum depth before quicksort switches to heapsort. */ |
288 | size_t depth = 2 * (sizeof (size_t) * CHAR_BIT - 1 | |
289 | - __builtin_clzl (total_elems)); | |
290 | ||
28f540f4 RM |
291 | if (total_elems > MAX_THRESH) |
292 | { | |
293 | char *lo = base_ptr; | |
294 | char *hi = &lo[size * (total_elems - 1)]; | |
28f540f4 | 295 | stack_node stack[STACK_SIZE]; |
274a46c9 | 296 | stack_node *top = push (stack, NULL, NULL, depth); |
28f540f4 | 297 | |
d097f3c7 | 298 | while (stack < top) |
28f540f4 | 299 | { |
274a46c9 AZ |
300 | if (depth == 0) |
301 | { | |
302 | heapsort_r (lo, hi, size, swap_type, cmp, arg); | |
303 | top = pop (top, &lo, &hi, &depth); | |
304 | continue; | |
305 | } | |
306 | ||
28f540f4 RM |
307 | char *left_ptr; |
308 | char *right_ptr; | |
309 | ||
28f540f4 | 310 | /* Select median value from among LO, MID, and HI. Rearrange |
6d52618b UD |
311 | LO and HI so the three values are sorted. This lowers the |
312 | probability of picking a pathological pivot value and | |
061d137b UD |
313 | skips a comparison for both the LEFT_PTR and RIGHT_PTR in |
314 | the while loops. */ | |
28f540f4 RM |
315 | |
316 | char *mid = lo + size * ((hi - lo) / size >> 1); | |
317 | ||
e458144c | 318 | if ((*cmp) ((void *) mid, (void *) lo, arg) < 0) |
21d30c77 | 319 | do_swap (mid, lo, size, swap_type); |
e458144c | 320 | if ((*cmp) ((void *) hi, (void *) mid, arg) < 0) |
21d30c77 | 321 | do_swap (mid, hi, size, swap_type); |
6d52618b | 322 | else |
28f540f4 | 323 | goto jump_over; |
e458144c | 324 | if ((*cmp) ((void *) mid, (void *) lo, arg) < 0) |
21d30c77 | 325 | do_swap (mid, lo, size, swap_type); |
28f540f4 | 326 | jump_over:; |
28f540f4 RM |
327 | |
328 | left_ptr = lo + size; | |
6d52618b | 329 | right_ptr = hi - size; |
28f540f4 | 330 | |
6d52618b UD |
331 | /* Here's the famous ``collapse the walls'' section of quicksort. |
332 | Gotta like those tight inner loops! They are the main reason | |
28f540f4 | 333 | that this algorithm runs much faster than others. */ |
6d52618b | 334 | do |
28f540f4 | 335 | { |
e458144c | 336 | while ((*cmp) ((void *) left_ptr, (void *) mid, arg) < 0) |
28f540f4 RM |
337 | left_ptr += size; |
338 | ||
e458144c | 339 | while ((*cmp) ((void *) mid, (void *) right_ptr, arg) < 0) |
28f540f4 RM |
340 | right_ptr -= size; |
341 | ||
6d52618b | 342 | if (left_ptr < right_ptr) |
28f540f4 | 343 | { |
21d30c77 | 344 | do_swap (left_ptr, right_ptr, size, swap_type); |
fa8d436c UD |
345 | if (mid == left_ptr) |
346 | mid = right_ptr; | |
347 | else if (mid == right_ptr) | |
348 | mid = left_ptr; | |
28f540f4 RM |
349 | left_ptr += size; |
350 | right_ptr -= size; | |
351 | } | |
6d52618b | 352 | else if (left_ptr == right_ptr) |
28f540f4 RM |
353 | { |
354 | left_ptr += size; | |
355 | right_ptr -= size; | |
356 | break; | |
357 | } | |
6d52618b | 358 | } |
28f540f4 RM |
359 | while (left_ptr <= right_ptr); |
360 | ||
361 | /* Set up pointers for next iteration. First determine whether | |
6d52618b | 362 | left and right partitions are below the threshold size. If so, |
28f540f4 RM |
363 | ignore one or both. Otherwise, push the larger partition's |
364 | bounds on the stack and continue sorting the smaller one. */ | |
365 | ||
366 | if ((size_t) (right_ptr - lo) <= max_thresh) | |
367 | { | |
368 | if ((size_t) (hi - left_ptr) <= max_thresh) | |
369 | /* Ignore both small partitions. */ | |
274a46c9 | 370 | top = pop (top, &lo, &hi, &depth); |
28f540f4 | 371 | else |
6d52618b | 372 | /* Ignore small left partition. */ |
28f540f4 RM |
373 | lo = left_ptr; |
374 | } | |
375 | else if ((size_t) (hi - left_ptr) <= max_thresh) | |
376 | /* Ignore small right partition. */ | |
377 | hi = right_ptr; | |
378 | else if ((right_ptr - lo) > (hi - left_ptr)) | |
6d52618b | 379 | { |
28f540f4 | 380 | /* Push larger left partition indices. */ |
274a46c9 | 381 | top = push (top, lo, right_ptr, depth - 1); |
28f540f4 RM |
382 | lo = left_ptr; |
383 | } | |
384 | else | |
6d52618b | 385 | { |
28f540f4 | 386 | /* Push larger right partition indices. */ |
274a46c9 | 387 | top = push (top, left_ptr, hi, depth - 1); |
28f540f4 RM |
388 | hi = right_ptr; |
389 | } | |
390 | } | |
391 | } | |
392 | ||
393 | /* Once the BASE_PTR array is partially sorted by quicksort the rest | |
6d52618b UD |
394 | is completely sorted using insertion sort, since this is efficient |
395 | for partitions below MAX_THRESH size. BASE_PTR points to the beginning | |
28f540f4 RM |
396 | of the array to sort, and END_PTR points at the very last element in |
397 | the array (*not* one beyond it!). */ | |
a035a985 AZ |
398 | insertion_sort_qsort_partitions (pbase, total_elems, size, swap_type, cmp, |
399 | arg); | |
28f540f4 | 400 | } |