数百万个UINT64 RGBZ图形像素的最快排序算法

Question

我正在对10多百万个带有RGB数据的uint64_t从.RAW文件进行排序，而我的C程序中79％的时间都花费在qsort上。我正在寻找针对这种特定数据类型的更快速的排序方法。

由于是原始图形数据，数字非常随机，并且大约80％是唯一的。不能期望部分排序或排好序的数据运行。 uint64_t内的4个uint16_t分别为R、G、B和零（可能是一个小计数<=〜20）。

我使用unsigned long long编写了我能想到的最简单的比较函数（您不能只是将它们相减）：

qsort(hpidx, num_pix, sizeof(uint64_t), comp_uint64); 
...
int comp_uint64(const void *a, const void *b)  {
    if(*((uint64_t *)a) > *((uint64_t *)b))  return(+1);
    if(*((uint64_t *)a) < *((uint64_t *)b))  return(-1);
    return(0);
}  // End Comp_uint64().

在StackExchange上有一个非常有趣的“编程谜题和代码高尔夫”，但是他们使用了float。然后有QSort、RecQuick、heap、stooge、tree、radix等。swenson/sort看起来很有趣，但没有（显而易见的）支持我的数据类型uint64_t。"快速排序"的时间最好。一些来源说系统qsort可以是任何东西，不一定是"快速排序"。C++ sort绕过了void指针的通用转换，在性能上比C实现了巨大的改进。必须有一种优化的方法将U8快速传递到64位处理器中。

---

系统/编译器信息：

我目前正在使用GCC和草莓珍珠岩。

gcc version 4.9.2 (x86_64-posix-sjlj, built by strawberryperl.com
Intel 2700K Sandy Bridge CPU, 32GB DDR3
windows 7/64 pro

gcc -D__USE_MINGW_ANSI_STDIO -O4 -ffast-math -m64 -Ofast -march=corei7-avx -mtune=corei7 -Ic:/bin/xxHash-master -Lc:/bin/xxHash-master c:/bin/stddev.c -o c:/bin/stddev.g6.exe 

尝试改进的qsort算法，QSORT()!

尝试使用Michael Tokarev的内联qsort算法。

“READY-TO-USE”？来自qsort.h文档。

-----------------------------
* Several ready-to-use examples:
 *
 * Sorting array of integers:
 * void int_qsort(int *arr, unsigned n) {
 * #define int_lt(a,b) ((*a)<(*b))
 *   QSORT(int, arr, n, int_lt);
--------------------------------

Change from type "int" to "uint64_t"
compile error on TYPE???
    
    c:/bin/bpbfct.c:586:8: error: expected expression before 'uint64_t'
      QSORT(uint64_t, hpidx, num_pix, islt);

我找不到一个真正的、编译的、可工作的示例程序，只有关于“一般概念”的注释。

#define QSORT_TYPE uint64_t 
#define islt(a,b) ((*a)<(*b))

uint64_t *QSORT_BASE; 
int QSORT_NELT;

hpidx=(uint64_t *) calloc(num_pix+2, sizeof(uint64_t));  // Hash . PIDX
QSORT_BASE = hpidx;
QSORT_NELT = num_pix;  // QSORT_LT is function QSORT_LT()
QSORT(uint64_t, hpidx, num_pix, islt);  
//QSORT(uint64_t *, hpidx, num_pix, QSORT_LT);  // QSORT_LT mal-defined?
//qsort(hpidx, num_pix, sizeof(uint64_t), comp_uint64); // << WORKS

这些“即用型”示例使用了int、char *和struct elt类型。但是uint64_t不是一种类型吗？尝试使用long long。

QSORT(long long, hpidx, num_pix, islt); 
c:/bin/bpbfct.c:586:8: error: expected expression before 'long'
 QSORT(long long, hpidx, num_pix, islt);

---

下一步尝试: RADIX_SORT:

结果: RADIX_SORT十分彻底！

  I:\br3\pf.249465>grep "Event" bb12.log | grep -i Sort       
 << 1.40 sec average
4) Time=1.411 sec    = 49.61%, Event RADIX_SORT        , hits=1
4) Time=1.396 sec    = 49.13%, Event RADIX_SORT        , hits=1
4) Time=1.392 sec    = 49.15%, Event RADIX_SORT        , hits=1
16) Time=1.414 sec    = 49.12%, Event RADIX_SORT        , hits=1

I:\br3\pf.249465>grep "Event" bb11.log | grep -i Sort 
 << 5.525 sec average  = 3.95 time slower
4) Time=5.538 sec    = 86.34%, Event QSort             , hits=1
4) Time=5.519 sec    = 79.41%, Event QSort             , hits=1
4) Time=5.519 sec    = 79.02%, Event QSort             , hits=1
4) Time=5.563 sec    = 79.49%, Event QSort             , hits=1
4) Time=5.684 sec    = 79.83%, Event QSort             , hits=1
4) Time=5.509 sec    = 79.30%, Event QSort             , hits=1

比起开箱即用的qsort排序算法，这个算法快了3.94倍！

更重要的是，这里有实际可行的代码，而不仅仅是某位大师给你80%的内容，假设你知道他们所知道的一切，并且可以填补另外20%的空缺。

绝妙的解决方案！感谢Louis Ricci！

- BrianP007

---

如果你所描述的数据是随机的，那么我会说qsort已经是其中性能最稳定的实现之一了。 - Jason Hu

你可以简单地使用C++的sort函数吗？你可以将它放入一个单独的.cpp文件中，并使用extern "C"来使得你的其余代码仍然可以保持在C语言中。 - Adam

@user3386109，您想到的是RGBA。我认为Z指的是深度信息。无论如何，这与问题无关。 - Adam

计数排序？也许可以。 - Nikos M.

RGBZ是为了RGB_ZERO_而设计的。我将3个UINT16压缩成一个UINT64，因为UINT32太小了，K和R都没有带有16位量子的数码相机，也没有UINT48。这是多么严重的疏忽啊！在加载RGB之前，我清除了所有8个字节。我还在最后2个字节中添加了计数和其他内容。 - BrianP007

Answer 1

我建议使用基数排序，基数为8位。对于64位的值，一个经过优化的基数排序将会遍历列表9次（1次用于预计算计数和偏移量，8次用于64位/8位）。时间复杂度为9*N，空间复杂度为2*N（使用一个阴影数组）。

下面是一个经过优化的基数排序的示例：

typedef union {
    struct {
        uint32_t c8[256];
        uint32_t c7[256];
        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 * 8];
} rscounts_t;

uint64_t * radixSort(uint64_t * array, uint32_t size) {
    rscounts_t counts;
    memset(&counts, 0, 256 * 8 * sizeof(uint32_t));
    uint64_t * cpy = (uint64_t *)malloc(size * sizeof(uint64_t));
    uint32_t o8=0, o7=0, o6=0, o5=0, o4=0, o3=0, o2=0, o1=0;
    uint32_t t8, t7, t6, t5, t4, t3, t2, t1;
    uint32_t x;
    // calculate counts
    for(x = 0; x < size; x++) {
        t8 = array[x] & 0xff;
        t7 = (array[x] >> 8) & 0xff;
        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.c8[t8]++;
        counts.c7[t7]++;
        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++) {
        t8 = o8 + counts.c8[x];
        t7 = o7 + counts.c7[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.c8[x] = o8;
        counts.c7[x] = o7;
        counts.c6[x] = o6;
        counts.c5[x] = o5;
        counts.c4[x] = o4;
        counts.c3[x] = o3;
        counts.c2[x] = o2;
        counts.c1[x] = o1;
        o8 = t8; 
        o7 = t7; 
        o6 = t6; 
        o5 = t5; 
        o4 = t4; 
        o3 = t3; 
        o2 = t2; 
        o1 = t1;
    }
    // radix
    for(x = 0; x < size; x++) {
        t8 = array[x] & 0xff;
        cpy[counts.c8[t8]] = array[x];
        counts.c8[t8]++;
    }
    for(x = 0; x < size; x++) {
        t7 = (cpy[x] >> 8) & 0xff;
        array[counts.c7[t7]] = cpy[x];
        counts.c7[t7]++;
    }
    for(x = 0; x < size; x++) {
        t6 = (array[x] >> 16) & 0xff;
        cpy[counts.c6[t6]] = array[x];
        counts.c6[t6]++;
    }
    for(x = 0; x < size; x++) {
        t5 = (cpy[x] >> 24) & 0xff;
        array[counts.c5[t5]] = cpy[x];
        counts.c5[t5]++;
    }
    for(x = 0; x < size; x++) {
        t4 = (array[x] >> 32) & 0xff;
        cpy[counts.c4[t4]] = array[x];
        counts.c4[t4]++;
    }
    for(x = 0; x < size; x++) {
        t3 = (cpy[x] >> 40) & 0xff;
        array[counts.c3[t3]] = cpy[x];
        counts.c3[t3]++;
    }
    for(x = 0; x < size; x++) {
        t2 = (array[x] >> 48) & 0xff;
        cpy[counts.c2[t2]] = array[x];
        counts.c2[t2]++;
    }
    for(x = 0; x < size; x++) {
        t1 = (cpy[x] >> 56) & 0xff;
        array[counts.c1[t1]] = cpy[x];
        counts.c1[t1]++;
    }
    free(cpy);
    return array;
}

编辑：此实现基于 JavaScript 版本，JavaScript 中排序 32 位有符号整数数组的最快方法？

以下是使用 C 实现基数排序的 IDEONE 链接：http://ideone.com/JHI0d9

- Louis Ricci

---

由于“size”在百万级别，因此“malloc”将非常巨大。 - user3629249

也许一个更简单的计数排序（类似于O(n)复杂度）也可以。 - Nikos M.

@user3629249 - 1000万个uint64 * 8字节 = 8000万字节 ~= 80MB，我认为现代计算机可以处理80MB的分配。有原地版本的基数排序可避免阴影数组分配，但这个版本似乎足够快。 - Louis Ricci

@NikosM. - 我认为你会发现计数排序会变成基数排序。被排序的元素是uint64值，因此对于计数排序，您需要一个大小为2^64的计数数组（对于这个特定问题，其中64位中的16位为零，因此您只需要2^48，但仍然非常巨大）。 - Louis Ricci

使用16位基数排序，尽管使用的步骤更少，但速度会变慢吗？ - Simd

Answer 2

我看到了几个选项，按照难易程度排序：

• 使用-flto开关启用链接时优化。这可能会让编译器内联您的比较函数。这么简单不尝试太可惜。
• 如果LTO没有效果，您可以使用内联qsort实现，例如Michael Tokarev's inline qsort。 此页面建议使用此方法可以提高2倍性能，仅因为编译器可以内联比较函数。
• 使用C++的std::sort。我知道你的代码是用C写的，但你可以创建一个只排序并提供C接口的小模块。您已经在使用具有出色C++支持的工具链。
• 尝试swenson/sort库。它实现了许多算法，因此您可以使用最适合您数据的算法。它似乎是可内联的，并且他们声称比qsort更快。
• 查找另一个排序库。能够执行Louis' Radix Sort的库是一个很好的建议。

请注意，您也可以使用单个分支进行比较，而不是两个。只需找出哪一个更大，然后减去即可。

- Adam

---

两个 uint64_t 的相减可能会生成一个不适合于 int 的值。如果没有这个问题，你可以消除所有的分支语句，直接返回相减的结果。 - Mark Ransom

减去两个uint64_t可能会生成一个不适合int的值。我发现一半的减法生成了负结果，这破坏了整个排序。双步骤过程旨在规避此问题。无论如何，我已经切换到RADIX_SORT;速度提高了4倍！ - BrianP007

@BrianP007 另一个无符号整数相减的问题是它们不能为负数。 - Adam

@BrianP007 如果你决定采用Louis的答案，你应该点击他的答案旁边的复选标记将其标记为“已接受”。 - Adam

Answer 3

对于一些编译器/平台，以下代码是无分支且更快的，尽管与 OP 原始代码没有太大区别。

int comp_uint64_b(const void *a, const void *b)  {
    return 
     (*((uint64_t *)a) > *((uint64_t *)b)) - 
     (*((uint64_t *)a) < *((uint64_t *)b));
}

- chux - Reinstate Monica

Answer 4

也许使用三目运算符而不是if语句会使事情更快一些。

- Lambder

---

微观优化通常是不好的 - 让编译器自己处理。 - tonysdg

这非常不可能。任何编写得当的编译器都将为传统的if ... else和相应的三元（即?:）表达式生成相同的代码。 - Jim Mischel