diff --git a/thiersort3.h b/thiersort3.h new file mode 100644 index 0000000..ecd039e --- /dev/null +++ b/thiersort3.h @@ -0,0 +1,109 @@ +#ifndef THIER_SORT3_H +#define THIER_SORT3_H +#include +#include "threepass_xbit.h" +/* A non-implace tricky float-hackz based bucket sort variant. Uses threepass_xbit and removes some compies! */ + +/* 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 + 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 thiersort3(uint32_t *arr, uint32_t *temparr, int n) { + 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_bucket3(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 */ + /* 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_bucket3(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; + } + + /* 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); + } +} + +#endif /* THIER_SORT3_H */ diff --git a/threepass_xbit.h b/threepass_xbit.h new file mode 100644 index 0000000..4298664 --- /dev/null +++ b/threepass_xbit.h @@ -0,0 +1,136 @@ +#ifndef THREE_PASS_XB_H +#define THREE_PASS_XB_H +/* Basically like threepass.h but but fewer bits - only for thiersort3! */ + +/* How the 28-30 bits gets separated? This is chosen because I saw that: */ +/* - in float trickery I saw at around 4 billion the buckets change around 200 million */ +/* - Which is smaller than 1/2 or 1/4 or 1/8 or 1/16 of it because 16*200 mil = 3.2 bill 28 bit is needed */ +/* - This means those bits in numbers DO get used when we are running as internal sort of thier3 */ +/* - But the numbers high bits (29..32th bits) stay the same inside each bucket so we spare it! */ +/* - If you know your data does not have all the range of 32 bits you can #define these by #define CUSTOM_TPBX_BITS */ +#ifndef CUSTOM_TPBX_BITS +#define TPBX1 10 // top +#define TPBX2 9 // mid +#define TPBX3 9 // bottom +#endif /* CUSTOM_TPBX_BITS */ + +static inline constexpr uint32_t min3u32_xb(uint32_t a, uint32_t b, uint32_t c) noexcept { + return (a <= b) ? + ((a <= c) ? a : c) : + ((b <= c) ? b : c); +} + +/** + * Simple three-pass (ok: 3 + 1) bottom-up internal radix sort writter for thiersort3 + * + * @param a The array to sort - will be changed too! + * @param buf Result array with the same size - result will be here + * @param n The number of elements + */ +static inline void threepass_xb(uint32_t *a, uint32_t *buf, int n) noexcept { + assert(buf != NULL); + constexpr int shr1 = TPBX3 + TPBX2; + constexpr int shr2 = TPBX3; + constexpr int shr3 = 0; + constexpr int mask1 = (1 << TPBX1) - 1; + constexpr int mask2 = (1 << TPBX2) - 1; + constexpr int mask3 = (1 << TPBX3) - 1; + + /* helper buffers. */ + int sz = n * sizeof(a[0]); + + static thread_local uint32_t bucket1[1 << TPBX1]; + memset(bucket1, 0, (1 << TPBX1) * sizeof(uint32_t)); + static thread_local uint32_t bucket2[1 << TPBX2]; + memset(bucket2, 0, (1 << TPBX2) * sizeof(uint32_t)); + static thread_local uint32_t bucket3[1 << TPBX3]; + memset(bucket3, 0, (1 << TPBX3) * sizeof(uint32_t)); + + /* Count occurences (can count together with good ILP) */ + #pragma GCC unroll 64 + for(uint32_t i = 0; i < n; ++i) { + ++bucket1[(a[i] >> shr1) & mask1]; + ++bucket2[(a[i] >> shr2) & mask2]; + ++bucket3[(a[i] >> shr3) & mask3]; + } + + /* Count prefix sums - try as much ILP as possible because bigger arrays than usual! */ + uint32_t prev1 = 0; + uint32_t prev2 = 0; + uint32_t prev3 = 0; + uint32_t common = min3u32_xb( + (1 << TPBX1), + (1 << TPBX2), + (1 << TPBX3) + ); + int i = 0; + #pragma GCC unroll 8 + for (; i < common; ++i) { + bucket1[i] += prev1; + prev1 = bucket1[i]; + bucket2[i] += prev2; + prev2 = bucket2[i]; + bucket3[i] += prev3; + prev3 = bucket3[i]; + } + /* Do remaining 1 */ + for (int j = i; j < (1 << TPBX1); ++j) { + bucket1[j] += prev1; + prev1 = bucket1[j]; + } + /* Do remaining 2 */ + for (int j = i; j< (1 << TPBX2); ++j) { + bucket2[j] += prev2; + prev2 = bucket2[j]; + } + /* Do remaining 3 */ + for (int j = i; j < (1 << TPBX3); ++j) { + bucket3[j] += prev3; + prev3 = bucket3[j]; + } + + // Bottom digit a->buf + // right-to-left to ensure already sorted digits order we keep for iterations + #pragma GCC unroll 48 + for(uint32_t i = n; i > 0; --i) { + // Prefetch caches + //__builtin_prefetch(&a[i-8]); + // Get num and its new offset / location + auto num = a[i - 1]; + auto bkeyni = (num >> shr3) & mask3; + auto offset = --bucket3[bkeyni]; + + // Add to the proper target location + buf[offset] = num; + } + // Mid digit buf->a + // right-to-left to ensure already sorted digits order we keep for iterations + #pragma GCC unroll 48 + for(uint32_t i = n; i > 0; --i) { + // Prefetch caches + //__builtin_prefetch(&buf[i-8]); + // Get num and its new offset / location + auto num = buf[i - 1]; + auto bkeyni = (num >> shr2) & mask2; + auto offset = --bucket2[bkeyni]; + + // Add to the proper target location + a[offset] = num; + } + // Top digit a->buf + // right-to-left to ensure already sorted digits order we keep for iterations + #pragma GCC unroll 48 + for(uint32_t i = n; i > 0; --i) { + // Prefetch caches + // __builtin_prefetch(&a[i-16]); + // Get num and its new offset / location + auto num = a[i - 1]; + auto bkeyni = (num >> shr1) & mask1; + auto offset = --bucket1[bkeyni]; + + // Add to the proper target location + buf[offset] = num; + } +} + +#endif /* THREE_PASS_XB_H */