stdlib/qsort.c

   1 /* Copyright (C) 1991-2023 Free Software Foundation, Inc.
   2    This file is part of the GNU C Library.
   3
   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.
   8
   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.
  13
  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/>.  */
  17
  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>
  24 #include <memswap.h>
  25 #include <stdlib.h>
  26 #include <string.h>
  27 #include <stdbool.h>
  28
  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 }
  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. */
  97 typedef struct
  98   {
  99     char *lo;
 100     char *hi;
 101   } stack_node;
 102
 103 /* The next 4 #defines implement a very fast in-line stack abstraction. */
 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 #define STACK_SIZE      (CHAR_BIT * sizeof (size_t))
 109 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
 110 #define POP(low, high)  ((void) (--top, (low = top->lo), (high = top->hi)))
 111 #define STACK_NOT_EMPTY (stack < top)
 112
 113
 114 /* Order size using quicksort.  This implementation incorporates
 115    four optimizations discussed in Sedgewick:
 116
 117    1. Non-recursive, using an explicit stack of pointer that store the
 118       next array partition to sort.  To save time, this maximum amount
 119       of space required to store an array of SIZE_MAX is allocated on the
 120       stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
 121       only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
 122       Pretty cheap, actually.
 123
 124    2. Chose the pivot element using a median-of-three decision tree.
 125       This reduces the probability of selecting a bad pivot value and
 126       eliminates certain extraneous comparisons.
 127
 128    3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
 129       insertion sort to order the MAX_THRESH items within each partition.
 130       This is a big win, since insertion sort is faster for small, mostly
 131       sorted array segments.
 132
 133    4. The larger of the two sub-partitions is always pushed onto the
 134       stack first, with the algorithm then concentrating on the
 135       smaller partition.  This *guarantees* no more than log (total_elems)
 136       stack size is needed (actually O(1) in this case)!  */
 137
 138 void
 139 _quicksort (void *const pbase, size_t total_elems, size_t size,
 140             __compar_d_fn_t cmp, void *arg)
 141 {
 142   char *base_ptr = (char *) pbase;
 143
 144   const size_t max_thresh = MAX_THRESH * size;
 145
 146   if (total_elems == 0)
 147     /* Avoid lossage with unsigned arithmetic below.  */
 148     return;
 149
 150   enum swap_type_t swap_type;
 151   if (is_aligned (pbase, size, 8))
 152     swap_type = SWAP_WORDS_64;
 153   else if (is_aligned (pbase, size, 4))
 154     swap_type = SWAP_WORDS_32;
 155   else
 156     swap_type = SWAP_BYTES;
 157
 158   if (total_elems > MAX_THRESH)
 159     {
 160       char *lo = base_ptr;
 161       char *hi = &lo[size * (total_elems - 1)];
 162       stack_node stack[STACK_SIZE];
 163       stack_node *top = stack;
 164
 165       PUSH (NULL, NULL);
 166
 167       while (STACK_NOT_EMPTY)
 168         {
 169           char *left_ptr;
 170           char *right_ptr;
 171
 172           /* Select median value from among LO, MID, and HI. Rearrange
 173              LO and HI so the three values are sorted. This lowers the
 174              probability of picking a pathological pivot value and
 175              skips a comparison for both the LEFT_PTR and RIGHT_PTR in
 176              the while loops. */
 177
 178           char *mid = lo + size * ((hi - lo) / size >> 1);
 179
 180           if ((*cmp) ((void *) mid, (void *) lo, arg) < 0)
 181             do_swap (mid, lo, size, swap_type);
 182           if ((*cmp) ((void *) hi, (void *) mid, arg) < 0)
 183             do_swap (mid, hi, size, swap_type);
 184           else
 185             goto jump_over;
 186           if ((*cmp) ((void *) mid, (void *) lo, arg) < 0)
 187             do_swap (mid, lo, size, swap_type);
 188         jump_over:;
 189
 190           left_ptr  = lo + size;
 191           right_ptr = hi - size;
 192
 193           /* Here's the famous ``collapse the walls'' section of quicksort.
 194              Gotta like those tight inner loops!  They are the main reason
 195              that this algorithm runs much faster than others. */
 196           do
 197             {
 198               while ((*cmp) ((void *) left_ptr, (void *) mid, arg) < 0)
 199                 left_ptr += size;
 200
 201               while ((*cmp) ((void *) mid, (void *) right_ptr, arg) < 0)
 202                 right_ptr -= size;
 203
 204               if (left_ptr < right_ptr)
 205                 {
 206                   do_swap (left_ptr, right_ptr, size, swap_type);
 207                   if (mid == left_ptr)
 208                     mid = right_ptr;
 209                   else if (mid == right_ptr)
 210                     mid = left_ptr;
 211                   left_ptr += size;
 212                   right_ptr -= size;
 213                 }
 214               else if (left_ptr == right_ptr)
 215                 {
 216                   left_ptr += size;
 217                   right_ptr -= size;
 218                   break;
 219                 }
 220             }
 221           while (left_ptr <= right_ptr);
 222
 223           /* Set up pointers for next iteration.  First determine whether
 224              left and right partitions are below the threshold size.  If so,
 225              ignore one or both.  Otherwise, push the larger partition's
 226              bounds on the stack and continue sorting the smaller one. */
 227
 228           if ((size_t) (right_ptr - lo) <= max_thresh)
 229             {
 230               if ((size_t) (hi - left_ptr) <= max_thresh)
 231                 /* Ignore both small partitions. */
 232                 POP (lo, hi);
 233               else
 234                 /* Ignore small left partition. */
 235                 lo = left_ptr;
 236             }
 237           else if ((size_t) (hi - left_ptr) <= max_thresh)
 238             /* Ignore small right partition. */
 239             hi = right_ptr;
 240           else if ((right_ptr - lo) > (hi - left_ptr))
 241             {
 242               /* Push larger left partition indices. */
 243               PUSH (lo, right_ptr);
 244               lo = left_ptr;
 245             }
 246           else
 247             {
 248               /* Push larger right partition indices. */
 249               PUSH (left_ptr, hi);
 250               hi = right_ptr;
 251             }
 252         }
 253     }
 254
 255   /* Once the BASE_PTR array is partially sorted by quicksort the rest
 256      is completely sorted using insertion sort, since this is efficient
 257      for partitions below MAX_THRESH size. BASE_PTR points to the beginning
 258      of the array to sort, and END_PTR points at the very last element in
 259      the array (*not* one beyond it!). */
 260
 261 #define min(x, y) ((x) < (y) ? (x) : (y))
 262
 263   {
 264     char *const end_ptr = &base_ptr[size * (total_elems - 1)];
 265     char *tmp_ptr = base_ptr;
 266     char *thresh = min(end_ptr, base_ptr + max_thresh);
 267     char *run_ptr;
 268
 269     /* Find smallest element in first threshold and place it at the
 270        array's beginning.  This is the smallest array element,
 271        and the operation speeds up insertion sort's inner loop. */
 272
 273     for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
 274       if ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, arg) < 0)
 275         tmp_ptr = run_ptr;
 276
 277     if (tmp_ptr != base_ptr)
 278       do_swap (tmp_ptr, base_ptr, size, swap_type);
 279
 280     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
 281
 282     run_ptr = base_ptr + size;
 283     while ((run_ptr += size) <= end_ptr)
 284       {
 285         tmp_ptr = run_ptr - size;
 286         while ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, arg) < 0)
 287           tmp_ptr -= size;
 288
 289         tmp_ptr += size;
 290         if (tmp_ptr != run_ptr)
 291           {
 292             char *trav;
 293
 294             trav = run_ptr + size;
 295             while (--trav >= run_ptr)
 296               {
 297                 char c = *trav;
 298                 char *hi, *lo;
 299
 300                 for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
 301                   *hi = *lo;
 302                 *hi = c;
 303               }
 304           }
 305       }
 306   }
 307 }