Rework block merge sort algorithm

This was previously broken for certain inputs, particularly two
concatenated sorted sequences.

1 files changed, 36 insertions, 28 deletions

diff --git a/grailsort.c b/grailsort.c
index 29bd694..6f63701 100644
--- a/grailsort.c
+++ b/grailsort.c
@@ -1,38 +1,42 @@
 #include <math.h>
 #include <stddef.h>
 #include <stdint.h>
+#include <stdlib.h>
+#include <limits.h>
 
 #include "sorters.h"
 #include "common.h"
 
-#define TAIL(base,nmel,width) ((base) + ((nmel) - 1) * (width))
+static inline void *tail(char *base, size_t nmel, size_t width)
+{
+	return base + (nmel-1) * width;
+}
 
 /* get index of last block whose head is less than the previous block's tail */
 static size_t last_overlap(char *base, size_t bcount, size_t bwidth, size_t width, cmpfun cmp)
 {
-	for (char *cur = TAIL(base, bcount, bwidth); --bcount; cur -= bwidth)
+	for (char *cur = tail(base, bcount, bwidth); --bcount; cur -= bwidth)
 		if (cmp(cur - width, cur) > 0) break;
 	return bcount;
 }
 
-size_t merge(char *base, size_t anmel, size_t bnmel, size_t bufnmel, size_t width, cmpfun cmp)
+size_t merge(char *buf, char *a, char *b, size_t bufnmel, size_t anmel, size_t bnmel, size_t width, cmpfun cmp)
 {
-	char *a_end   = base  +   anmel * width;
-	char *b_end   = a_end +   bnmel * width;
-	char *buf_end = b_end + bufnmel * width;
-	size_t count = 0;
-
 	while (bnmel) {
-		char *a   =   a_end - width;
-		char *b   =   b_end - width;
-		char *buf = buf_end - width;
-
-		if (cmp(a, b) > 0) { swap(a, buf, width); a_end = a; anmel--; }
-		else               { swap(b, buf, width); b_end = b; bnmel--; }
-		count++;
-		buf_end = buf;
+		char *a_tail   = tail(a,   anmel,   width);
+		char *b_tail   = tail(b,   bnmel,   width);
+		char *buf_tail = tail(buf, bufnmel, width);
+
+		if (cmp(a_tail, b_tail) <= 0) {
+			swap(buf_tail, b_tail, width);
+			bnmel--;
+		} else {
+			swap(buf_tail, a_tail, width);
+			anmel--;
+		}
+		bufnmel--;
 	}
-	return count;
+	return anmel;
 }
 
 void grailsort(void *unsorted, size_t nmel, size_t width, cmpfun cmp)
@@ -51,8 +55,10 @@ void grailsort(void *unsorted, size_t nmel, size_t width, cmpfun cmp)
 	size_t acount  = blocks / 2;
 	size_t bcount  = blocks - acount;
 
-	char *a = base;
-	char *b = a + acount * bwidth;
+	char *a   = base;
+	char *b   = a + acount * bwidth;
+	char *buf = b + bcount * bwidth;
+
 	grailsort(a, acount * blknmel, width, cmp);
 	grailsort(b, bcount * blknmel, width, cmp);
 
@@ -60,20 +66,22 @@ void grailsort(void *unsorted, size_t nmel, size_t width, cmpfun cmp)
 	grailsort(base, blocks, bwidth, cmp);
 
 	/* merge, starting from the end and working towards the beginning */
-	size_t bufidx = blocks * blknmel;
+	size_t last_nmel = blknmel;
 	while (blocks > 1) {
 		size_t overlap = last_overlap(base, blocks, bwidth, width, cmp);
-		if (overlap) {
-			bufidx -= merge(TAIL(base, overlap, bwidth), blknmel, bufidx - overlap*blknmel, bufnmel, width, cmp);
-		}
-		blocks = overlap;
-	}
 
-	if (bufidx) {
-		rotate(base, (bufidx + bufnmel) * width, bufidx * width);
+		if (!overlap) break;
+
+		char *block_a = base + bwidth * (overlap-1);
+		char *block_b = base + bwidth *  overlap;
+		size_t bnmel  = (blocks - overlap - 1) * blknmel + last_nmel;
+
+		last_nmel = merge(buf, block_a, block_b, bufnmel, blknmel, bnmel, width, cmp);
+		buf       = block_a + last_nmel * width;
+		blocks    = overlap;
 	}
 
-distribute:
+	rotate(base, buf-base + bufnmel*width, buf-base);
 	grailsort(base, bufnmel, width, cmp);
 	distribute_buffer(base, bufnmel, nmel - bufnmel, width, cmp);
 }