【发布时间】：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 of
malloc。同样要避免的是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 会有所帮助。在基数排序期间，一次只使用一行，所以不是问题。
标签： c sorting radix-sort