1 /* Copyright (C) 1991-2023 Free Software Foundation, Inc.
2 This file is part of the GNU C Library.
4 The GNU C Library is free software; you can redistribute it and/or
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.
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
12 Lesser General Public License for more details.
14 You should have received a copy of the GNU Lesser General Public
15 License along with the GNU C Library; if not, see
16 <https://www.gnu.org/licenses/>. */
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. */
29 /* Swap SIZE bytes between addresses A and B. These helpers are provided
30 along the generic one as an optimization. */
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
47 is_aligned (const void *base
, size_t size
, size_t wordsize
)
49 return (((uintptr_t) base
| size
) & (wordsize
- 1)) == 0;
53 swap_words_64 (void * restrict a
, void * restrict b
, size_t n
)
55 typedef uint64_t __attribute__ ((__may_alias__
)) u64_alias_t
;
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
;
66 swap_words_32 (void * restrict a
, void * restrict b
, size_t n
)
68 typedef uint32_t __attribute__ ((__may_alias__
)) u32_alias_t
;
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
;
78 /* Replace the indirect call with a serie of if statements. It should help
79 the branch predictor. */
81 do_swap (void * restrict a
, void * restrict b
, size_t size
,
82 enum swap_type_t swap_type
)
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
);
89 __memswap (a
, b
, size
);
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. */
96 /* Stack node declarations used to store unfulfilled partition obligations. */
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). */
108 enum { STACK_SIZE
= CHAR_BIT
* sizeof (size_t) };
110 static inline stack_node
*
111 push (stack_node
*top
, char *lo
, char *hi
, size_t depth
)
119 static inline stack_node
*
120 pop (stack_node
*top
, char **lo
, char **hi
, size_t *depth
)
129 /* NB: N is inclusive bound for BASE. */
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
)
137 if (j
< n
&& cmp (base
+ (j
* size
), base
+ ((j
+ 1) * size
), arg
) < 0)
140 if (cmp (base
+ (k
* size
), base
+ (j
* size
), arg
) >= 0)
143 do_swap (base
+ (size
* j
), base
+ (k
* size
), size
, swap_type
);
149 heapify (void *base
, size_t size
, size_t n
, enum swap_type_t swap_type
,
150 __compar_d_fn_t cmp
, void *arg
)
155 siftdown (base
, size
, k
, n
, swap_type
, cmp
, arg
);
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
168 heapsort_r (void *base
, void *end
, size_t size
, enum swap_type_t swap_type
,
169 __compar_d_fn_t cmp
, void *arg
)
171 const size_t count
= ((uintptr_t) end
- (uintptr_t) base
) / size
;
176 size_t n
= count
- 1;
178 /* Build the binary heap, largest value at the base[0]. */
179 heapify (base
, size
, n
, swap_type
, cmp
, arg
);
181 /* On each iteration base[0:n] is the binary heap, while base[n:count]
185 do_swap (base
, base
+ (n
* size
), size
, swap_type
);
187 siftdown (base
, size
, 0, n
, swap_type
, cmp
, arg
);
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
)
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
);
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. */
208 for (run_ptr
= tmp_ptr
+ size
; run_ptr
<= thresh
; run_ptr
+= size
)
209 if (cmp (run_ptr
, tmp_ptr
, arg
) < 0)
212 if (tmp_ptr
!= base_ptr
)
213 do_swap (tmp_ptr
, base_ptr
, size
, swap_type
);
215 /* Insertion sort, running from left-hand-side up to right-hand-side. */
217 run_ptr
= base_ptr
+ size
;
218 while ((run_ptr
+= size
) <= end_ptr
)
220 tmp_ptr
= run_ptr
- size
;
221 while (cmp (run_ptr
, tmp_ptr
, arg
) < 0)
225 if (tmp_ptr
!= run_ptr
)
229 trav
= run_ptr
+ size
;
230 while (--trav
>= run_ptr
)
235 for (hi
= lo
= trav
; (lo
-= size
) >= tmp_ptr
; hi
= lo
)
243 /* Order size using quicksort. This implementation incorporates
244 four optimizations discussed in Sedgewick:
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
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.
253 2. Chose the pivot element using a median-of-three decision tree.
254 This reduces the probability of selecting a bad pivot value and
255 eliminates certain extraneous comparisons.
257 3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
258 insertion sort to order the MAX_THRESH items within each partition.
259 This is a big win, since insertion sort is faster for small, mostly
260 sorted array segments.
262 4. The larger of the two sub-partitions is always pushed onto the
263 stack first, with the algorithm then concentrating on the
264 smaller partition. This *guarantees* no more than log (total_elems)
265 stack size is needed (actually O(1) in this case)! */
268 _quicksort (void *const pbase
, size_t total_elems
, size_t size
,
269 __compar_d_fn_t cmp
, void *arg
)
271 char *base_ptr
= (char *) pbase
;
273 const size_t max_thresh
= MAX_THRESH
* size
;
275 if (total_elems
<= 1)
276 /* Avoid lossage with unsigned arithmetic below. */
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
;
285 swap_type
= SWAP_BYTES
;
287 /* Maximum depth before quicksort switches to heapsort. */
288 size_t depth
= 2 * (sizeof (size_t) * CHAR_BIT
- 1
289 - __builtin_clzl (total_elems
));
291 if (total_elems
> MAX_THRESH
)
294 char *hi
= &lo
[size
* (total_elems
- 1)];
295 stack_node stack
[STACK_SIZE
];
296 stack_node
*top
= push (stack
, NULL
, NULL
, depth
);
302 heapsort_r (lo
, hi
, size
, swap_type
, cmp
, arg
);
303 top
= pop (top
, &lo
, &hi
, &depth
);
310 /* Select median value from among LO, MID, and HI. Rearrange
311 LO and HI so the three values are sorted. This lowers the
312 probability of picking a pathological pivot value and
313 skips a comparison for both the LEFT_PTR and RIGHT_PTR in
316 char *mid
= lo
+ size
* ((hi
- lo
) / size
>> 1);
318 if ((*cmp
) ((void *) mid
, (void *) lo
, arg
) < 0)
319 do_swap (mid
, lo
, size
, swap_type
);
320 if ((*cmp
) ((void *) hi
, (void *) mid
, arg
) < 0)
321 do_swap (mid
, hi
, size
, swap_type
);
324 if ((*cmp
) ((void *) mid
, (void *) lo
, arg
) < 0)
325 do_swap (mid
, lo
, size
, swap_type
);
328 left_ptr
= lo
+ size
;
329 right_ptr
= hi
- size
;
331 /* Here's the famous ``collapse the walls'' section of quicksort.
332 Gotta like those tight inner loops! They are the main reason
333 that this algorithm runs much faster than others. */
336 while ((*cmp
) ((void *) left_ptr
, (void *) mid
, arg
) < 0)
339 while ((*cmp
) ((void *) mid
, (void *) right_ptr
, arg
) < 0)
342 if (left_ptr
< right_ptr
)
344 do_swap (left_ptr
, right_ptr
, size
, swap_type
);
347 else if (mid
== right_ptr
)
352 else if (left_ptr
== right_ptr
)
359 while (left_ptr
<= right_ptr
);
361 /* Set up pointers for next iteration. First determine whether
362 left and right partitions are below the threshold size. If so,
363 ignore one or both. Otherwise, push the larger partition's
364 bounds on the stack and continue sorting the smaller one. */
366 if ((size_t) (right_ptr
- lo
) <= max_thresh
)
368 if ((size_t) (hi
- left_ptr
) <= max_thresh
)
369 /* Ignore both small partitions. */
370 top
= pop (top
, &lo
, &hi
, &depth
);
372 /* Ignore small left partition. */
375 else if ((size_t) (hi
- left_ptr
) <= max_thresh
)
376 /* Ignore small right partition. */
378 else if ((right_ptr
- lo
) > (hi
- left_ptr
))
380 /* Push larger left partition indices. */
381 top
= push (top
, lo
, right_ptr
, depth
- 1);
386 /* Push larger right partition indices. */
387 top
= push (top
, left_ptr
, hi
, depth
- 1);
393 /* Once the BASE_PTR array is partially sorted by quicksort the rest
394 is completely sorted using insertion sort, since this is efficient
395 for partitions below MAX_THRESH size. BASE_PTR points to the beginning
396 of the array to sort, and END_PTR points at the very last element in
397 the array (*not* one beyond it!). */
398 insertion_sort_qsort_partitions (pbase
, total_elems
, size
, swap_type
, cmp
,