Index: sys/libkern/qsort.c
===================================================================
--- sys/libkern/qsort.c
+++ sys/libkern/qsort.c
@@ -1,6 +1,7 @@
 /*-
  * Copyright (c) 1992, 1993
  *	The Regents of the University of California.  All rights reserved.
+ * Copyright (c) 2016 Hans Petter Selasky. All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
@@ -41,10 +42,15 @@
 static __inline char	*med3(char *, char *, char *, cmp_t *, void *);
 static __inline void	 swapfunc(char *, char *, int, int);
 
-#define min(a, b)	(a) < (b) ? (a) : (b)
+#define min(a, b)	((a) < (b) ? (a) : (b))
 
 /*
  * Qsort routine from Bentley & McIlroy's "Engineering a Sort Function".
+ *
+ * NOTE: This implementation of qsort() was designed to not have the
+ * worst case complexity of N**2, as seen with the regular quick sort
+ * functions. The worst case complexity of this sorting function is
+ * log2(N) * log2(N) * N.
  */
 #define swapcode(TYPE, parmi, parmj, n) { 		\
 	long i = (n) / sizeof (TYPE); 			\
@@ -97,6 +103,74 @@
               :(CMP(thunk, b, c) > 0 ? b : (CMP(thunk, a, c) < 0 ? a : c ));
 }
 
+/*
+ * The following sorting algorithm has a worst case complexity of
+ * log2(N) * log2(N) * N and does not require any working memory to
+ * sort the input. It replaces the insertion sort algorithm which has
+ * a worst case complexity of (N * N). This alorithm requires that N
+ * is a powers of two.
+ */
+void
+#ifdef I_AM_QSORT_R
+hpsort_r(void *a, size_t n, size_t es, void *thunk, cmp_t *cmp)
+#define HPSORT(a,n,es,thunk,cmp) hpsort_r(a,n,es,thunk,cmp)
+#else
+hpsort(void *a, size_t n, size_t es, cmp_t *cmp)
+#define HPSORT(a,n,es,thunk,cmp) hpsort(a,n,es,cmp)
+#endif
+{
+	char *pa;
+	char *pb;
+	size_t x;
+	size_t y;
+	size_t z;
+	int swaptype;
+	int tog;
+
+	if (n <= 1)
+		return;
+
+	SWAPINIT(a, es);
+
+	while (1) {
+		/* detect when array is fully sorted */
+		pa = (char *)a + es;
+		pb = (char *)a + n * es;
+		for (; pa != pb; pa += es) {
+			if (CMP(thunk, pa, pa - es) < 0)
+				break;
+		}
+		if (pa == pb)
+			break;
+
+		/* sort array in toggle order */
+		for (tog = 1, x = 1; x != n; x *= 2) {
+			for (y = 0; y != n; y += (2 * x), tog = -tog) {
+				pa = (char *)a + y * es;
+				pb = pa + x * es;
+				for (z = 0; z != x; z++, pa += es, pb += es) {
+					if ((tog * CMP(thunk, pa, pb)) > 0) {
+						swap(pa, pb);
+					}
+				}
+			}
+		}
+
+		/* sort array in rising order */
+		for (x = 1; x != n; x *= 2) {
+			for (y = 0; y != n; y += (2 * x)) {
+				pa = (char *)a + y * es;
+				pb = pa + x * es;
+				for (z = 0; z != x; z++, pa += es, pb += es) {
+					if (CMP(thunk, pa, pb) > 0) {
+						swap(pa, pb);
+					}
+				}
+			}
+		}
+	}
+}
+
 #ifdef I_AM_QSORT_R
 void
 qsort_r(void *a, size_t n, size_t es, void *thunk, cmp_t *cmp)
@@ -158,13 +232,65 @@
 		pb += es;
 		pc -= es;
 	}
-	if (swap_cnt == 0) {  /* Switch to insertion sort */
+	if (swap_cnt == 0 && n < 16) {   /* Switch to insertion sort */
 		for (pm = (char *)a + es; pm < (char *)a + n * es; pm += es)
 			for (pl = pm; pl > (char *)a && CMP(thunk, pl - es, pl) > 0;
 			     pl -= es)
 				swap(pl, pl - es);
 		return;
 	}
+	if (swap_cnt == 0) {  /* Fallback to hpsort */
+		size_t msb = ((size_t)1) << (flsl(n) - 1);
+		size_t p;
+		size_t x;
+
+		/* sort upper partition */
+		HPSORT((char *)a + (n - msb) * es, msb, es, thunk, cmp);
+
+		/* check for lower partition */
+		if (n == msb)
+			return;
+
+		/* compute size of lower partition */
+		p = ((size_t)1) << (flsl(n ^ msb) - 1);
+
+		/* round up */
+		if (n & (p - 1))
+			p *= 2;
+
+		/* sort lower partition */
+		HPSORT(a, p, es, thunk, cmp);
+
+		/* merge the two sorted partitions */
+		pa = (char *)a;
+		pb = (char *)a + (p * es);
+		pc = (char *)a + (p * es);
+		pd = (char *)a + (n * es);
+		for (x = p; x--; ) {
+			if (pa == pb)
+				break;
+			if (pc == pd)
+				break;
+			if (CMP(thunk, pa, pc) < 0)
+				pa += es;
+			else
+				pc += es;
+		}
+
+		/* figure out the smallest partition */
+		x = min(pc - pb, pb - pa);
+
+		pa = (char *)a + (p * es) - x;
+		pb = (char *)a + (p * es);
+
+		vecswap(pa, pb, x);
+		HPSORT(a, p, es, thunk, cmp);
+
+		/* loop on biggest remainder */
+		a = (char *)a + (n - msb) * es;
+		n = msb;
+		goto loop;
+	}
 
 	pn = (char *)a + n * es;
 	r = min(pa - (char *)a, pb - pa);
Index: sys/sys/libkern.h
===================================================================
--- sys/sys/libkern.h
+++ sys/sys/libkern.h
@@ -125,6 +125,10 @@
 void	*memcchr(const void *s, int c, size_t n);
 int	 memcmp(const void *b1, const void *b2, size_t len);
 void	*memmem(const void *l, size_t l_len, const void *s, size_t s_len);
+void	 hpsort(void *base, size_t nmemb, size_t size,
+	    int (*compar)(const void *, const void *));
+void	 hpsort_r(void *base, size_t nmemb, size_t size, void *thunk,
+	    int (*compar)(void *, const void *, const void *));
 void	 qsort(void *base, size_t nmemb, size_t size,
 	    int (*compar)(const void *, const void *));
 void	 qsort_r(void *base, size_t nmemb, size_t size, void *thunk,