magyarsort/thiersort3.h

#ifndef THIER_SORT3_H
#define THIER_SORT3_H
#include <stdint.h>
/* A non-implace tricky float-hackz based bucket sort variant. Uses threepass_xbit and removes some compies! */

#ifndef NO_MLOCK
#include <sys/mman.h> // mlock & munlock
#endif // !NO_MLOCK

/* Disables extra 1-bit split processing before threepass_xb */
#define NO_EXTRA_BIT /* I measure extra split to be slightly slower because overhead - despite less cache misses */

#ifndef NO_EXTRA_BIT
#define CUSTOM_TPBX_BITS
#define TPBX1 9 /* top - not 10 because our extra bit split / partitioning trickery! */
#define TPBX2 9 /* mid */
#define TPBX3 9 /* bottom */
#endif /* NO_EXTRA_BIT */
#include "threepass_xbit.h"

/* Float and unsigned32 reinterpreter */
union th3_fu {
	float f;
	uint32_t u;
};
typedef union th3_fu th3_fu;

/** Tells from the key which bucket it is in. */
static inline uint32_t witch_bucket3(uint32_t key) {
	/* This is hackz to misuse int->float converter HEAVILY and IEE for bucketing */
	/* https://en.wikipedia.org/wiki/Single-precision_floating-point_format */

	/* Old Hungarian ASM trick I know from Tomcat/Abaddon mailing list and prog.hu */
	/* "A következő nagyon fontos gondolat, hogy két lebegőpontos számot, */
	/* amennyiben pozitív, "simán" is összehasonlíthatunk" */
	/* See: https://prog.hu/cikkek/100239/fpu-gems */

	/* This approach uses 12 bits from a 32 bit float to map onto a byte bucket index */
	th3_fu as;
	as.f = (float) key;

	uint32_t witch_base = (key <= 2) ? 0 : (as.u >> 23) - 128; // 0, [127..159] -> [0..31]
	return witch_base * 8 + ((as.u >> (23 - 3)) & 7); // 0..255
	/* All below ones were worse, because three-ways 28 bit sorting would just make unnecessary work at cost of some locality */
	//return witch_base * 128 + ((as.u >> (23 - 7)) & 127); // 0..4095
	//return witch_base * 64 + ((as.u >> (23 - 6)) & 63); // 0..2047
	//return witch_base * 32 + ((as.u >> (23 - 5)) & 31); // 0..1023
	//return witch_base; // not enough!
}

/**
 * Sort the array using the temporary array of the same size with fast bucket sort thiersort.
 *
 * @param arr The array to sort, will contain result afterwards
 * @param temparr The temporary array with same size
 * @param n Number of elements in arr and temparr
 * @param rstate Create with sch_rand_state rstate = schwab_rand_state(junk_uint32_t);
 */
static inline void thiersort3(uint32_t *arr, uint32_t *temparr, int n) {
	int bucket[256]; /* Inclusive */
	int bucket_end[256]; /* Not inclusive */

#ifndef NO_EXTRA_BIT
	int bucket_leftend[256]; /* for extra 1bit split processing */
	int bucket_left[256]; /* for extra 1bit split processing */
#endif /* NO_EXTRA_BIT */

	/* Check if need to sort at all - needed for invariants later */
	if(n < 2) {
		return;
	}

#ifndef NO_MLOCK
	mlock(arr, n * sizeof(uint32_t));
	mlock(temparr, n * sizeof(uint32_t));
#endif // !NO_MLOCK

	/* Count */
	#pragma GCC unroll 64
	for(int i = 0; i < 256; ++i) {
		bucket[i] = 0;
	}
	#pragma GCC unroll 128
	for(int i = 0; i < n; ++i) {
		++bucket[witch_bucket3(arr[i])];
	}

	/* Prefix sum (like in Magyarsort) */
	uint32_t prev = 0;
	#pragma GCC unroll 4
	for (int i = 0; i < 256; i++) {
		bucket[i] += prev;
		prev = bucket[i];
	}

	/* Save end-offsets */
	#pragma GCC unroll 64
	for(int i = 0; i < 256; ++i) {
		bucket_end[i] = bucket[i];
	}

#ifndef NO_EXTRA_BIT
	/* Save left-offsets */
	bucket_left[0] = 0;
	bucket_leftend[0] = 0;
	#pragma GCC unroll 64
	for(int i = 0; i < 255; ++i) {
		bucket_left[1 + i] = bucket[i];
		bucket_leftend[1 + i] = bucket[i];
	}
#endif /* NO_EXTRA_BIT */

	/* arr -> temparr */
#ifdef NO_EXTRA_BIT
	/* Move to the buckets */
	/* Rem.: This also changes bucket[i] so they will point to bucket beginnings */
	#pragma GCC unroll 128
	for(int i = 0; i < n; ++i) {
		uint32_t num = arr[i];
		uint32_t witch = witch_bucket3(num);
		int offset = (--bucket[witch]);
		temparr[offset] = num;
	}

	/* temparr -> arr each bucket and sort them in-place */
	#pragma GCC unroll 2
	for(int b = 0; b < 256; ++b) {
		int begin = bucket[b];
		int end = bucket_end[b];

		/* Ensure exists */
		if(begin >= end) {
			continue;
		}

		/* Use our specialized threepass */
		/* - which does not use up all the bits */
		/* - which does not allocate and copy back to arr here */
		uint32_t n = end - begin;
		uint32_t *a = &(temparr[begin]);
		uint32_t *buf = &(arr[begin]);
		threepass_xb(a, buf, n);
	}
#else
	/* Move to the buckets by split-partitioning on the 28th bit */
	/* Rem.: This also changes bucket[i] so they will point to bucket.right beginnings */
	/* Rem.: This also changes bucket_leftend[i] so they will point to bucket.left endings (needed in-process only) */
	#pragma GCC unroll 64
	for(int i = 0; i < n; ++i) {
		uint32_t num = arr[i];
		uint32_t witch = witch_bucket3(num);
		int offset = (num & (1 << 27)) ?
			(--bucket[witch]) :
			(bucket_leftend[witch]++);

		temparr[offset] = num;
	}

	/* temparr -> arr each bucket and sort them in-place */
	#pragma GCC unroll 2
	for(int b = 0; b < 256; ++b) {
		assert(bucket_leftend[b] == bucket[b]);
		int lbegin = bucket_left[b];
		int lend = bucket[b]; /* non-inclusive */
		int rbegin = bucket[b];
		int rend = bucket_end[b]; /* non-inclusive */

		/* Ensure exists and process left part */
		if(lbegin < lend) {
			/* Use our specialized threepass */
			/* - which does not use up all the bits */
			/* - which does not allocate and copy back to arr here */
			uint32_t n = lend - lbegin;
			uint32_t *a = &(temparr[lbegin]);
			uint32_t *buf = &(arr[lbegin]);
			threepass_xb(a, buf, n);
		}

		/* Ensure exists and process right part */
		if(rbegin < rend) {
			/* Use our specialized threepass */
			/* - which does not use up all the bits */
			/* - which does not allocate and copy back to arr here */
			uint32_t n = rend - rbegin;
			uint32_t *a = &(temparr[rbegin]);
			uint32_t *buf = &(arr[rbegin]);
			threepass_xb(a, buf, n);
		}
	}
#endif /* NO_EXTRA_BIT */
#ifndef NO_MLOCK
	munlock(arr, n * sizeof(uint32_t));
	munlock(temparr, n * sizeof(uint32_t));
#endif // !NO_MLOCK
}

#endif /* THIER_SORT3_H */