Initial stab at wikisort

4 files changed, 227 insertions, 1 deletions

diff --git a/Makefile b/Makefile
index 6271d7c..70095f4 100644
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,5 @@
 CFLAGS  = -O3 -pipe -std=c99 -D_XOPEN_SOURCE=500
-LDFLAGS =
+LDFLAGS = -lm
 
 SRCS    = $(sort $(wildcard *.c))
 OBJS    = $(SRCS:.c=.o)
diff --git a/sorters.c b/sorters.c
index 33f85c2..99e000e 100644
--- a/sorters.c
+++ b/sorters.c
@@ -5,5 +5,6 @@ const struct sorter sorters[] = {
 	{ .name = "glibc quicksort", .func = glibc_quicksort },
 	{ .name = "glibc mergesort", .func = glibc_mergesort },
 	{ .name = "musl",            .func = musl_qsort      },
+	{ .name = "wikisort",        .func = wikisort        },
 	{ 0 }
 };
diff --git a/sorters.h b/sorters.h
index 2e1c4ec..b683e91 100644
--- a/sorters.h
+++ b/sorters.h
@@ -7,6 +7,7 @@ void musl_qsort(void *, size_t, size_t, cmpfun);
 void glibc_quicksort(void *, size_t, size_t, cmpfun);
 void glibc_mergesort(void *, size_t, size_t, cmpfun);
 void freebsd_qsort(void *, size_t, size_t, cmpfun);
+void wikisort(void *, size_t, size_t, cmpfun);
 
 extern const struct sorter {
 	const char *name;
diff --git a/wikisort.c b/wikisort.c
new file mode 100644
index 0000000..8b27c57
--- /dev/null
+++ b/wikisort.c
@@ -0,0 +1,224 @@
+#include <math.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <string.h>
+
+#include "sorters.h"
+
+struct blockarray {
+	char  *base;
+	size_t blocks;
+};
+
+static void swap(char *a, char *b, size_t width)
+{
+	while (width--) {
+		char tmp = *a;
+		*a++ = *b;
+		*b++ = tmp;
+	}
+}
+
+static void reverse(char *base, size_t nmel, size_t width)
+{
+	char *last = base + (nmel - 1) * width;
+	while (last > base) {
+		swap(base, last, width);
+		base += width;
+		last -= width;
+	}
+}
+
+/* rotates left */
+static void rotate(char *base, size_t nmel, size_t width, size_t shift)
+{
+	if (shift == 0 || shift == nmel) return;
+	reverse(base, shift, width);
+	reverse(base + shift * width, nmel - shift, width);
+	reverse(base, nmel, width);
+}
+
+/* buf must be large enough to hold a_count elements */
+static void merge_blocks(char *buf, char *base, size_t as, size_t bs, size_t width, cmpfun cmp)
+{
+	char *a = buf;
+	char *b = base + as * width;
+	char *target = base;
+
+	swap(a, base, as * width);
+
+	while (as > 0 && bs > 0) {
+		if (cmp(a, b) <= 0) {
+			swap(target, a, width);
+			a += width; as--;
+		} else {
+			swap(target, b, width);
+			b += width; bs--;
+		}
+		target += width;
+	}
+
+	swap(target, a, as * width);
+}
+
+static char *find_min_block(char *base, size_t count, size_t bwidth, size_t width, cmpfun cmp)
+{
+	char *min = base;
+	for (base += bwidth, count--; count--; base += bwidth) {
+		if (cmp(min + bwidth - width, base) > 0) min = base;
+	}
+	return min;
+}
+
+/* If a matching element exists, returns index of one of them.
+   If none exists, returns index where such an element should be inserted. */
+static size_t binary_search(const char *needle, char *haystack,
+	size_t nmel, size_t width, cmpfun cmp)
+{
+	size_t baseidx = 0;
+	while (nmel) {
+		size_t tryidx = baseidx + nmel / 2;
+		char *try = haystack + tryidx * width;
+		int sign = cmp(try, needle);
+		if (!sign) return tryidx;
+		else if (sign < 0) {
+			baseidx = tryidx + 1;
+			nmel -= nmel / 2 + 1;
+		}
+		else {
+			nmel /= 2;
+		}
+	}
+	return baseidx;
+}
+
+static size_t roll(struct blockarray *a, struct blockarray *b,
+	size_t bwidth, size_t width, char **mina, cmpfun cmp)
+{
+	size_t iters = 0;
+	char  *prevb;
+
+	do {
+		swap(a->base, b->base, bwidth);
+		if (*mina == a->base) *mina = b->base;
+		prevb = a->base;
+		a->base += bwidth;
+		b->base += bwidth;
+		b->blocks--;
+		iters++;
+	} while (b->blocks && cmp(prevb + bwidth - width, *mina) < 0);
+
+	return iters;
+}
+
+static void distribute_buffer(char *base, size_t bufnmel, size_t sortnmel, size_t width, cmpfun cmp)
+{
+	char *sorted = base + bufnmel * width;
+	while (bufnmel) {
+		char *last = sorted - width;
+		size_t equals;
+		for (equals = 1; equals < bufnmel; equals++) {
+			if (cmp(last, last - equals * width)) break;
+		}
+
+		size_t insertpos = binary_search(last, sorted, sortnmel, width, cmp);
+		rotate(sorted - equals * width, insertpos + equals, width, equals);
+
+		bufnmel  -= equals;
+		sortnmel += equals;
+		sorted   -= equals * width;
+	}
+}
+
+static const uint8_t sortnet[][2] = {
+	/* 0: index =  0 */
+	/* 1: index =  0 */
+	/* 2: index =  0 */
+	{0,1},
+	/* 3: index =  1 */
+	{0,1}, {0,2}, {1,2},
+	/* 4: index =  4 */
+	{0,1}, {2,3}, {0,2}, {1,3}, {1,2},
+	/* 5: index =  9 */
+	{0,1}, {3,4}, {2,4}, {2,3}, {1,4}, {0,3}, {0,2}, {1,3}, {1,2},
+	/* 6: index = 18 */
+	{1,2}, {4,5}, {0,2}, {3,5}, {0,1}, {3,4}, {2,5}, {0,3}, {1,4}, {2,4},
+	{1,3}, {2,3},
+	/* 7: index = 30 */
+	{1,2}, {3,4}, {5,6}, {0,2}, {3,5}, {4,6}, {0,1}, {4,5}, {2,6}, {0,4},
+	{1,5}, {0,3}, {2,5}, {1,3}, {2,4}, {2,3},
+	/* 8: index = 46 */
+	{0,1}, {2,3}, {4,5}, {6,7}, {0,2}, {1,3}, {4,6}, {5,7}, {1,2}, {5,6},
+	{0,4}, {3,7}, {1,5}, {2,6}, {1,4}, {3,6}, {2,4}, {3,5}, {3,4},
+	/* 9: index = 65, if we used a sorting network this large */
+};
+static const uint8_t sortnet_index[] = { 0, 0, 0, 1, 4, 9, 18, 30, 46, 65 };
+
+void wikisort(void *unsorted, size_t nmel, size_t width, cmpfun cmp)
+{
+	char *base = unsorted;
+
+	/* use a sorting network for small inputs */
+	if (nmel <= 8) {
+		for (int i = sortnet_index[nmel]; i < sortnet_index[nmel+1]; i++) {
+			char *elem1 = base + sortnet[i][0] * width;
+			char *elem2 = base + sortnet[i][1] * width;
+			if (cmp(elem1, elem2) > 0) swap(elem1, elem2, width);
+		}
+		return;
+	}
+
+	size_t bnmel  = sqrt(nmel);   /* elements in a block */
+	size_t bcount = nmel / bnmel;
+	size_t bwidth = bnmel * width;
+
+	/* temporary buffer is at base. carve remainder of input into a and b subarrays.
+	 * the nmels for these subarrays is specified in terms of blocks. */
+	size_t bufnmel = bnmel + nmel % bnmel;
+	struct blockarray a = { .base =   base +  width * bufnmel,  .blocks = (bcount - 1) / 2 };
+	struct blockarray b = { .base = a.base + bwidth * a.blocks, .blocks =  bcount - 1 - a.blocks };
+
+	/* sort the buffer as if it were part of the a array */
+	wikisort(  base, nmel - b.blocks * bnmel, width, cmp);
+	wikisort(b.base,        b.blocks * bnmel, width, cmp);
+
+	size_t nmel_pend_a = 0, nmel_pend_b = 0;
+	char *mina = a.base, *penda;
+	while (a.blocks) {
+		size_t rolled_blocks = b.blocks ? roll(&a, &b, bwidth, width, &mina, cmp) : 0;
+		size_t offset;
+
+		/* move mina block to immediately after rolled b blocks */
+		if (mina != a.base) swap(mina, a.base, bwidth);
+
+		/* rotate mina into correct spot in the last rolled b */
+		nmel_pend_b += rolled_blocks * bnmel;
+		if (nmel_pend_b) {
+			offset = nmel_pend_b - binary_search(a.base, a.base - nmel_pend_b * width, nmel_pend_b, width, cmp);
+			rotate(a.base - offset * width, bnmel + offset, width, offset);
+
+			/* merge any previous pending as and bs */
+			nmel_pend_b -= offset;
+			if (nmel_pend_a) {
+				merge_blocks(base, penda, nmel_pend_a, nmel_pend_b, width, cmp);
+			}
+		} else {
+			offset = 0;
+		}
+
+		nmel_pend_a = bnmel;
+		nmel_pend_b = offset;
+		penda = a.base - offset * width;
+		a.base += bwidth;
+		if (--a.blocks) {
+			mina = find_min_block(a.base, a.blocks, bwidth, width, cmp);
+		}
+	}
+
+	if (nmel_pend_a) {
+		merge_blocks(base, penda, nmel_pend_a, b.blocks * bnmel + nmel_pend_b, width, cmp);
+	}
+
+	wikisort(base, bufnmel, width, cmp);
+	distribute_buffer(base, bufnmel, nmel - bufnmel, width, cmp);
+}