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

Posted

技术标签:

【中文标题】数百万 UINT64 RGBZ 图形像素的最快排序算法【英文标题】:Fastest sort algorithm for millions of UINT64 RGBZ graphics pixels 【发布时间】:2015-11-18 05:38:15 【问题描述】:

我正在使用 .RAW 文件中的 RGB 数据对 10+ 百万个 uint64_ts 进行排序,并且我 79% 的 C 程序时间都花在了 qsort 上。我正在为这种特定的数据类型寻找一种更快的排序方式。

作为 RAW 图形数据,这些数字非常随机,大约 80% 是唯一的。预计不会对已排序的数据进行部分排序或运行。 uint64_t 中的 4 个 uint16_ts 是 R、G、B 和零(可能是一个小计数

我有我能想到的最简单的比较函数,使用unsigned long longs(你不能只减去它们):

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 上有一个非常有趣的“Programming Puzzles & Code Golf”,但他们使用了floats。然后是QSort、RecQuick、heap、stooge、tree、radix...

swenson/sort 看起来很有趣,但对我的数据类型uint64_t 没有(明显的)支持。 “快速排序”时间是最好的。一些消息来源说系统qsort 可以是任何东西,不一定是“快速排序”。

C++ 排序绕过了 void 指针的通用转换,并实现了相对于 C 的性能的巨大改进。必须有一种优化的方法来通过 64 位处理器以扭曲速度猛击 U8。


系统/编译器信息:

我目前正在将 GCC 与 Strawberry Perl 一起使用

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 

第一次尝试更好的qsortQSORT()

尝试使用 Michael Tokarev 的内联 qsort

“准备使用”?来自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);

我找不到一个真正的、编译的、工作的示例程序,只是带有“一般概念”的 cmets

#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

“即用型”示例使用 intchar *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);

下一次尝试:RADIXSORT:

结果:RADIX_SORT 是 RADICAL!

  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%。

绝妙的解决方案!谢谢路易斯·里奇!

【问题讨论】:

如果您的数据像您描述的那样随机,那么我会说 qsort 将是性能最稳定的实现之一。 你能简单地使用 C++ 排序吗?您可以将其放入带有 extern "C" 的单独 .cpp 文件中,以便您的其余代码可以保留在 C 中。 @user3386109 您正在考虑 RGBA。我认为 Z 指的是深度信息。无论哪种方式,它都与问题无关。 counting sort?也许 RGBZ 用于 RGB_ZERO_ 我将 3 个 UINT16 粉碎成 UINT64,因为 UINT32 太小,而且 K 和 R 都没有 16 位量子的数码相机,没有 UINT48。多么严重的疏忽!在加载 RGB 之前,我清除了所有 8 个字节。我还在最后 2 个字节中添加了计数和其他内容。 【参考方案1】:

我会使用带有 8 位基数的基数排序。对于 64 位值,一个经过良好优化的基数排序将必须遍历列表 9 次(一次用于预先计算计数和偏移量,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 版本 Fastest way to sort 32bit signed integer arrays in JavaScript?

这是用于 C 基数排序的 IDEONE http://ideone.com/JHI0d9

【讨论】:

基数排序没有太多用处,但对于像素排序来说几乎是完美的。不幸的是,有一次我必须对像素进行排序,直到我实现了其他东西时才意识到这一点。 “大小”以百万计,“malloc”将是巨大的。 也许更简单的计数排序(类似O(n) 复杂度)也可以 @user3629249 - 1000 万个 uint64 * 8 字节 = 8000 万字节 ~= 80MB,我认为现代计算机可以处理 80MB 的分配。有基数排序的就地版本来避免影子数组分配,但这个版本似乎足够快。 @NikosM。 - 我想你会发现计数排序会变成基数排序。被排序的元素是 uint64 值,因此对于计数排序,您需要一个大小为 2^64 的计数数组(对于这个特定问题,这 64 位中的 16 位为零,因此您只需要 2^48,这仍然很大)。 【参考方案2】:

我看到了一些选项,大致按从易到难的顺序排列。

使用-flto 开关启用链接时间优化。这可能让编译器内联您的比较函数。太容易不尝试了。 如果 LTO 没有效果,您可以使用像 Michael Tokarev's inline qsort 这样的内联 qsort 实现。 This page 建议将性能提高 2 倍,这也是因为编译器能够内联比较函数。 使用 C++ std::sort。我知道您的代码是用 C 编写的,但您可以制作一个仅排序并提供 C 接口的小模块。您已经在使用具有出色 C++ 支持的工具链。 试试swenson/sort 的图书馆。它实现了许多算法,因此您可以使用最适合您的数据的算法。它似乎是可内联的,并且他们声称比 qsort 更快。 查找另一个排序库。可以做 Louis 的 Radix Sort 的东西是一个很好的建议。

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

【讨论】:

减去两个uint64_t 可能会生成一个不适合int 的值。如果不是这样,您可以完全消除分支并返回减法的结果。 >> 减去两个 uint64_t 可能会生成一个不适合 int 的值。我发现一半的减法产生了负面结果,这搞砸了整个排序。双步骤过程旨在回避这个问题。无论如何,我已经切换到 RADIX_SORT;快 4 倍! @BrianP007 减去无符号整数的另一个问题是它们不能为负数。 @BrianP007 如果您决定接受 Louis 的回答,您应该通过单击他的回答旁边的复选标记将其标记为“已接受”。【参考方案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));

【讨论】:

【参考方案4】:

也许一些?:而不是 ifs 会让事情变得更快。

【讨论】:

像这样的微优化通常是不好的——让编译器自己做吧。 这不太可能。任何构造得当的编译器都会为传统的if ... else 和相应的三元(即?:)表达式生成相同的代码。

以上是关于数百万 UINT64 RGBZ 图形像素的最快排序算法的主要内容,如果未能解决你的问题,请参考以下文章

对于具有数百万像素的 2D 未装箱像素阵列,建议使用哪种 Haskell 表示?

当您在 LAMP 服务器上拥有数百万用户时,存储和获取图像的最快和最有效的方法是啥?

具有数千、数百万和数十亿的问题排序数字(货币)

一百万个结构数组,根据其中一项值排序,用双链表还是数组排序效率更好,请给出最快C或C++算法代码。

在 uint64 中对 38 位进行基数排序?

在 OpenGL 中实例化数百万个对象:提高每秒帧数