guile18/libguile/quicksort.i.c

   1 /* The routine quicksort was extracted from the GNU C Library qsort.c
   2    written by Douglas C. Schmidt (schmidt@ics.uci.edu)
   3    and adapted to guile by adding an extra pointer less
   4    to quicksort by Roland Orre <orre@nada.kth.se>.
   5
   6    The reason to do this instead of using the library function qsort
   7    was to avoid dependency of the ANSI-C extensions for local functions
   8    and also to avoid obscure pool based solutions.
   9
  10    This sorting routine is not much more efficient than the stable
  11    version but doesn't consume extra memory.
  12  */
  13
  14 #define SWAP(a, b) do { const SCM _tmp = a; a = b; b = _tmp; } while (0)
  15
  16
  17 /* Order using quicksort.  This implementation incorporates four
  18    optimizations discussed in Sedgewick:
  19
  20    1. Non-recursive, using an explicit stack of pointer that store the next
  21    array partition to sort.  To save time, this maximum amount of space
  22    required to store an array of MAX_SIZE_T is allocated on the stack.
  23    Assuming a bit width of 32 bits for size_t, this needs only
  24    32 * sizeof (stack_node) == 128 bytes.  Pretty cheap, actually.
  25
  26    2. Chose the pivot element using a median-of-three decision tree.  This
  27    reduces the probability of selecting a bad pivot value and eliminates
  28    certain extraneous comparisons.
  29
  30    3. Only quicksorts NR_ELEMS / MAX_THRESH partitions, leaving insertion sort
  31    to order the MAX_THRESH items within each partition.  This is a big win,
  32    since insertion sort is faster for small, mostly sorted array segments.
  33
  34    4. The larger of the two sub-partitions is always pushed onto the
  35    stack first, with the algorithm then concentrating on the
  36    smaller partition.  This *guarantees* no more than log (n)
  37    stack size is needed (actually O(1) in this case)!  */
  38
  39
  40 /* Discontinue quicksort algorithm when partition gets below this size.
  41  * This particular magic number was chosen to work best on a Sun 4/260. */
  42 #define MAX_THRESH 4
  43
  44
  45 /* Inline stack abstraction:  The stack size for quicksorting at most as many
  46  * elements as can be given by a value of type size_t is, as described above,
  47  * log (MAX_SIZE_T), which is the number of bits of size_t.  More accurately,
  48  * we would only need ceil (log (MAX_SIZE_T / MAX_THRESH)), but this is
  49  * ignored below. */
  50
  51 #define STACK_SIZE       (8 * sizeof (size_t))  /* assume 8 bit char */
  52 #define PUSH(low, high)  ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
  53 #define POP(low, high)   ((void) (--top, (low = top->lo), (high = top->hi)))
  54 #define STACK_NOT_EMPTY  (stack < top)
  55
  56 static void
  57 NAME (SCM *const base_ptr, size_t nr_elems, INC_PARAM
  58       scm_t_trampoline_2 cmp, SCM less)
  59 {
  60   /* Stack node declarations used to store unfulfilled partition obligations. */
  61   typedef struct {
  62     size_t lo;
  63     size_t hi;
  64   } stack_node;
  65
  66   static const char s_buggy_less[] = "buggy less predicate used when sorting";
  67
  68 #define ELT(i) base_ptr[(i)*INC]
  69
  70   if (nr_elems == 0)
  71     /* Avoid lossage with unsigned arithmetic below.  */
  72     return;
  73
  74   if (nr_elems > MAX_THRESH)
  75     {
  76       size_t lo = 0;
  77       size_t hi = nr_elems-1;
  78
  79       stack_node stack[STACK_SIZE];
  80       stack_node *top = stack + 1;
  81
  82       while (STACK_NOT_EMPTY)
  83         {
  84           size_t left;
  85           size_t right;
  86           size_t mid = lo + (hi - lo) / 2;
  87           SCM pivot;
  88
  89           /* Select median value from among LO, MID, and HI. Rearrange
  90              LO and HI so the three values are sorted. This lowers the
  91              probability of picking a pathological pivot value and
  92              skips a comparison for both the left and right. */
  93
  94           SCM_TICK;
  95
  96           if (scm_is_true ((*cmp) (less, ELT(mid), ELT(lo))))
  97             SWAP (ELT(mid), ELT(lo));
  98           if (scm_is_true ((*cmp) (less, ELT(hi), ELT(mid))))
  99             SWAP (ELT(mid), ELT(hi));
 100           else
 101             goto jump_over;
 102           if (scm_is_true ((*cmp) (less, ELT(mid), ELT(lo))))
 103             SWAP (ELT(mid), ELT(lo));
 104         jump_over:;
 105
 106           pivot = ELT(mid);
 107           left = lo + 1;
 108           right = hi - 1;
 109
 110           /* Here's the famous ``collapse the walls'' section of quicksort.
 111              Gotta like those tight inner loops!  They are the main reason
 112              that this algorithm runs much faster than others. */
 113           do
 114             {
 115               while (scm_is_true ((*cmp) (less, ELT(left), pivot)))
 116                 {
 117                   left += 1;
 118                   /* The comparison predicate may be buggy */
 119                   if (left > hi)
 120                     scm_misc_error (NULL, s_buggy_less, SCM_EOL);
 121                 }
 122
 123               while (scm_is_true ((*cmp) (less, pivot, ELT(right))))
 124                 {
 125                   right -= 1;
 126                   /* The comparison predicate may be buggy */
 127                   if (right < lo)
 128                     scm_misc_error (NULL, s_buggy_less, SCM_EOL);
 129                 }
 130
 131               if (left < right)
 132                 {
 133                   SWAP (ELT(left), ELT(right));
 134                   left += 1;
 135                   right -= 1;
 136                 }
 137               else if (left == right)
 138                 {
 139                   left += 1;
 140                   right -= 1;
 141                   break;
 142                 }
 143             }
 144           while (left <= right);
 145
 146           /* Set up pointers for next iteration.  First determine whether
 147              left and right partitions are below the threshold size.  If so,
 148              ignore one or both.  Otherwise, push the larger partition's
 149              bounds on the stack and continue sorting the smaller one. */
 150
 151           if ((size_t) (right - lo) <= MAX_THRESH)
 152             {
 153               if ((size_t) (hi - left) <= MAX_THRESH)
 154                 /* Ignore both small partitions. */
 155                 POP (lo, hi);
 156               else
 157                 /* Ignore small left partition. */
 158                 lo = left;
 159             }
 160           else if ((size_t) (hi - left) <= MAX_THRESH)
 161             /* Ignore small right partition. */
 162             hi = right;
 163           else if ((right - lo) > (hi - left))
 164             {
 165               /* Push larger left partition indices. */
 166               PUSH (lo, right);
 167               lo = left;
 168             }
 169           else
 170             {
 171               /* Push larger right partition indices. */
 172               PUSH (left, hi);
 173               hi = right;
 174             }
 175         }
 176     }
 177
 178   /* Once the BASE_PTR array is partially sorted by quicksort the rest is
 179      completely sorted using insertion sort, since this is efficient for
 180      partitions below MAX_THRESH size. BASE_PTR points to the beginning of the
 181      array to sort, and END idexes the very last element in the array (*not*
 182      one beyond it!). */
 183
 184   {
 185     size_t tmp = 0;
 186     size_t end = nr_elems-1;
 187     size_t thresh = min (end, MAX_THRESH);
 188     size_t run;
 189
 190     /* Find smallest element in first threshold and place it at the
 191        array's beginning.  This is the smallest array element,
 192        and the operation speeds up insertion sort's inner loop. */
 193
 194     for (run = tmp + 1; run <= thresh; run += 1)
 195       if (scm_is_true ((*cmp) (less, ELT(run), ELT(tmp))))
 196         tmp = run;
 197
 198     if (tmp != 0)
 199       SWAP (ELT(tmp), ELT(0));
 200
 201     /* Insertion sort, running from left-hand-side up to right-hand-side.  */
 202
 203     run = 1;
 204     while (++run <= end)
 205       {
 206         SCM_TICK;
 207
 208         tmp = run - 1;
 209         while (scm_is_true ((*cmp) (less, ELT(run), ELT(tmp))))
 210           {
 211             /* The comparison predicate may be buggy */
 212             if (tmp == 0)
 213               scm_misc_error (NULL, s_buggy_less, SCM_EOL);
 214
 215             tmp -= 1;
 216           }
 217
 218         tmp += 1;
 219         if (tmp != run)
 220           {
 221             SCM to_insert = ELT(run);
 222             size_t hi, lo;
 223
 224             for (hi = lo = run; --lo >= tmp; hi = lo)
 225               ELT(hi) = ELT(lo);
 226             ELT(hi) = to_insert;
 227           }
 228       }
 229   }
 230 }
 231
 232 #undef SWAP
 233 #undef MAX_THRESH
 234 #undef STACK_SIZE
 235 #undef PUSH
 236 #undef POP
 237 #undef STACK_NOT_EMPTY
 238 #undef ELT
 239
 240 #undef NAME
 241 #undef INC_PARAM
 242 #undef INC
 243