Use David Musser's introsort algorithm to fall back to heapsort(3)

when the recursion depth reaches 2*lg(n + 1).  This avoids quicksort's
quadratic behavior for pathological input without appreciably
changing the average run time.

2 files changed, 45 insertions, 14 deletions

diff --git a/src/lib/libc/stdlib/heapsort.c b/src/lib/libc/stdlib/heapsort.c
index 878567729e..f1db2205b0 100644
--- a/src/lib/libc/stdlib/heapsort.c
+++ b/src/lib/libc/stdlib/heapsort.c
@@ -1,4 +1,4 @@
-/*      $OpenBSD: heapsort.c,v 1.10 2016/10/22 19:19:34 tb Exp $ */
+/*      $OpenBSD: heapsort.c,v 1.11 2017/05/20 12:48:56 millert Exp $ */
 /*-
  * Copyright (c) 1991, 1993
  *      The Regents of the University of California.  All rights reserved.
@@ -172,3 +172,4 @@ heapsort(void *vbase, size_t nmemb, size_t size,
         free(k);
         return (0);
 }
+DEF_WEAK(heapsort);
diff --git a/src/lib/libc/stdlib/qsort.c b/src/lib/libc/stdlib/qsort.c
index c97979ef01..4fe54eb5ff 100644
--- a/src/lib/libc/stdlib/qsort.c
+++ b/src/lib/libc/stdlib/qsort.c
@@ -1,4 +1,4 @@
-/*      $OpenBSD: qsort.c,v 1.15 2017/05/17 16:58:20 millert Exp $ */
+/*      $OpenBSD: qsort.c,v 1.16 2017/05/20 12:48:56 millert Exp $ */
 /*-
  * Copyright (c) 1992, 1993
  *      The Regents of the University of California.  All rights reserved.
@@ -38,6 +38,18 @@ static __inline void	 swapfunc(char *, char *, size_t, int);
 
 /*
  * Qsort routine from Bentley & McIlroy's "Engineering a Sort Function".
+ *
+ * This version differs from Bentley & McIlroy in the following ways:
+ *   1. The partition value is swapped into a[0] instead of being
+ *      stored out of line.
+ *
+ *   2. It uses David Musser's introsort algorithm to fall back to
+ *      heapsort(3) when the recursion depth reaches 2*lg(n + 1).
+ *      This avoids quicksort's quadratic behavior for pathological
+ *      input without appreciably changing the average run time.
+ *
+ *   3. Tail recursion is eliminated when sorting the larger of two
+ *      subpartitions to save stack space.
  */
 #define swapcode(TYPE, parmi, parmj, n) {               \
         size_t i = (n) / sizeof (TYPE);                 \
@@ -80,15 +92,20 @@ med3(char *a, char *b, char *c, int (*cmp)(const void *, const void *))
               :(cmp(b, c) > 0 ? b : (cmp(a, c) < 0 ? a : c ));
 }
 
-void
-qsort(void *aa, size_t n, size_t es, int (*cmp)(const void *, const void *))
+static void
+introsort(char *a, size_t n, size_t es, size_t maxdepth,
+    int (*cmp)(const void *, const void *))
 {
         char *pa, *pb, *pc, *pd, *pl, *pm, *pn;
         int cmp_result, swaptype;
-        size_t d, r, s;
-        char *a = aa;
+        size_t r, s;
 
-loop:   SWAPINIT(a, es);
+loop:   if (maxdepth == 0) {
+                if (heapsort(a, n, es, cmp) == 0)
+                        return;
+        }
+        maxdepth--;
+        SWAPINIT(a, es);
         if (n < 7) {
                 for (pm = a + es; pm < a + n * es; pm += es)
                         for (pl = pm; pl > a && cmp(pl - es, pl) > 0;
@@ -101,16 +118,15 @@ loop:	SWAPINIT(a, es);
                 pl = a;
                 pn = a + (n - 1) * es;
                 if (n > 40) {
-                        d = (n / 8) * es;
-                        pl = med3(pl, pl + d, pl + 2 * d, cmp);
-                        pm = med3(pm - d, pm, pm + d, cmp);
-                        pn = med3(pn - 2 * d, pn - d, pn, cmp);
+                        s = (n / 8) * es;
+                        pl = med3(pl, pl + s, pl + 2 * s, cmp);
+                        pm = med3(pm - s, pm, pm + s, cmp);
+                        pn = med3(pn - 2 * s, pn - s, pn, cmp);
                 }
                 pm = med3(pl, pm, pn, cmp);
         }
         swap(a, pm);
         pa = pb = a + es;
-
         pc = pd = a + (n - 1) * es;
         for (;;) {
                 while (pb <= pc && (cmp_result = cmp(pb, a)) <= 0) {
@@ -149,7 +165,7 @@ loop:	SWAPINIT(a, es);
                 /* Recurse for 1st side, iterate for 2nd side. */
                 if (s > es) {
                         if (r > es)
-                                qsort(a, r / es, es, cmp);
+                                introsort(a, r / es, es, maxdepth, cmp);
                         a = pn - s;
                         n = s / es;
                         goto loop;
@@ -158,10 +174,24 @@ loop:	SWAPINIT(a, es);
                 /* Recurse for 2nd side, iterate for 1st side. */
                 if (r > es) {
                         if (s > es)
-                                qsort(pn - s, s / es, es, cmp);
+                                introsort(pn - s, s / es, es, maxdepth, cmp);
                         n = r / es;
                         goto loop;
                 }
         }
 }
+
+void
+qsort(void *a, size_t n, size_t es, int (*cmp)(const void *, const void *))
+{
+        size_t i, maxdepth = 0;
+
+        /* Approximate 2*ceil(lg(n + 1)) */
+        for (i = n; i > 0; i >>= 1)
+                maxdepth++;
+        maxdepth *= 2;
+
+        introsort(a, n, es, maxdepth, cmp);
+}
+
 DEF_STRONG(qsort);