在 uint64 中对 38 位进行基数排序？
Posted
技术标签:
【中文标题】在 uint64 中对 38 位进行基数排序？【英文标题】：Radix sort on 38 bits in a uint64? 【发布时间】：2020-01-28 17:34:46 【问题描述】：问题：是否有“Dial_A_Bit”基数排序可用于对数据类型中的位子集进行排序？
uint64 *Radix_Sort_Dial_A_Bit(uint64_t *a, int num_a, int sort_bits);
具体来说，是否可以对 64 位数据进行 38 位排序，并且其速度介于如下所示的 32/64 和 48/64 之间？
uint64 *Radix_Sort_ui64_38MSB(uint64_t *a, int num_a);
请注意，与对 uint64_t[] 中的所有 64 位进行排序相比，对 48 位和 32 位排序的研究如何验证了速度和正确性的提高。
似乎对数据包大小的子集进行排序的 Radix_Sort() 通常有用且高效，只对需要的内容进行排序。
在某些情况下，结果是针对每个 Pixel 计算并需要排序的。一个 uint64_t[] 用于保存计算结果和 XY 位置。
总共需要 26 位（X 为 13，Y 为 13，最大分辨率为 8192）来保存像素 XY 坐标，剩下 38 位用于数据排序。
可以使用 Radix_Sort_Uint64() 对整个 64 位包进行排序。
一种更快的方法是使用 Radix_Sort_Uint48()（见下文），因此在排序中不考虑最后 16 位。这对所有数据进行了正确排序，并对 13 个 X 坐标位中的 10 个进行了排序，这是不需要的。
由于性能几乎与排序的位成线性比例，因此在最佳情况下，排序中只考虑 38 个最高有效位。
即使是 40 位基数排序也比使用 48 位要好。我试图将工作的 48 位基数排序概括为在 40 位上运行，但它不能正确排序。
QSort_uint64_38_msb()：
static inline int comp_uint64_38_msb(const void *a, const void *b)
register int64_t ca, cb;
ca = (*((uint64_t *)a)) >> 26; // Send 26 LSBs to oblivion
cb = (*((uint64_t *)b)) >> 26; // Send 26 LSBs to oblivion
return((ca > cb) - (ca < cb)); // Calcs to [+1, 0 -1]
如您所见，48 位排序比完整的 64 位排序要快得多。 32 位排序几乎是完整 64 位排序的两倍。而 qsort 远远落后：
Time= 2.136 sec = 3.390%, RADIX_SORT_FFFF0000 , hits=4, 0.534 sec each
Time= 2.944 sec = 4.672%, RADIX_SORT_FFFFFF00 , hits=4, 0.736 sec each
Time= 4.691 sec = 7.444%, RADIX_SORT_64 , hits=5, 0.938 sec each
Time= 25.209 sec = 40.010%, QSORT_UINT64_ARRAY , hits=4, 6.302 sec each
Time= 26.191 sec = 41.569%, QSORT_UINT64_38_ARRAY, hits=4, 6.548 sec each
从 64、48 和 32 位结果线性化到 38 位：
lsf 64 0.938 48 0.736 32 0.534 38 -> 0.6500
38 位 Radix_Sort 比 64 位排序快 35%，比 48 位排序快 17%。
即使是 40 位也会更快，5 个字节与每个 uint64 处理 6 个字节。
=========
目前最快，6 个 8 字节的 uint64[]，概括自： 32 MSbits used to sort uint64
// #############################################################################
// From: http://ideone.com/JHI0d9
// RadixSort--- for 48 MSB of uint64
typedef union
struct
uint32_t c6[256];
uint32_t c5[256];
uint32_t c4[256];
uint32_t c3[256];
uint32_t c2[256];
uint32_t c1[256];
;
uint32_t counts[256 * 6];
rscounts6_t;
// #############################################################################
// Patterned off of Radix_Sort_64 but looks only at the 48 MostSigBits
// 0XFFFF-FFFF-FFFF-0000 << Ignore the zeros, sort on 3 MostSigBytes
// Made for RGB48 stuffed into uint64 with 2 LeastSig bytes zero
// Get rid of the 7 and 8 level comps
uint64_t *radix_sort_48_msb(uint64_t *arrayA, uint32_t asize)
register uint64_t *array=arrayA; // Slam arg into Register!
register int ii; // Loop control
rscounts6_t counts;
memset(&counts, 0, 256 * 6 * sizeof(uint32_t));
uint64_t *cpy = (uint64_t *)malloc(asize * sizeof(uint64_t));
uint32_t o6=0, o5=0, o4=0, o3=0, o2=0, o1=0;
uint32_t t6, t5, t4, t3, t2, t1;
register uint32_t x;
// calculate counts
for(x = 0; x < asize; x++)
t6 = (array[x] >> 16) & 0xff;
t5 = (array[x] >> 24) & 0xff;
t4 = (array[x] >> 32) & 0xff;
t3 = (array[x] >> 40) & 0xff;
t2 = (array[x] >> 48) & 0xff;
t1 = (array[x] >> 56) & 0xff;
counts.c6[t6]++;
counts.c5[t5]++;
counts.c4[t4]++;
counts.c3[t3]++;
counts.c2[t2]++;
counts.c1[t1]++;
// convert counts to offsets
for(x = 0; x < 256; x++)
t6 = o6 + counts.c6[x];
t5 = o5 + counts.c5[x];
t4 = o4 + counts.c4[x];
t3 = o3 + counts.c3[x];
t2 = o2 + counts.c2[x];
t1 = o1 + counts.c1[x];
counts.c6[x] = o6;
counts.c5[x] = o5;
counts.c4[x] = o4;
counts.c3[x] = o3;
counts.c2[x] = o2;
counts.c1[x] = o1;
o6 = t6;
o5 = t5;
o4 = t4;
o3 = t3;
o2 = t2;
o1 = t1;
// radix
for(x = 0; x < asize; x++)
t6 = (array[x] >> 16) & 0xff;
cpy[counts.c6[t6]] = array[x];
counts.c6[t6]++;
for(x = 0; x < asize; x++)
t5 = (cpy[x] >> 24) & 0xff;
array[counts.c5[t5]] = cpy[x];
counts.c5[t5]++;
for(x = 0; x < asize; x++)
t4 = (array[x] >> 32) & 0xff;
cpy[counts.c4[t4]] = array[x];
counts.c4[t4]++;
for(x = 0; x < asize; x++)
t3 = (cpy[x] >> 40) & 0xff;
array[counts.c3[t3]] = cpy[x];
counts.c3[t3]++;
for(x = 0; x < asize; x++)
t2 = (array[x] >> 48) & 0xff;
cpy[counts.c2[t2]] = array[x];
counts.c2[t2]++;
for(x = 0; x < asize; x++)
t1 = (cpy[x] >> 56) & 0xff;
array[counts.c1[t1]] = cpy[x];
counts.c1[t1]++;
free(cpy);
return array;
// End radix_sort_48_msb().
===================================
再次感谢 Rcgldr 提供创新的编程建议！ 我没有使用 10、10、9、9，而是使用带有 [4][10]
的快速 32 位模式它可以工作，但它比 48 MSB 排序慢很多，40 MSB 为 737 毫秒，48 MSB 为 588 毫秒。 :(
也许我编码不好。
Time= 6.108 sec = 33.668%, QSORT_UINT64_ARRAY , hits=1
Time= 3.060 sec = 16.866%, RADIX_SORT_UINT64_REG, hits=4, 0.765 sec each
Time= 2.947 sec = 16.241%, RADIX_SORT_UINT64_40R, hits=4, 0.737 sec each < SLOW
Time= 2.354 sec = 12.973%, RADIX_SORT_UINT64_48R, hits=4, 0.588 sec each
Time= 1.542 sec = 8.498%, RADIX_SORT_UINT64_32R, hits=4, 0.385 sec each
Time= 0.769 sec = 4.236%, RADIX_SORT_64 , hits=1
测试：
创建一个随机的 uint64_t[36M] 主数组 使用 qsort 和已知良好的 Radix Sort 进行排序，radixSort() 以创建标准数组 比较 Qsort 和 radixSort() 结果的相同性。 使用 32、40、48 和 64 MSB 基数排序对主副本进行排序 在屏蔽掉忽略的 LSB 后，将每个测试排序与标准进行比较代码如下：
//=============================================================================
// From code submitted by rcgldr, Feb 8 2020
// Optimized to use Registers and to sort on 40 MSBs, ignoring 24 LSBs
void radix_sort_r64_40(uint64_t *pData, uint64_t *pTemp, size_t count,
EV_TIME_STR *tsa)
size_t mIndex[4][1024] = 0 ; /* index matrix */
size_t * pmIndex; /* ptr to row of matrix */
size_t i, j, m, n;
uint64_t u;
if(tsa) time_event(E_RADIX_SORT_UINT64_40R, tsa, E_TIME_EVENT, 1, 0);
for (i = 0; i < count; i++) /* generate histograms */
u = pData[i];
mIndex[3][(u >> 24) & 0x3ff]++;
mIndex[2][(u >> 34) & 0x3ff]++;
mIndex[1][(u >> 44) & 0x3ff]++;
mIndex[0][(u >> 54) & 0x3ff]++;
for (j = 0; j < 4; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 1024; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (i = 0; i < count; i++) /* radix sort */
u = pData[i];
pTemp[mIndex[3][(u >> 24) & 0x3ff]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[2][(u >> 34) & 0x3ff]++] = u;
for (i = 0; i < count; i++)
u = pData[i];
pTemp[mIndex[1][(u >> 44) & 0x3ff]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[0][(u >> 54) & 0x3ff]++] = u;
// End Radix_Sort_R64_40().
这是 FAST 32 MSB 版本和克隆的 40 MSB 慢速版本之间的唯一差异。
Unique lines from "~/tmp/radix.sort.32.c":
02) void radix_sort_r64_32(uint64_t *pData, uint64_t *pTemp, size_t count,
05) size_t mIndex[4][256] = 0 ; /* index matrix */
09) if(tsa) time_event(E_RADIX_SORT_UINT64_32R, tsa, E_TIME_EVENT, 1, 0);
13) mIndex[3][(u >> 32) & 0xff]++; // B4
14) mIndex[2][(u >> 40) & 0xff]++; // B5
15) mIndex[1][(u >> 48) & 0xff]++; // B6
16) mIndex[0][(u >> 56) & 0xff]++; // B7
22) for (i = 0; i < 256; i++)
31) pTemp[mIndex[3][(u >> 32) & 0xff]++] = u;
35) pData[mIndex[2][(u >> 40) & 0xff]++] = u;
39) pTemp[mIndex[1][(u >> 48) & 0xff]++] = u;
43) pData[mIndex[0][(u >> 56) & 0xff]++] = u;
Unique lines from "~/tmp/radix.sort.40.c":
01) void radix_sort_r64_40(uint64_t *pData, uint64_t *pTemp, size_t count,
04) size_t mIndex[4][1024] = 0 ; /* index matrix */
08) if(tsa) time_event(E_RADIX_SORT_UINT64_40R, tsa, E_TIME_EVENT, 1, 0);
12) mIndex[3][(u >> 24) & 0x3ff]++;
13) mIndex[2][(u >> 34) & 0x3ff]++;
14) mIndex[1][(u >> 44) & 0x3ff]++;
15) mIndex[0][(u >> 54) & 0x3ff]++;
21) for (i = 0; i < 1024; i++)
30) pTemp[mIndex[3][(u >> 24) & 0x3ff]++] = u;
34) pData[mIndex[2][(u >> 34) & 0x3ff]++] = u;
38) pTemp[mIndex[1][(u >> 44) & 0x3ff]++] = u;
42) pData[mIndex[0][(u >> 54) & 0x3ff]++] = u;
【问题讨论】：
这是个问题吗？ 一些不同（但不相关）的点：首先，在 C 中你不需要而且真的 shouldn't cast the result ofmalloc。同样要避免的是magic numbers（尤其是当sizeof 运算符可以得到“幻数”时，如sizeof counts）。而register 只是一个建议，如果编译器愿意，编译器根本不必费心（除了它必须遵循链接语义，比如没有获取地址等）。
@rcgldr - Skylake i7-6700K CPU @ 4.00GHz, cache_alignment : 64, AVX2, 64 GB DDR4 3.2GHz, CL16 (F4-3200C16-16GVK) 设置为 Interleaved 以实现更快的大块顺序。 4x32k inst/dat L1、4x256k L2、8M L3
@rcgldr - 您提到缓存一致性是 10 位 bin 与 8 位的性能下降的可能原因。写入是在总数组大小范围内从 A 点到 B 点完成的，所以它不依赖于数组大小并且独立于 bin 位吗？为了测试这一点，可以将阵列大小从 1M 更改为 32M，然后查看性能下降的地方。这可能会显示一个性能悬崖，指示缓存被淹没的位置。新的 AMD 3950x 具有 64 MB L3 和总共 8 MB L2！希望我有一个可以测试！
@BrianP007 - 缓存一致性 - 与内核或处理器之间的共享内存有关，这里不是问题。我指的是 L1 缓存的大小，即每个内核 32KB（6700K 和 3770K 的缓存大小相同）。 [8][256] 64 位索引的矩阵占用 16KB 空间，留出空间 L1 缓存用于数组读取。 [4][1024] 64 位索引的矩阵占用 32KB 空间，在初始计数过程中从 size_t 切换到 uint32_t 会有所帮助。在基数排序期间，一次只使用一行，所以不是问题。
【参考方案1】：
可以使用 10 bit, 10 bit, 9 bit, 9 bit 计数 == 偏移在 4 遍中完成 38 msb 的基数排序。更改您的代码以使用大小：c1[1024]、c2[1024]、c3[512]、c4[512] 和计数初始化以使用移位和掩码 (...>>54)&0x3ff ,(.. .>>44)&0x3ff, (...>>35)&0x1ff,(...>>26)&0x1ff，并对其余代码进行类似的更改。
尝试使用此代码进行比较。更改您的代码或将计时器内容添加到此代码中：
void RadixSort(uint64_t *a, uint64_t *b, size_t count)
uint32_t mIndex[4][1024] = 0 ; /* index matrix */
uint32_t * pmIndex; /* ptr to row of matrix */
uint32_t i, j, m, n;
uint64_t u;
for (i = 0; i < count; i++) /* generate histograms */
u = a[i];
mIndex[3][(u >> 26) & 0x1ff]++;
mIndex[2][(u >> 35) & 0x1ff]++;
mIndex[1][(u >> 44) & 0x3ff]++;
mIndex[0][(u >> 54) & 0x3ff]++;
for (j = 0; j < 2; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 1024; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (j = 2; j < 4; j++)
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 512; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
pmIndex = mIndex[3];
for (i = 0; i < count; i++) /* radix sort */
u = a[i];
b[pmIndex[(u >> 26) & 0x1ff]++] = u;
pmIndex = mIndex[2];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 35) & 0x1ff]++] = u;
pmIndex = mIndex[1];
for (i = 0; i < count; i++)
u = a[i];
b[pmIndex[(u >> 44) & 0x3ff]++] = u;
pmIndex = mIndex[0];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 54) & 0x3ff]++] = u;
更快的方法是先按最高有效位 10 位排序，创建 1024 个 bin，然后对 1024 个 bin 排序，10,9,9 位字段，最低有效位在前。这将加快排序速度，因为 1024 个 bin 中的每一个都适合缓存，从而减少了所有随机访问写入的开销。注意 - aIndex 的大小为 1025，用于跟踪最后一个 bin 的大小。
void RadixSort3(uint64_t *, uint64_t *, size_t);
/* split array into 1024 bins according to most significant 10 bits */
void RadixSort(uint64_t *a, uint64_t *b, size_t count)
uint32_t aIndex[1025] = 0; /* index array */
uint32_t i, m, n;
for(i = 0; i < count; i++) /* generate histogram */
aIndex[(a[i] >> 54)]++;
n = 0; /* convert to indices */
for (i = 0; i < 1025; i++)
m = aIndex[i];
aIndex[i] = n;
n += m;
for(i = 0; i < count; i++) /* sort by ms 10 bits */
b[aIndex[a[i]>>54]++] = a[i];
for(i = 1024; i; i--) /* restore aIndex */
aIndex[i] = aIndex[i-1];
aIndex[0] = 0;
for(i = 0; i < 1024; i++) /* radix sort the 1024 bins */
RadixSort3(&b[aIndex[i]], &a[aIndex[i]], aIndex[i+1]-aIndex[i]);
void RadixSort3(uint64_t *a, uint64_t *b, size_t count)
uint32_t mIndex[3][1024] = 0 ; /* index matrix */
uint32_t * pmIndex; /* ptr to row of matrix */
uint32_t i, j, m, n;
uint64_t u;
for (i = 0; i < count; i++) /* generate histograms */
u = a[i];
mIndex[2][(u >> 26) & 0x1ff]++;
mIndex[1][(u >> 35) & 0x1ff]++;
mIndex[0][(u >> 44) & 0x3ff]++;
for (j = 0; j < 1; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 1024; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (j = 1; j < 3; j++)
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 512; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
pmIndex = mIndex[2];
for (i = 0; i < count; i++) /* radix sort */
u = a[i];
b[pmIndex[(u >> 26) & 0x1ff]++] = u;
pmIndex = mIndex[1];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 35) & 0x1ff]++] = u;
pmIndex = mIndex[0];
for (i = 0; i < count; i++)
u = a[i];
b[pmIndex[(u >> 44) & 0x3ff]++] = u;
【讨论】：
【参考方案2】：40 MSB 排序的性能很差，比 48 MSB 版本慢很多。
所以我尝试了一个 [6][64]，36 MSB 的版本，它以相对较快的 48 MSB 为蓝本。
我知道我想要 38 位，而不是 36。优化的方法是如何将 XY 位置和该点的属性打包到 uint64_t 中。使用 8k X 和 Y 坐标，XY 占用 13+13 位，数据不超过 64-26=38 位。
当前数据最大为 34 或 35 位，因此 36 应该可以工作。
表演如下：
Time= 6.104 sec = 30.673%, QSORT_UINT64_ARRAY , hits=1
Time= 3.117 sec = 15.663%, RADIX_SORT_UINT64_REG, hits=4, 0.779 sec each
Time= 2.931 sec = 14.731%, RADIX_SORT_UINT64_40R, hits=4, 0.733 sec each
Time= 2.269 sec = 11.401%, RADIX_SORT_UINT64_48R, hits=4, 0.567 sec each
Time= 1.663 sec = 8.359%, RADIX_SORT_UINT64_36R, hits=4, 0.416 sec each < FAST
Time= 1.516 sec = 7.620%, RADIX_SORT_UINT64_32R, hits=4, 0.379 sec each
Time= 0.734 sec = 3.689%, RADIX_SORT_64 , hits=1
它比 48 位代码快 27%。
而且，如果 36 位太紧，它看起来应该能够扩展到 [6][7] 以对 42 个 MSB 进行排序！
这是完整的代码：
void radix_sort_r64_36(uint64_t *pData, uint64_t *pTemp, size_t count,
EV_TIME_STR *tsa)
size_t mIndex[6][64] = 0 ; /* index matrix */
size_t * pmIndex; /* ptr to row of matrix */
size_t i, j, m, n;
uint64_t u;
if(tsa) time_event(E_RADIX_SORT_UINT64_36R, tsa, E_TIME_EVENT, 1, 0);
// 64 -- 56 48 40 32 24 16 -- 8 bits each
// 64 -- 58 52 46 40 34 28 -- 6 bits each
for (i = 0; i < count; i++) /* generate histograms */
u = pData[i]; // Igonores Nibbles 0, 1 & 2
mIndex[5][(u >> 28) & 0x3F]++; // N2
mIndex[4][(u >> 34) & 0x3F]++; // N3
mIndex[3][(u >> 40) & 0x3F]++; // N4
mIndex[2][(u >> 46) & 0x3F]++; // N5
mIndex[1][(u >> 52) & 0x3F]++; // N6
mIndex[0][(u >> 58) & 0x3F]++; // N7
for (j = 0; j < 6; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 64; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (i = 0; i < count; i++) /* radix sort */
u = pData[i];
pTemp[mIndex[5][(u >> 28) & 0x3F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[4][(u >> 34) & 0x3F]++] = u;
for (i = 0; i < count; i++)
u = pData[i];
pTemp[mIndex[3][(u >> 40) & 0x3F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[2][(u >> 46) & 0x3F]++] = u;
for (i = 0; i < count; i++)
u = pData[i];
pTemp[mIndex[1][(u >> 52) & 0x3F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[0][(u >> 58) & 0x3F]++] = u;
// End Radix_Sort_R64_36().
具有 48 MSB 功能的唯一差异：
Unique lines from "/home/brianp/tmp/radix.sort.36.c":
01) void radix_sort_r64_36(uint64_t *pData, uint64_t *pTemp, size_t count,
04) size_t mIndex[6][64] = 0 ; /* index matrix */
08) if(tsa) time_event(E_RADIX_SORT_UINT64_36R, tsa, E_TIME_EVENT, 1, 0);
11) mIndex[5][(u >> 28) & 0x3F]++; // N2
12) mIndex[4][(u >> 34) & 0x3F]++; // N3
13) mIndex[3][(u >> 40) & 0x3F]++; // N4
14) mIndex[2][(u >> 46) & 0x3F]++; // N5
15) mIndex[1][(u >> 52) & 0x3F]++; // N6
16) mIndex[0][(u >> 58) & 0x3F]++; // N7
22) for (i = 0; i < 64; i++)
31) pTemp[mIndex[5][(u >> 28) & 0x3F]++] = u;
35) pData[mIndex[4][(u >> 34) & 0x3F]++] = u;
39) pTemp[mIndex[3][(u >> 40) & 0x3F]++] = u;
43) pData[mIndex[2][(u >> 46) & 0x3F]++] = u;
47) pTemp[mIndex[1][(u >> 52) & 0x3F]++] = u;
51) pData[mIndex[0][(u >> 58) & 0x3F]++] = u;
Unique lines from "/home/brianp/tmp/radix.sort.48.c":
01) void radix_sort_r64_48(uint64_t *pData, uint64_t *pTemp, size_t count,
04) size_t mIndex[6][256] = 0 ; /* index matrix */
08) if(tsa) time_event(E_RADIX_SORT_UINT64_48R, tsa, E_TIME_EVENT, 1, 0);
14) mIndex[5][(u >> 16) & 0xff]++; // B2
15) mIndex[4][(u >> 24) & 0xff]++; // B3
16) mIndex[3][(u >> 32) & 0xff]++; // B4
17) mIndex[2][(u >> 40) & 0xff]++; // B5
18) mIndex[1][(u >> 48) & 0xff]++; // B6
19) mIndex[0][(u >> 56) & 0xff]++; // B7
25) for (i = 0; i < 256; i++)
34) pTemp[mIndex[5][(u >> 16) & 0xff]++] = u;
38) pData[mIndex[4][(u >> 24) & 0xff]++] = u;
42) pTemp[mIndex[3][(u >> 32) & 0xff]++] = u;
46) pData[mIndex[2][(u >> 40) & 0xff]++] = u;
50) pTemp[mIndex[1][(u >> 48) & 0xff]++] = u;
54) pData[mIndex[0][(u >> 56) & 0xff]++] = u;
【讨论】：
【参考方案3】：为了完整性，6 bin，7 bit/bin 42 MSB Radix Sort 可以正常工作，性能与 48 MSB 和 36 MSB 版本完全一致。
Time= 6.334 sec = 25.435%, QSORT_UINT64_ARRAY , hits=1
Time= 3.519 sec = 14.131%, RADIX_SORT_UINT64_REG, hits=4, 0.880 sec each
Time= 3.273 sec = 13.145%, RADIX_SORT_UINT64_40R, hits=4, 0.818 sec each < anomaly
Time= 2.680 sec = 10.764%, RADIX_SORT_UINT64_48R, hits=4, 0.670 sec each
Time= 2.302 sec = 9.246%, RADIX_SORT_UINT64_42R, hits=4, 0.576 sec each < NEW
Time= 2.025 sec = 8.132%, RADIX_SORT_UINT64_36R, hits=4, 0.506 sec each
Time= 1.767 sec = 7.094%, RADIX_SORT_UINT64_32R, hits=4, 0.442 sec each
Time= 0.955 sec = 3.835%, RADIX_SORT_64 , hits=1
谁能解释为什么 40 MSB 排序比 48 MSB 排序慢得多？
完整代码：
void radix_sort_r64_42(uint64_t *pData, uint64_t *pTemp, size_t count,
EV_TIME_STR *tsa)
size_t mIndex[6][128] = 0 ; /* index matrix */
size_t * pmIndex; /* ptr to row of matrix */
size_t i, j, m, n;
uint64_t u;
if(tsa) time_event(E_RADIX_SORT_UINT64_42R, tsa, E_TIME_EVENT, 1, 0);
// 64 -- 56 48 40 32 24 16 -- 8 bits each
// 64 -- 57 50 43 36 29 22 -- 7 bits each
// 64 -- 58 52 46 40 34 28 -- 6 bits each
for (i = 0; i < count; i++) /* generate histograms */
u = pData[i]; // Igonores Nibbles 0, 1 & 2
mIndex[5][(u >> 22) & 0x7F]++; // N2
mIndex[4][(u >> 29) & 0x7F]++; // N3
mIndex[3][(u >> 36) & 0x7F]++; // N4
mIndex[2][(u >> 43) & 0x7F]++; // N5
mIndex[1][(u >> 50) & 0x7F]++; // N6
mIndex[0][(u >> 57) & 0x7F]++; // N7
for (j = 0; j < 6; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 128; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (i = 0; i < count; i++) /* radix sort */
u = pData[i];
pTemp[mIndex[5][(u >> 22) & 0x7F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[4][(u >> 29) & 0x7F]++] = u;
for (i = 0; i < count; i++)
u = pData[i];
pTemp[mIndex[3][(u >> 36) & 0x7F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[2][(u >> 43) & 0x7F]++] = u;
for (i = 0; i < count; i++)
u = pData[i];
pTemp[mIndex[1][(u >> 50) & 0x7F]++] = u;
for (i = 0; i < count; i++)
u = pTemp[i];
pData[mIndex[0][(u >> 57) & 0x7F]++] = u;
// End Radix_Sort_R64_42().
添加 36 MSB 与 42 MSB 的差异版本
Unique lines from "~/tmp/radix.sort.36.c":
01) void radix_sort_r64_36(uint64_t *pData, uint64_t *pTemp, size_t count,
04) size_t mIndex[6][64] = 0 ; /* index matrix */
11) mIndex[5][(u >> 28) & 0x3F]++; // N2
12) mIndex[4][(u >> 34) & 0x3F]++; // N3
13) mIndex[3][(u >> 40) & 0x3F]++; // N4
14) mIndex[2][(u >> 46) & 0x3F]++; // N5
15) mIndex[1][(u >> 52) & 0x3F]++; // N6
16) mIndex[0][(u >> 58) & 0x3F]++; // N7
22) for (i = 0; i < 64; i++)
31) pTemp[mIndex[5][(u >> 28) & 0x3F]++] = u;
35) pData[mIndex[4][(u >> 34) & 0x3F]++] = u;
39) pTemp[mIndex[3][(u >> 40) & 0x3F]++] = u;
43) pData[mIndex[2][(u >> 46) & 0x3F]++] = u;
47) pTemp[mIndex[1][(u >> 52) & 0x3F]++] = u;
51) pData[mIndex[0][(u >> 58) & 0x3F]++] = u;
19 Unique lines from "~/tmp/radix.sort.42.c":
01) void radix_sort_r64_42(uint64_t *pData, uint64_t *pTemp, size_t count,
04) size_t mIndex[6][128] = 0 ; /* index matrix */
10) // 64 -- 56 48 40 32 24 16 -- 8 bits each
11) // 64 -- 57 50 43 36 29 22 -- 7 bits each
12) // 64 -- 58 52 46 40 34 28 -- 6 bits each
15) mIndex[5][(u >> 22) & 0x7F]++; // N2
16) mIndex[4][(u >> 29) & 0x7F]++; // N3
17) mIndex[3][(u >> 36) & 0x7F]++; // N4
18) mIndex[2][(u >> 43) & 0x7F]++; // N5
19) mIndex[1][(u >> 50) & 0x7F]++; // N6
20) mIndex[0][(u >> 57) & 0x7F]++; // N7
26) for (i = 0; i < 128; i++)
35) pTemp[mIndex[5][(u >> 22) & 0x7F]++] = u;
39) pData[mIndex[4][(u >> 29) & 0x7F]++] = u;
43) pTemp[mIndex[3][(u >> 36) & 0x7F]++] = u;
47) pData[mIndex[2][(u >> 43) & 0x7F]++] = u;
51) pTemp[mIndex[1][(u >> 50) & 0x7F]++] = u;
55) pData[mIndex[0][(u >> 57) & 0x3F]++] = u;
【讨论】：
【参考方案4】：另外2个排序算法测试结果：
Time= 6.208 sec = 21.838%, QSORT_UINT64_ARRAY , hits=1
Time= 3.358 sec = 11.813%, RADIX_SORT_UINT64_REG, hits=4, 0.840 sec each
Time= 2.525 sec = 8.884%, RADIX_SORT_UI64_AA99 , hits=4, 0.631 sec each <NEW
Time= 2.509 sec = 8.825%, RADIX_SORT_UINT64_48R, hits=4, 0.627 sec each
Time= 2.461 sec = 8.658%, RADIX_SORT_UI64_1024 , hits=4, 0.615 sec each <NEW
Time= 2.223 sec = 7.822%, RADIX_SORT_UINT64_42R, hits=4, 0.556 sec each
Time= 2.215 sec = 7.791%, RADIX_SORT_UI64_40_85, hits=4, 0.554 sec each
Time= 1.930 sec = 6.788%, RADIX_SORT_UINT64_36R, hits=4, 0.482 sec each
Time= 1.710 sec = 6.014%, RADIX_SORT_UINT64_32R, hits=4, 0.427 sec each
Time= 0.915 sec = 3.220%, COMP_UINT64_ARRAYS , hits=32, 0.029 sec each
另一个没有 Firefox 占用系统的运行：
Time= 6.156 sec = 23.199%, QSORT_UINT64_ARRAY , hits=1
Time= 2.993 sec = 11.277%, RADIX_SORT_UINT64_REG, hits=4, 0.748 sec each
Time= 2.409 sec = 9.077%, RADIX_SORT_UI64_AA99 , hits=4, 0.602 sec each < NEW
Time= 2.330 sec = 8.778%, RADIX_SORT_UI64_1024 , hits=4, 0.582 sec each < NEW
Time= 2.241 sec = 8.443%, RADIX_SORT_UINT64_48R, hits=4, 0.560 sec each
Time= 2.124 sec = 8.002%, RADIX_SORT_UI64_40_85, hits=4, 0.531 sec each
Time= 1.982 sec = 7.468%, RADIX_SORT_UINT64_42R, hits=4, 0.495 sec each
Time= 1.725 sec = 6.499%, RADIX_SORT_UINT64_36R, hits=4, 0.431 sec each
Time= 1.507 sec = 5.677%, RADIX_SORT_UINT64_32R, hits=4, 0.377 sec each
Time= 0.889 sec = 3.348%, COMP_UINT64_ARRAYS , hits=32, 0.028 sec each
RADIX_SORT_UI64_AA99 在 4 个通道中使用 bin 位 [10, 10, 9, 9]。
RADIX_SORT_UI64_1024 首先按最高有效 10 位排序。
RADIX_SORT_UINT64_42R 使用 6 个 bin，每个 bin 有 7 位
注意决斗 40 MSB、8 bin 5 bit RADIX_SORT_UI64_40_85 和 42 MSB、7 bin 6 bit RADIX_SORT_UINT64_42R 的相对性能。 42 MSB 排序通常要快得多，尽管 40 位有时会将其排除在外。
编译器是 GCC，针对速度进行了优化：
gcc -Ofast -ffast-math -m64 -march=native -funroll-loops -fopenmp -flto -finline-functions -Wuninitialized ~/bin/pb.c -lm -o ~/bin/pb_a
有更好的编译器选项吗？
gcc 版本 7.4.1 20190905 [gcc-7-branch 修订版 275407] (SUSE Linux)
===========================
完整代码：
// =============================================================================
// Sort with bin bits 10, 10, 9,9
// From code by rcgldr *** Feb 8, 2020
void RadixSort_aa99(uint64_t *a, uint64_t *b, size_t count, EV_TIME_STR *tsa)
uint32_t mIndex[4][1024] = 0 ; /* index matrix */
uint32_t * pmIndex; /* ptr to row of matrix */
uint32_t i, j, m, n;
uint64_t u;
if(tsa) time_event(E_RADIX_SORT_UI64_AA99, tsa, E_TIME_EVENT, 1, 0);
for (i = 0; i < count; i++) /* generate histograms */
u = a[i];
mIndex[3][(u >> 26) & 0x1ff]++;
mIndex[2][(u >> 35) & 0x1ff]++;
mIndex[1][(u >> 44) & 0x3ff]++;
mIndex[0][(u >> 54) & 0x3ff]++;
for (j = 0; j < 2; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 1024; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (j = 2; j < 4; j++)
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 512; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
pmIndex = mIndex[3];
for (i = 0; i < count; i++) /* radix sort */
u = a[i];
b[pmIndex[(u >> 26) & 0x1ff]++] = u;
pmIndex = mIndex[2];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 35) & 0x1ff]++] = u;
pmIndex = mIndex[1];
for (i = 0; i < count; i++)
u = a[i];
b[pmIndex[(u >> 44) & 0x3ff]++] = u;
pmIndex = mIndex[0];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 54) & 0x3ff]++] = u;
// =============================================================================
/* split array into 1024 bins according to most significant 10 bits */
void RadixSort_1024(uint64_t *a, uint64_t *b, size_t count, EV_TIME_STR *tsa)
uint32_t aIndex[1025] = 0; /* index array */
uint32_t i, m, n;
if(tsa) time_event(E_RADIX_SORT_UI64_1024, tsa, E_TIME_EVENT, 1, 0);
for(i = 0; i < count; i++) /* generate histogram */
aIndex[(a[i] >> 54)]++;
n = 0; /* convert to indices */
for (i = 0; i < 1025; i++)
m = aIndex[i];
aIndex[i] = n;
n += m;
for(i = 0; i < count; i++) /* sort by ms 10 bits */
b[aIndex[a[i]>>54]++] = a[i];
for(i = 1024; i; i--) /* restore aIndex */
aIndex[i] = aIndex[i-1];
aIndex[0] = 0;
for(i = 0; i < 1024; i++) /* radix sort the 1024 bins */
RadixSort3(&b[aIndex[i]], &a[aIndex[i]], aIndex[i+1]-aIndex[i]);
void RadixSort3(uint64_t *a, uint64_t *b, size_t count)
uint32_t mIndex[3][1024] = 0 ; /* index matrix */
uint32_t * pmIndex; /* ptr to row of matrix */
uint32_t i, j, m, n;
uint64_t u;
for (i = 0; i < count; i++) /* generate histograms */
u = a[i];
mIndex[2][(u >> 26) & 0x1ff]++;
mIndex[1][(u >> 35) & 0x1ff]++;
mIndex[0][(u >> 44) & 0x3ff]++;
for (j = 0; j < 1; j++) /* convert to indices */
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 1024; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
for (j = 1; j < 3; j++)
pmIndex = mIndex[j];
n = 0;
for (i = 0; i < 512; i++)
m = pmIndex[i];
pmIndex[i] = n;
n += m;
pmIndex = mIndex[2];
for (i = 0; i < count; i++) /* radix sort */
u = a[i];
b[pmIndex[(u >> 26) & 0x1ff]++] = u;
pmIndex = mIndex[1];
for (i = 0; i < count; i++)
u = b[i];
a[pmIndex[(u >> 35) & 0x1ff]++] = u;
pmIndex = mIndex[0];
for (i = 0; i < count; i++)
u = a[i];
b[pmIndex[(u >> 44) & 0x3ff]++] = u;
【讨论】：
设置大小为 36.2M uint64_t 和 100% 不同，使用 mrand48() 生成两次，以上是关于在 uint64 中对 38 位进行基数排序？的主要内容，如果未能解决你的问题，请参考以下文章