202 lines
6.4 KiB
C
202 lines
6.4 KiB
C
#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 */
|
|
|
|
/* Only define these to be smaller than these in case your input data shows you can - never change orders or make bigger! */
|
|
#ifndef CUSTOM_TPBX_BITS
|
|
#define TPBX1 10 // top
|
|
#define TPBX2 9 // mid
|
|
#define TPBX3 9 // bottom
|
|
#endif /* CUSTOM_TPBX_BITS */
|
|
|
|
|
|
#ifndef TPBX_NO_DEF_WCACHE
|
|
#define TPBX_WCACHE 2 /* depth of the write-cache - I let the user */
|
|
#endif /* TPBX_NO_DEF_WCACHE */
|
|
|
|
#ifdef __cplusplus
|
|
#define TPBX_NOEXCEPT noexcept
|
|
#define TPBX_CONSTEXPR constexpr
|
|
#else
|
|
#define TPBX_NOEXCEPT
|
|
#define TPBX_CONSTEXPR const
|
|
#endif /* __cplusplus */
|
|
|
|
static inline TPBX_CONSTEXPR uint32_t min3u32_xb(uint32_t a, uint32_t b, uint32_t c) TPBX_NOEXCEPT {
|
|
return (a <= b) ?
|
|
((a <= c) ? a : c) :
|
|
((b <= c) ? b : c);
|
|
}
|
|
|
|
static inline TPBX_CONSTEXPR uint32_t max3u32_xb(uint32_t a, uint32_t b, uint32_t c) TPBX_NOEXCEPT {
|
|
return (a >= b) ?
|
|
((a >= c) ? a : c) :
|
|
((b >= c) ? b : c);
|
|
}
|
|
|
|
/** Copy the elements to their respective radics-place (f->t copy) bit should be the TPBX1,2,3 */
|
|
static inline void copy_radics_tpxp(
|
|
uint32_t *f,
|
|
uint32_t *t,
|
|
uint32_t *bucket,
|
|
uint32_t shr,
|
|
uint32_t mask,
|
|
uint32_t n,
|
|
uint32_t bit) {
|
|
#ifndef TPBX_WCACHE
|
|
// 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
|
|
uint32_t num = f[i - 1];
|
|
uint32_t bkeyni = (num >> shr) & mask;
|
|
uint32_t offset = --bucket[bkeyni];
|
|
|
|
// Add to the proper target location
|
|
t[offset] = num;
|
|
}
|
|
#else
|
|
int stride = (1 << bit);
|
|
static __thread uint32_t write_cache[1024 * TPBX_WCACHE];
|
|
static __thread uint8_t write_qc[1024];
|
|
|
|
for(uint32_t i = 0; i < 1024; ++i) {
|
|
write_qc[i] = 0;
|
|
}
|
|
|
|
/* right-to-left to ensure already sorted digits order we keep for iterations */
|
|
#pragma GCC unroll 4
|
|
for(uint32_t i = n; i > 0; --i) {
|
|
/* calculate which bucket this would belong and where the queue is */
|
|
uint32_t num = f[i - 1];
|
|
uint32_t bkeyni = (num >> shr) & mask;
|
|
uint8_t qc = write_qc[bkeyni];
|
|
|
|
if(qc < TPBX_WCACHE) {
|
|
/* Cache is not full yet */
|
|
write_cache[bkeyni + stride * qc] = num;
|
|
++write_qc[bkeyni];
|
|
} else {
|
|
/* Cache is full, make a real writeout - all from this queue */
|
|
|
|
/* We also directly move the current element too - micro optimized */
|
|
uint32_t offset = (bucket[bkeyni] -= (TPBX_WCACHE + 1));
|
|
t[offset++] = num;
|
|
/* Copy from the cache - order counts here because keeping opposites! */
|
|
for(uint32_t j = 0; j < TPBX_WCACHE; ++j) {
|
|
uint32_t cached = write_cache[bkeyni + stride * (TPBX_WCACHE - 1 - j)];
|
|
t[offset + j] = cached;
|
|
}
|
|
|
|
/* Reset this queue */
|
|
write_qc[bkeyni] = 0;
|
|
}
|
|
}
|
|
|
|
/* Write out any data remaining in the write cache */
|
|
for(uint32_t bkeyni = 0; bkeyni < stride; ++bkeyni) {
|
|
uint32_t qc = write_qc[bkeyni];
|
|
uint32_t offset = (bucket[bkeyni] -= qc);
|
|
if(qc > 0) for(uint32_t j = 0; j < qc; ++j) {
|
|
uint32_t cached = write_cache[bkeyni + stride * (qc - 1 - j)];
|
|
t[offset + j] = cached;
|
|
}
|
|
}
|
|
#endif /* TPBX_WCACHE */
|
|
}
|
|
|
|
/* I pulled these out only for better flame graph support */
|
|
/** Count occurences (can count together with good ILP) */
|
|
static inline void count_occurences_tpxp(uint32_t *bucket1, uint32_t *bucket2, uint32_t *bucket3, const uint32_t shr1, const uint32_t shr2, const uint32_t shr3, const uint32_t mask1, const uint32_t mask2, const uint32_t mask3, uint32_t *a, uint32_t n) TPBX_NOEXCEPT {
|
|
#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];
|
|
}
|
|
}
|
|
|
|
/**
|
|
* 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, uint32_t n) TPBX_NOEXCEPT {
|
|
assert(buf != NULL);
|
|
TPBX_CONSTEXPR uint32_t shr1 = TPBX3 + TPBX2;
|
|
TPBX_CONSTEXPR uint32_t shr2 = TPBX3;
|
|
TPBX_CONSTEXPR uint32_t shr3 = 0;
|
|
TPBX_CONSTEXPR uint32_t mask1 = (1 << TPBX1) - 1;
|
|
TPBX_CONSTEXPR uint32_t mask2 = (1 << TPBX2) - 1;
|
|
TPBX_CONSTEXPR uint32_t mask3 = (1 << TPBX3) - 1;
|
|
|
|
/* helper buffers. */
|
|
uint32_t sz = n * sizeof(a[0]);
|
|
|
|
static __thread uint32_t bucket1[1 << TPBX1];
|
|
memset(bucket1, 0, (1 << TPBX1) * sizeof(uint32_t));
|
|
static __thread uint32_t bucket2[1 << TPBX2];
|
|
memset(bucket2, 0, (1 << TPBX2) * sizeof(uint32_t));
|
|
static __thread uint32_t bucket3[1 << TPBX3];
|
|
memset(bucket3, 0, (1 << TPBX3) * sizeof(uint32_t));
|
|
|
|
count_occurences_tpxp(bucket1, bucket2, bucket3, shr1, shr2, shr3, mask1, mask2, mask3, a, n);
|
|
|
|
/* 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)
|
|
);
|
|
uint32_t 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 (uint32_t j = i; j < (1 << TPBX1); ++j) {
|
|
bucket1[j] += prev1;
|
|
prev1 = bucket1[j];
|
|
}
|
|
/* Do remaining 2 */
|
|
for (uint32_t j = i; j< (1 << TPBX2); ++j) {
|
|
bucket2[j] += prev2;
|
|
prev2 = bucket2[j];
|
|
}
|
|
/* Do remaining 3 */
|
|
for (uint32_t j = i; j < (1 << TPBX3); ++j) {
|
|
bucket3[j] += prev3;
|
|
prev3 = bucket3[j];
|
|
}
|
|
|
|
// Bottom digit a->buf
|
|
copy_radics_tpxp(a, buf, bucket3, shr3, mask3, n, TPBX3);
|
|
|
|
// Mid digit buf->a
|
|
copy_radics_tpxp(buf, a, bucket2, shr2, mask2, n, TPBX2);
|
|
|
|
// Top digit a->buf
|
|
copy_radics_tpxp(a, buf, bucket1, shr1, mask1, n, TPBX1);
|
|
}
|
|
|
|
#endif /* THREE_PASS_XB_H */
|