magyarsort/thiersort2.h

#ifndef THIER_SORT2_H
#define THIER_SORT2_H
#include <stdint.h>
#include "qsort/schwab_sort.h"
#include "threepass.h"
/* A non-implace tricky float-hackz based bucket sort variant. Uses schwabsort or threepass! */

#ifdef _MSC_VER
#define KM_PREFETCH(x)
#define LIKELY(x)
#define UNLIKELY(x)
#define KM_NOINLINE __declspec(noinline)
#define KM_ALWAYS_INLINE __forceinline
#else
#define KM_PREFETCH(x) __builtin_prefetch(x)
#define LIKELY(x) __builtin_expect((x),1)
/* alternative: #define LIKELY(x) __builtin_expect(((x) != 0),1) */
#define UNLIKELY(x) __builtin_expect((x),0)
#define KM_NOINLINE __attribute__ ((noinline))
#define KM_ALWAYS_INLINE __attribute__ ((always_inline))
#endif

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

/** Tells from the key which bucket it is in. */
static inline uint32_t witch_bucket(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 */
	th2_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
	return witch_base * 128 + ((as.u >> (23 - 7)) & 127); // 0..2047
	/* Alternative (but I measure it being worse):
	return (as.u >> 23);
	*/
}

/**
 * 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 thiersort2(uint32_t *arr, uint32_t *temparr, int n, sch_rand_state *rstate) {
	int bucket[4096]; /* Inclusive */
	int bucket_end[4096]; /* Not inclusive */

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

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

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

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

	/* arr -> temparr */
	/* Move to the buckets - backwards going save a few cache miss */
	/* Rem.: This also changes bucket[i] so they will point to bucket beginnings */
	#pragma GCC unroll 64
	for(int i = 0; i < n; ++i) {
		uint32_t num = arr[i];
		uint32_t witch = witch_bucket(num);
		int offset = (--bucket[witch]);
		temparr[offset] = num;
	}

	/* temparr -> arr each bucket and sort them in-place */
	#pragma GCC unroll 64
	for(int b = 0; b < 4096; ++b) {
		int begin = bucket[b];
		int end = bucket_end[b];
		
		/* Ensure exists */
		if(begin >= end) {
			continue;
		}

		/* We make a three-way FAST quicksort partitioning with first elem pivot: */
		/* Basically a Lomuto-like unidirectional partitioning for pivot and two-way for small/big */
		int smalli = begin;
		int biggi = end - 1; /* always exists */
		int i = begin;
		uint32_t pivot = temparr[i];
		#pragma GCC unroll 4
		for(int j = begin + 1; j < end; ++j) {
			if(LIKELY(temparr[j] != pivot)) {
				/* Branchless partitioning */
				/* copy to left */
				/* Rem.: Because we overwrite by copy we can simplify this line */
				/* arr[smalli] = (temparr[j] < pivot) ? temparr[j] : arr[smalli]; */
				arr[smalli] = temparr[j];
				smalli += (temparr[j] < pivot);
				/* copy to right */
				/* Rem.: Because we overwrite by copy we can simplify this line */
				/* arr[biggi] = (temparr[j] > pivot) ? temparr[j] : arr[biggi]; */
				arr[biggi] = temparr[j];
				biggi -= (temparr[j] > pivot);
			} else {
				/* swap to front partition */
				++i;
				uint32_t tmp = temparr[i];
				temparr[i] = temparr[j];
				temparr[j] = tmp;
			}
		}
		/* Copy the mid elements back */
		int target = smalli;
		#pragma GCC unroll 64
		for(int j = begin; j < i + 1; ++j) {
			arr[target++] = temparr[j];
		}

		/* Call schwabsort - only to [begin..smalli) and (biggie..end) */
#ifdef USE_SCHWAB
		schwab_sort(arr, begin, smalli - 1, rstate);
		schwab_sort(arr, biggi + 1, end - 1, rstate);
#else
		threepass(&arr[begin], smalli - begin);
		threepass(&arr[biggi + 1], end - (biggi + 1));
#endif /* USE_SCHWAB */
	}
}

#endif /* THIER_SORT2_H */