固定长度为6的整型数组中最快的排序方法

418

回答另一个Stack Overflow的问题(this one),我偶然发现了一个有趣的子问题。如何最快地对包含6个整数的数组进行排序?

由于这个问题非常低级:

  • 我们不能假定库是可用的(而且调用本身也有成本),只能使用纯C。
  • 为了避免清空指令流水线(这是一种非常高的成本),我们应该尽可能减少分支、跳转和其他任何控制流中断(例如在&&||中隐藏的那些)。
  • 空间受限,最小化寄存器和内存使用是一个问题,最好采用原地排序。

实际上,这个问题是一种高尔夫游戏,目标不是使源长度最小化,而是执行时间最短。我称之为“禅”的代码,正如Michael Abrash的书《Zen of Code optimization》及其sequels中所使用的标题。

至于为什么它很有趣,有几个层次:

  • 这个例子简单易懂,容易衡量,不需要太多的C语言技巧
  • 它展示了选择一个好的算法对问题的影响,也展示了编译器和底层硬件的影响。

这是我的参考(天真的,未经优化的)实现和测试集。

#include <stdio.h>

static __inline__ int sort6(int * d){

    char j, i, imin;
    int tmp;
    for (j = 0 ; j < 5 ; j++){
        imin = j;
        for (i = j + 1; i < 6 ; i++){
            if (d[i] < d[imin]){
                imin = i;
            }
        }
        tmp = d[j];
        d[j] = d[imin];
        d[imin] = tmp;
    }
}

static __inline__ unsigned long long rdtsc(void)
{
  unsigned long long int x;
     __asm__ volatile (".byte 0x0f, 0x31" : "=A" (x));
     return x;
}

int main(int argc, char ** argv){
    int i;
    int d[6][5] = {
        {1, 2, 3, 4, 5, 6},
        {6, 5, 4, 3, 2, 1},
        {100, 2, 300, 4, 500, 6},
        {100, 2, 3, 4, 500, 6},
        {1, 200, 3, 4, 5, 600},
        {1, 1, 2, 1, 2, 1}
    };

    unsigned long long cycles = rdtsc();
    for (i = 0; i < 6 ; i++){
        sort6(d[i]);
        /*
         * printf("d%d : %d %d %d %d %d %d\n", i,
         *  d[i][0], d[i][6], d[i][7],
         *  d[i][8], d[i][9], d[i][10]);
        */
    }
    cycles = rdtsc() - cycles;
    printf("Time is %d\n", (unsigned)cycles);
}

原始结果

由于变量数量越来越多,我将它们全部收集到一个测试套件中,可以在这里找到。实际使用的测试比上面展示的要聪明一些,感谢Kevin Stock。您可以在自己的环境中编译和执行它。我对不同目标架构/编译器的行为非常感兴趣。(好的,伙计们,把它放在答案中,我会+1新结果集的每个贡献者)。

我一年前向Daniel Stutzbach(为了高尔夫)给出了答案,因为他当时是最快解决方案的源头(排序网络)。

Linux 64位,gcc 4.6.1 64位,Intel Core 2 Duo E8400,-O2

  • 直接调用qsort库函数:689.38
  • 朴素实现(插入排序):285.70
  • 插入排序(Daniel Stutzbach):142.12
  • 展开的插入排序:125.47
  • 排名排序:102.26
  • 寄存器排名排序:58.03
  • 排序网络(Daniel Stutzbach):111.68
  • 排序网络(Paul R):66.36
  • 带快速交换的排序网络12:58.86
  • 重排的排序网络12交换:53.74
  • 重排的简单交换排序网络12:31.54
  • 带快速交换的重排排序网络:31.54
  • 带快速交换的重排排序网络V2:33.63
  • 内联冒泡排序(Paolo Bonzini):48.85
  • 展开的插入排序(Paolo Bonzini):75.30

Linux 64位,gcc 4.6.1 64位,Intel Core 2 Duo E8400,-O1

  • 调用qsort库函数:705.93
  • 朴素实现(插入排序):135.60
  • 插入排序(Daniel Stutzbach):142.11
  • 展开的插入排序:126.75
  • 等级排序:46.42
  • 带寄存器的等级排序:43.58
  • 排序网络(Daniel Stutzbach):115.57
  • 排序网络(Paul R):64.44
  • 使用快速交换的12位排序网络:61.98
  • 重新排序的12位排序网络交换:54.67
  • 重新排序的12位排序网络简单交换:31.54
  • 带快速交换的重新排序排序网络:31.24
  • 带快速交换的重新排序排序网络V2:33.07
  • 内联冒泡排序(Paolo Bonzini):45.79
  • 展开的插入排序(Paolo Bonzini):80.15

我包含了-O1和-O2的结果,因为令人惊讶的是对于几个程序,O2比O1效率更低。我想知道具体哪种优化会产生这种影响?

关于提出的解决方案的评论

插入排序(Daniel Stutzbach)

如预期的那样,最小化分支确实是一个好主意。

排序网络(Daniel Stutzbach)

比插入排序更好。我想知道主要效果是否不是通过避免外部循环获得的。我尝试了展开插入排序以进行检查,确实我们得到了大致相同的数据(代码在这里)。
排序网络(Paul R)
到目前为止最好的。我用来测试的实际代码在这里。还不知道为什么它几乎比其他排序网络实现快两倍。参数传递?快速最大值?
带快速交换的12交换排序网络
正如Daniel Stutzbach建议的那样,我将他的12个交换排序网络与无分支快速交换组合在一起(代码在这里)。确实更快,迄今为止具有小幅度的优势(大约5%),因为使用1次较少的交换是可以预期的。
另外,有趣的是,分支交换似乎比在PPC架构上使用if的简单交换效率低得多(4倍)。
调用库qsort
为了再给一个参考点,我也试着按建议调用库函数qsort(代码在这里)。正如预期的那样,它慢得多:10到30倍……随着新的测试套件的出现,主要问题似乎是第一次调用后库的初始加载,而且与其他版本相比表现并不差。在我的Linux上,它只比其他版本慢3到20倍。在其他人测试时使用的某些架构上,它甚至似乎更快(我对这个结果感到非常惊讶,因为库函数qsort使用了更复杂的API)。 排名顺序 Rex Kerr提出了另一种完全不同的方法:针对数组的每个项直接计算其最终位置。这很有效,因为计算排名顺序不需要分支。但该方法的缺点是需要数组大小三倍的内存(一个数组副本和变量来存储排名顺序)。性能结果非常令人惊讶(也很有趣)。在我的32位操作系统和Intel Core2 Quad E8300参考架构上,循环计数略低于1000(像带分支交换的排序网络一样)。但是,当在我的64位系统(Intel Core2 Duo)上编译和执行时,它表现得更好:到目前为止,它成为最快的。我最终找到了真正的原因。我的32位系统使用gcc 4.4.1,而我的64位系统使用gcc 4.4.3,后者在优化这个特定代码方面要好得多(对其他提案几乎没有什么区别)。 更新: 如上面发布的数字表明,这个效果仍然受到后来版本的gcc的增强,排名顺序变得比任何其他替代方案都快两倍。

重新排序的交换排序网络12

Rex Kerr提议使用gcc 4.4.3实现的惊人效率让我想知道:一个内存使用量是无分支排序网络三倍的程序怎么可能更快呢?我的假设是它具有更少的读写依赖关系,从而允许更好地利用x86的超标量指令调度器。这给了我一个想法:重新排序交换以最小化读写依赖关系。简单来说,当你执行SWAP(1,2); SWAP(0,2);时,必须等待第一个交换完成后才能执行第二个,因为两个都访问同一内存单元。当你执行SWAP(1,2); SWAP(4,5);时,处理器可以并行执行两个操作。我尝试了一下,效果符合预期,排序网络运行速度提高了约10%。

简单交换排序网络12

在原始帖子发表一年后,Steinar H. Gunderson建议我们不要试图聪明地超越编译器,保持交换代码的简单性。这确实是个好主意,因为结果代码大约快了40%!他还提出了一种优化手动使用x86内联汇编代码的交换方法,可以节省更多的时间。最令人惊讶的是(它表明了程序员心理的很多东西),一年前我们都没有尝试过那个版本的交换。我用来测试的代码在这里。其他人建议以其他方式编写C快速交换,但使用合适的编译器时,它们与简单的方法产生相同的性能。

现在,“最佳”代码如下:

static inline void sort6_sorting_network_simple_swap(int * d){
#define min(x, y) (x<y?x:y)
#define max(x, y) (x<y?y:x) 
#define SWAP(x,y) { const int a = min(d[x], d[y]); \
                    const int b = max(d[x], d[y]); \
                    d[x] = a; d[y] = b; }
    SWAP(1, 2);
    SWAP(4, 5);
    SWAP(0, 2);
    SWAP(3, 5);
    SWAP(0, 1);
    SWAP(3, 4);
    SWAP(1, 4);
    SWAP(0, 3);
    SWAP(2, 5);
    SWAP(1, 3);
    SWAP(2, 4);
    SWAP(2, 3);
#undef SWAP
#undef min
#undef max
}

如果我们相信我们的测试集(是的,它相当差,它的唯一好处是简短、简单并且易于理解我们正在测量的内容),则一个排序的结果代码的平均循环次数低于40个周期(执行了6个测试)。这使得每个交换的平均时间为4个周期。我认为这非常快。还有其他改进可能吗?

2
你对整数有一些限制吗?例如,我们可以假设任何两个x、y的x-yx+y不会导致下溢或上溢吗? - Matthieu M.
3
你应该尝试将我的12次交换排序网络与Paul的无分支交换函数结合起来。他的解决方案将所有参数作为独立元素传递到堆栈上,而不是传递一个指向数组的单个指针。这可能也会有所不同。 - Daniel Stutzbach
2
请注意,在64位上正确实现rdtsc的方法是__asm__ volatile (".byte 0x0f, 0x31; shlq $32, %%rdx; orq %%rdx, %0" : "=a" (x) : : "rdx");,因为rdtsc将答案放在EDX:EAX中,而GCC希望它在一个64位寄存器中。您可以通过在-O3下编译来查看错误。另请参见我对Paul R关于更快速的SWAP的评论。 - Paolo Bonzini
3
@Tyler: 在汇编语言中,如果没有使用分支指令,你如何实现它? - Loren Pechtel
4
@Loren说:CMP EAX, EBX; SBB EAX, EAX会根据EAXEBX的大小关系,在EAX中放入0或0xFFFFFFFF。 SBB是"带借位减法",与ADC("带进位加法")相对应; 您所指的状态位 就是 进位位。然而,我记得在 Pentium 4 上,ADCSBB的延迟和吞吐量很差,与ADDSUB相比仍然慢两倍,并且在Core CPU上也是如此。自从80386以来,还有SETcc条件存储和CMOVcc条件移动指令,但它们也很慢。 - j_random_hacker
显示剩余37条评论
25个回答

166

任何优化都需要反复测试。我建议至少尝试排序网络和插入排序。基于过去的经验,我会把钱押在插入排序上。

您了解输入数据吗?某些算法在特定类型的数据上表现更好。例如,插入排序在排序或几乎排序的数据上表现更好,因此如果几乎排好序的数据出现的概率高于平均水平,则它将是更好的选择。

您发布的算法类似于插入排序,但看起来您最小化了交换次数,代价是进行更多的比较。不过,与交换相比,比较要昂贵得多,因为分支可能导致指令流水线停滞。

下面是一个插入排序实现:

static __inline__ int sort6(int *d){
        int i, j;
        for (i = 1; i < 6; i++) {
                int tmp = d[i];
                for (j = i; j >= 1 && tmp < d[j-1]; j--)
                        d[j] = d[j-1];
                d[j] = tmp;
        }
}

这是我建立排序网络的方法。首先,使用这个网站来生成适当长度的网络所需的最小SWAP宏集。将其封装在一个函数中:

static __inline__ int sort6(int * d){
#define SWAP(x,y) if (d[y] < d[x]) { int tmp = d[x]; d[x] = d[y]; d[y] = tmp; }
    SWAP(1, 2);
    SWAP(0, 2);
    SWAP(0, 1);
    SWAP(4, 5);
    SWAP(3, 5);
    SWAP(3, 4);
    SWAP(0, 3);
    SWAP(1, 4);
    SWAP(2, 5);
    SWAP(2, 4);
    SWAP(1, 3);
    SWAP(2, 3);
#undef SWAP
}

9
+1:很好,你只用了12个步骤完成了任务,而不是我手工编码和经验推导网络中的13个步骤。如果可以的话,我会给你另外一个+1,因为你提供了一个能够生成网络的网址,我已经将其添加到书签中了。 - Paul R
9
如果您期望大多数请求都是小型数组,那么这是一个通用排序函数的绝妙想法。对于您想要优化的情况,请使用 switch 语句,按照此过程进行操作;让默认情况使用库排序函数。 - Mark Ransom
6
一种好的库排序函数通常会针对小数组使用快速路径。许多现代库会使用递归的快速排序或归并排序,在递归到n < SMALL_CONSTANT后切换到插入排序。 - Daniel Stutzbach
3
@Mark,一个C语言库排序函数需要通过函数指针指定比较操作。对于每次比较都调用函数的开销很大。通常情况下,这仍然是最简洁的方法,因为这很少成为程序中的关键路径。但是,如果这是关键路径,我们确实可以更快地完成排序,如果我们知道我们正在对整数进行排序,并且恰好有6个数。 :) - Daniel Stutzbach
8
XOR 交换几乎总是不明智的做法。 - Paul R
显示剩余6条评论

66

这里是使用排序网络实现的代码:

inline void Sort2(int *p0, int *p1)
{
    const int temp = min(*p0, *p1);
    *p1 = max(*p0, *p1);
    *p0 = temp;
}

inline void Sort3(int *p0, int *p1, int *p2)
{
    Sort2(p0, p1);
    Sort2(p1, p2);
    Sort2(p0, p1);
}

inline void Sort4(int *p0, int *p1, int *p2, int *p3)
{
    Sort2(p0, p1);
    Sort2(p2, p3);
    Sort2(p0, p2);  
    Sort2(p1, p3);  
    Sort2(p1, p2);  
}

inline void Sort6(int *p0, int *p1, int *p2, int *p3, int *p4, int *p5)
{
    Sort3(p0, p1, p2);
    Sort3(p3, p4, p5);
    Sort2(p0, p3);  
    Sort2(p2, p5);  
    Sort4(p1, p2, p3, p4);  
}

针对这个问题,你需要非常高效的无分支的 minmax 实现,因为本质上这段代码可以抽象成一系列的 minmax 操作(总共有13个)。我将这留给读者作为练习。

请注意,这个实现方法很容易适用于向量化(例如,SIMD - 大多数 SIMD 指令集都有向量 min/max 操作),也很容易实现于 GPU 上(例如,CUDA - 由于是无分支的,所以不存在翻转等问题)。

同时参见:快速算法实现用于排序非常小的列表


1
有关最小值/最大值的一些位操作技巧请参见:http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax - Rubys
1
@Paul:在实际的CUDA使用环境中,这肯定是最好的答案。我会检查它在高尔夫x64环境中是否也是(以及有多少),并公布结果。 - kriss
1
如果您注意到(a+b+c)-(min+max)是中间数字,那么Sort3在大多数架构上会更快。 - Rex Kerr
1
@Rex:我看明白了 - 看起来不错。对于像AltiVec和SSE这样的SIMD架构,指令周期数是相同的(最大值和最小值都是单周期指令,如加/减),但对于普通标量CPU,您的方法看起来更好。 - Paul R
2
如果我让GCC使用条件移动指令来优化min函数,我可以获得33%的速度提升:#define SWAP(x,y) { int dx = d[x], dy = d[y], tmp; tmp = d[x] = dx < dy ? dx : dy; d[y] ^= dx ^ tmp; }。这里我没有在d[y]中使用?:是因为它的性能稍微差一些,但几乎可以忽略不计。 - Paolo Bonzini
显示剩余2条评论

49

既然这些是整数,比较速度很快,为什么不直接计算每个的排名顺序:

inline void sort6(int *d) {
  int e[6];
  memcpy(e,d,6*sizeof(int));
  int o0 = (d[0]>d[1])+(d[0]>d[2])+(d[0]>d[3])+(d[0]>d[4])+(d[0]>d[5]);
  int o1 = (d[1]>=d[0])+(d[1]>d[2])+(d[1]>d[3])+(d[1]>d[4])+(d[1]>d[5]);
  int o2 = (d[2]>=d[0])+(d[2]>=d[1])+(d[2]>d[3])+(d[2]>d[4])+(d[2]>d[5]);
  int o3 = (d[3]>=d[0])+(d[3]>=d[1])+(d[3]>=d[2])+(d[3]>d[4])+(d[3]>d[5]);
  int o4 = (d[4]>=d[0])+(d[4]>=d[1])+(d[4]>=d[2])+(d[4]>=d[3])+(d[4]>d[5]);
  int o5 = 15-(o0+o1+o2+o3+o4);
  d[o0]=e[0]; d[o1]=e[1]; d[o2]=e[2]; d[o3]=e[3]; d[o4]=e[4]; d[o5]=e[5];
}

1
@Rex:抱歉,我第一眼漏掉了> vs >=模式。它适用于所有情况。 - kriss
3
@kriss: 哦,这并不完全令人惊讶——有很多变量在浮动,它们必须被精心排序并缓存在寄存器中等等。 - Rex Kerr
这种方法的效率在很大程度上取决于在特定机器上执行(int)(a>=b)有多容易。Pentiums提供了“setcc”指令,允许用两个指令完成此操作,但在某些处理器上,可能需要更复杂的代码。虽然在几乎所有现代设备上都可以不使用分支实现,但编译器并不总是遵从。如果您可以假设减法没有溢出,则((a-b)>>31)会给出0或-1,并且可以代替使用。 - greggo
这里的一个重要优点是主要部分没有“存储器”,因此优化器可以轻松地在寄存器中工作,并且仅读取每个d[i]一次。如果编译器能够确定,例如e[0]已经被读取为d[0](并且您有足够的寄存器),那就更好了。这可以通过使用6个本地变量d0..d5而不是将memcpy用于e来明确表示。 - greggo
2
@SSpoke 0+1+2+3+4+5=15 其中一个数字缺失,用15减去其余数字的和即可得到缺失的数字。 - Glenn Teitelbaum
显示剩余8条评论

35

看起来我来晚了一年,但是我们还是开始吧...

通过观察由gcc 4.5.2生成的汇编代码,我发现每次交换时都会进行加载和存储,这实际上是不必要的。最好将6个值加载到寄存器中,对它们进行排序,然后再将它们存回内存。我将加载和存储的顺序调整为最接近需要第一次使用和最后一次使用寄存器的位置。我还使用了Steinar H. Gunderson的SWAP宏。更新:我转而使用Paolo Bonzini的SWAP宏,gcc将其转换为类似于Gunderson的内容,但gcc能够更好地排序指令,因为它们没有作为显式汇编给出。

我使用与重新排序的交换网络给出的最佳性能相同的交换顺序,尽管可能有更好的顺序。如果我有更多时间,我会生成和测试一堆排列。

我将测试代码改为考虑超过4000个数组,并显示每个数组所需的平均周期数。在i5-650上(使用-O3),我得到约34.1个周期/排序,而原始的重新排序排序网络(使用-O1,胜过-O2和-O3)则得到约65.3个周期/排序。

#include <stdio.h>

static inline void sort6_fast(int * d) {
#define SWAP(x,y) { int dx = x, dy = y, tmp; tmp = x = dx < dy ? dx : dy; y ^= dx ^ tmp; }
    register int x0,x1,x2,x3,x4,x5;
    x1 = d[1];
    x2 = d[2];
    SWAP(x1, x2);
    x4 = d[4];
    x5 = d[5];
    SWAP(x4, x5);
    x0 = d[0];
    SWAP(x0, x2);
    x3 = d[3];
    SWAP(x3, x5);
    SWAP(x0, x1);
    SWAP(x3, x4);
    SWAP(x1, x4);
    SWAP(x0, x3);
    d[0] = x0;
    SWAP(x2, x5);
    d[5] = x5;
    SWAP(x1, x3);
    d[1] = x1;
    SWAP(x2, x4);
    d[4] = x4;
    SWAP(x2, x3);
    d[2] = x2;
    d[3] = x3;

#undef SWAP
#undef min
#undef max
}

static __inline__ unsigned long long rdtsc(void)
{
    unsigned long long int x;
    __asm__ volatile ("rdtsc; shlq $32, %%rdx; orq %%rdx, %0" : "=a" (x) : : "rdx");
    return x;
}

void ran_fill(int n, int *a) {
    static int seed = 76521;
    while (n--) *a++ = (seed = seed *1812433253 + 12345);
}

#define NTESTS 4096
int main() {
    int i;
    int d[6*NTESTS];
    ran_fill(6*NTESTS, d);

    unsigned long long cycles = rdtsc();
    for (i = 0; i < 6*NTESTS ; i+=6) {
        sort6_fast(d+i);
    }
    cycles = rdtsc() - cycles;
    printf("Time is %.2lf\n", (double)cycles/(double)NTESTS);

    for (i = 0; i < 6*NTESTS ; i+=6) {
        if (d[i+0] > d[i+1] || d[i+1] > d[i+2] || d[i+2] > d[i+3] || d[i+3] > d[i+4] || d[i+4] > d[i+5])
            printf("d%d : %d %d %d %d %d %d\n", i,
                    d[i+0], d[i+1], d[i+2],
                    d[i+3], d[i+4], d[i+5]);
    }
    return 0;
}
我更改了测试套件,以便报告每个排序和运行的时钟,并运行更多的测试(也更新了cmp函数以处理整数溢出),这里是在一些不同架构上的结果。我尝试在 AMD CPU 上进行测试,但我可用的 X6 1100T 上 rdtsc 不可靠。修改了测试套件
Clarkdale (i5-650)
==================
Direct call to qsort library function      635.14   575.65   581.61   577.76   521.12
Naive implementation (insertion sort)      538.30   135.36   134.89   240.62   101.23
Insertion Sort (Daniel Stutzbach)          424.48   159.85   160.76   152.01   151.92
Insertion Sort Unrolled                    339.16   125.16   125.81   129.93   123.16
Rank Order                                 184.34   106.58   54.74    93.24    94.09
Rank Order with registers                  127.45   104.65   53.79    98.05    97.95
Sorting Networks (Daniel Stutzbach)        269.77   130.56   128.15   126.70   127.30
Sorting Networks (Paul R)                  551.64   103.20   64.57    73.68    73.51
Sorting Networks 12 with Fast Swap         321.74   61.61    63.90    67.92    67.76
Sorting Networks 12 reordered Swap         318.75   60.69    65.90    70.25    70.06
Reordered Sorting Network w/ fast swap     145.91   34.17    32.66    32.22    32.18

Kentsfield (Core 2 Quad)
========================
Direct call to qsort library function      870.01   736.39   723.39   725.48   721.85
Naive implementation (insertion sort)      503.67   174.09   182.13   284.41   191.10
Insertion Sort (Daniel Stutzbach)          345.32   152.84   157.67   151.23   150.96
Insertion Sort Unrolled                    316.20   133.03   129.86   118.96   105.06
Rank Order                                 164.37   138.32   46.29    99.87    99.81
Rank Order with registers                  115.44   116.02   44.04    116.04   116.03
Sorting Networks (Daniel Stutzbach)        230.35   114.31   119.15   110.51   111.45
Sorting Networks (Paul R)                  498.94   77.24    63.98    62.17    65.67
Sorting Networks 12 with Fast Swap         315.98   59.41    58.36    60.29    55.15
Sorting Networks 12 reordered Swap         307.67   55.78    51.48    51.67    50.74
Reordered Sorting Network w/ fast swap     149.68   31.46    30.91    31.54    31.58

Sandy Bridge (i7-2600k)
=======================
Direct call to qsort library function      559.97   451.88   464.84   491.35   458.11
Naive implementation (insertion sort)      341.15   160.26   160.45   154.40   106.54
Insertion Sort (Daniel Stutzbach)          284.17   136.74   132.69   123.85   121.77
Insertion Sort Unrolled                    239.40   110.49   114.81   110.79   117.30
Rank Order                                 114.24   76.42    45.31    36.96    36.73
Rank Order with registers                  105.09   32.31    48.54    32.51    33.29
Sorting Networks (Daniel Stutzbach)        210.56   115.68   116.69   107.05   124.08
Sorting Networks (Paul R)                  364.03   66.02    61.64    45.70    44.19
Sorting Networks 12 with Fast Swap         246.97   41.36    59.03    41.66    38.98
Sorting Networks 12 reordered Swap         235.39   38.84    47.36    38.61    37.29
Reordered Sorting Network w/ fast swap     115.58   27.23    27.75    27.25    26.54

Nehalem (Xeon E5640)
====================
Direct call to qsort library function      911.62   890.88   681.80   876.03   872.89
Naive implementation (insertion sort)      457.69   236.87   127.68   388.74   175.28
Insertion Sort (Daniel Stutzbach)          317.89   279.74   147.78   247.97   245.09
Insertion Sort Unrolled                    259.63   220.60   116.55   221.66   212.93
Rank Order                                 140.62   197.04   52.10    163.66   153.63
Rank Order with registers                  84.83    96.78    50.93    109.96   54.73
Sorting Networks (Daniel Stutzbach)        214.59   220.94   118.68   120.60   116.09
Sorting Networks (Paul R)                  459.17   163.76   56.40    61.83    58.69
Sorting Networks 12 with Fast Swap         284.58   95.01    50.66    53.19    55.47
Sorting Networks 12 reordered Swap         281.20   96.72    44.15    56.38    54.57
Reordered Sorting Network w/ fast swap     128.34   50.87    26.87    27.91    28.02

1
@cdunn2001,我刚刚测试了一下,除了在-O0和-Os模式下有几个周期的改进之外,我没有看到任何改进。从汇编代码来看,gcc似乎已经成功地使用寄存器并消除了对memcpy的调用。 - Kevin Stock
1
你的代码仍然使用 Gunderson 的交换算法,而我的则是 #define SWAP(x,y) { int oldx = x; x = x < y ? x : y; y ^= oldx ^ x; } - Paolo Bonzini
在我看来,装载的顺序很可能会导致停顿。由于所有的装载都发生在第一次使用之前,交换必须等待变量被加载,从而导致停顿。交错的存储是可以的,你应该希望在数据存储和下一个交换发生之间有一些重叠。无论如何,对于一个足够好的编译器来说,这可能都是一个无关紧要的问题,它会重新排序事物,以便在等待从RAM获取数据时发生工作。在GPU上,将所有的装载都放在一个块中可能会非常有益。 - Axman6
有六个元素的排列方式恰好有720种。从某种意义上说,“最佳”排序应该是在这720种排序中具有最佳的最坏情况性能的排序。 - Yakk - Adam Nevraumont
@cdunn2001 “Rank Order”解决方案通常会在寄存器中,因为它会推迟存储操作,直到完成所有主要计算。使用6个本地变量进行读取而不是使用memcpy到“e”可能会提高性能。我不会感到惊讶如果gcc可以自己做到这一点(它*可以将小型本地数组视为单独的int变量; 在这里,它还需要将memcpy转换为6个单独的操作e [0] = d [0]; e [1] = d [1]..以进行优化。但最好在代码中明确表示)。 - greggo
显示剩余7条评论

15
我几天前从谷歌上发现了这个问题,因为我也需要快速排序一个长度固定的6个整数数组。然而,在我的情况下,我的整数只有8位(而不是32位),并且我没有严格要求仅使用C语言。无论如何,我想分享一下我的发现,以防它们可能对某人有用...
我实现了一种变体的网络排序,使用SSE向量化比较和交换操作,尽可能地进行优化。完全排序该数组需要六个“passes”。我使用了一种新颖的机制,直接将PCMPGTB(向量化比较)的结果转换为PSHUFB(向量化交换)的洗牌参数,仅使用PADDB(向量化加法)和在某些情况下还有PAND(按位与)指令。
这种方法还具有一个副作用,即产生了一个真正的无分支函数。没有任何跳转指令。
看起来,这个实现比当前标记为最快选项的实现(“Sorting Networks 12 with Simple Swap”)快约38%。在我的测试中,我修改了那个实现,使用char数组元素,以使比较公平。
我应该注意到,这种方法可以应用于任何最多16个元素的数组大小。我预计相对于其他替代方案的速度优势会随着数组大小的增加而变得更大。
代码是使用SSSE3的x86_64处理器的MASM编写的。该函数使用“新”的Windows x64调用约定。这就是它...
PUBLIC simd_sort_6

.DATA

ALIGN 16

pass1_shuffle   OWORD   0F0E0D0C0B0A09080706040503010200h
pass1_add       OWORD   0F0E0D0C0B0A09080706050503020200h
pass2_shuffle   OWORD   0F0E0D0C0B0A09080706030405000102h
pass2_and       OWORD   00000000000000000000FE00FEFE00FEh
pass2_add       OWORD   0F0E0D0C0B0A09080706050405020102h
pass3_shuffle   OWORD   0F0E0D0C0B0A09080706020304050001h
pass3_and       OWORD   00000000000000000000FDFFFFFDFFFFh
pass3_add       OWORD   0F0E0D0C0B0A09080706050404050101h
pass4_shuffle   OWORD   0F0E0D0C0B0A09080706050100020403h
pass4_and       OWORD   0000000000000000000000FDFD00FDFDh
pass4_add       OWORD   0F0E0D0C0B0A09080706050403020403h
pass5_shuffle   OWORD   0F0E0D0C0B0A09080706050201040300h
pass5_and       OWORD 0000000000000000000000FEFEFEFE00h
pass5_add       OWORD   0F0E0D0C0B0A09080706050403040300h
pass6_shuffle   OWORD   0F0E0D0C0B0A09080706050402030100h
pass6_add       OWORD   0F0E0D0C0B0A09080706050403030100h

.CODE

simd_sort_6 PROC FRAME

    .endprolog

    ; pxor xmm4, xmm4
    ; pinsrd xmm4, dword ptr [rcx], 0
    ; pinsrb xmm4, byte ptr [rcx + 4], 4
    ; pinsrb xmm4, byte ptr [rcx + 5], 5
    ; The benchmarked 38% faster mentioned in the text was with the above slower sequence that tied up the shuffle port longer.  Same on extract
    ; avoiding pins/extrb also means we don't need SSE 4.1, but SSSE3 CPUs without SSE4.1 (e.g. Conroe/Merom) have slow pshufb.
    movd    xmm4, dword ptr [rcx]
    pinsrw  xmm4,  word ptr [rcx + 4], 2  ; word 2 = bytes 4 and 5


    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass1_shuffle]
    pcmpgtb xmm5, xmm4
    paddb xmm5, oword ptr [pass1_add]
    pshufb xmm4, xmm5

    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass2_shuffle]
    pcmpgtb xmm5, xmm4
    pand xmm5, oword ptr [pass2_and]
    paddb xmm5, oword ptr [pass2_add]
    pshufb xmm4, xmm5

    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass3_shuffle]
    pcmpgtb xmm5, xmm4
    pand xmm5, oword ptr [pass3_and]
    paddb xmm5, oword ptr [pass3_add]
    pshufb xmm4, xmm5

    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass4_shuffle]
    pcmpgtb xmm5, xmm4
    pand xmm5, oword ptr [pass4_and]
    paddb xmm5, oword ptr [pass4_add]
    pshufb xmm4, xmm5

    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass5_shuffle]
    pcmpgtb xmm5, xmm4
    pand xmm5, oword ptr [pass5_and]
    paddb xmm5, oword ptr [pass5_add]
    pshufb xmm4, xmm5

    movdqa xmm5, xmm4
    pshufb xmm5, oword ptr [pass6_shuffle]
    pcmpgtb xmm5, xmm4
    paddb xmm5, oword ptr [pass6_add]
    pshufb xmm4, xmm5

    ;pextrd dword ptr [rcx], xmm4, 0    ; benchmarked with this
    ;pextrb byte ptr [rcx + 4], xmm4, 4 ; slower version
    ;pextrb byte ptr [rcx + 5], xmm4, 5
    movd   dword ptr [rcx], xmm4
    pextrw  word ptr [rcx + 4], xmm4, 2  ; x86 is little-endian, so this is the right order

    ret

simd_sort_6 ENDP

END

您可以将此编译为可执行对象并将其链接到您的C项目中。有关在Visual Studio中如何执行此操作的说明,请阅读本文。 您可以使用以下C原型从您的C代码调用该函数:

void simd_sort_6(char *values);

将您的汇编级别提案与其他提案进行比较会很有趣。实现的性能比较不包括它们。无论如何,使用SSE听起来不错。 - kriss
未来研究的另一个领域将是将新的Intel AVX指令应用于此问题。更大的256位向量足够大,可以容纳8个DWORDs。 - Joe Crivello
1
不要使用 pxor / pinsrd xmm4, mem, 0,直接使用 movd - Peter Cordes

15
测试代码非常糟糕;它溢出了初始数组(这里的人们不阅读编译器警告吗?),printf 打印的是错误的元素,使用 .byte 进行 rdtsc 没有好的理由,只有一次运行(!),没有检查最终结果是否正确(因此很容易将其“优化”成一些微妙错误的东西),包含的测试非常简单(没有负数?)并且没有阻止编译器将整个函数作为无用代码进行丢弃。

话虽如此,改进比特网络解决方案也很容易;只需将 min/max/SWAP 更改为

#define SWAP(x,y) { int tmp; asm("mov %0, %2 ; cmp %1, %0 ; cmovg %1, %0 ; cmovg %2, %1" : "=r" (d[x]), "=r" (d[y]), "=r" (tmp) : "0" (d[x]), "1" (d[y]) : "cc"); }

对我来说,它大约快了65%(使用Debian gcc 4.4.5和-O2、amd64、Core i7)。


好的,测试代码很差。请随意改进它。是的,您可以使用汇编代码。为什么不走到底,完全使用x86汇编器编码呢?这可能会稍微不太便携,但为什么要费心呢? - kriss
感谢您注意到数组溢出问题,我已经进行了更正。其他人可能没有注意到这个问题,因为他们点击链接复制/粘贴代码时,并没有发现溢出情况。 - kriss
4
实际上,甚至不需要汇编器;如果你放弃所有聪明的技巧,GCC会识别这个序列并为你插入条件移动:#define min(a, b) ((a < b) ? a : b) #define max(a, b) ((a < b) ? b : a) #define SWAP(x,y) { int a = min(d[x], d[y]); int b = max(d[x], d[y]); d[x] = a; d[y] = b; }这可能比内联汇编慢几个百分点,但由于缺乏适当的基准测试,很难说。 - Steinar H. Gunderson
3
最后,如果你的数字是浮点数,也不用担心NaN等问题,GCC可以将其转换为minss / maxss SSE指令,这样速度会快约25%。趣味:放弃聪明的位操作技巧,让编译器发挥作用。 :-) - Steinar H. Gunderson

13

虽然我很喜欢提供的swap宏:

#define min(x, y) (y ^ ((x ^ y) & -(x < y)))
#define max(x, y) (x ^ ((x ^ y) & -(x < y)))
#define SWAP(x,y) { int tmp = min(d[x], d[y]); d[y] = max(d[x], d[y]); d[x] = tmp; }

我看到了一些改进(好的编译器可能会做到):

#define SWAP(x,y) { int tmp = ((x ^ y) & -(y < x)); y ^= tmp; x ^= tmp; }

我们注意到min和max函数是如何工作的,并明确提取公共子表达式。这将完全消除min和max宏。


这会让它们反过来,注意到 d[y] 得到了最大值,即 x^(公共子表达式)。 - Kevin Stock
我注意到了同样的事情;我认为为了使你的实现正确,你想要使用 d[x] 而不是 x(对于 y 也是一样),并且在这里的不等式中应该是 d[y] < d[x](是的,与最小/最大代码不同)。 - Tyler
我试过用你的交换方法,但是在更大的级别上进行局部优化会产生负面影响(我猜它会引入依赖)。而且结果比另一个交换方法还要慢。但是,正如你所看到的,通过优化交换,的确可以获得更高的性能。 - kriss

12

在优化min/max之前,请务必进行基准测试并查看实际生成的编译器汇编代码。如果我允许GCC使用条件移动指令来优化min,我可以获得33%的加速:

#define SWAP(x,y) { int dx = d[x], dy = d[y], tmp; tmp = d[x] = dx < dy ? dx : dy; d[y] ^= dx ^ tmp; }

在测试代码中,使用 ?: 进行最大值比较与使用上述方法相差不大,但上述方法稍微更快一些(280 vs. 420 cycles)。这种 SWAP 方法在 GCC 和 Clang 中都更快。

编译器在寄存器分配和别名分析方面也做得非常出色,有效地将 d[x] 移入本地变量中,并且只在最后才复制回内存。事实上,它们甚至比完全使用本地变量(如 d0 = d[0], d1 = d[1], d2 = d[2], d3 = d[3], d4 = d[4], d5 = d[5])还要好。我写这篇文章是因为你假设了强大的优化,却试图在 min/max 上超越编译器。:)

顺便说一下,我尝试过 Clang 和 GCC。它们都进行了相同的优化,但由于调度差异,两者的结果有所不同,不能真正确定哪个更快或更慢。GCC 在排序网络上更快,Clang 在二次排序上更快。

为了完整起见,展开的冒泡排序和插入排序也是可能的。这里是冒泡排序:

SWAP(0,1); SWAP(1,2); SWAP(2,3); SWAP(3,4); SWAP(4,5);
SWAP(0,1); SWAP(1,2); SWAP(2,3); SWAP(3,4);
SWAP(0,1); SWAP(1,2); SWAP(2,3);
SWAP(0,1); SWAP(1,2);
SWAP(0,1);

这里是插入排序算法:

//#define ITER(x) { if (t < d[x]) { d[x+1] = d[x]; d[x] = t; } }
//Faster on x86, probably slower on ARM or similar:
#define ITER(x) { d[x+1] ^= t < d[x] ? d[x] ^ d[x+1] : 0; d[x] = t < d[x] ? t : d[x]; }
static inline void sort6_insertion_sort_unrolled_v2(int * d){
    int t;
    t = d[1]; ITER(0);
    t = d[2]; ITER(1); ITER(0);
    t = d[3]; ITER(2); ITER(1); ITER(0);
    t = d[4]; ITER(3); ITER(2); ITER(1); ITER(0);
    t = d[5]; ITER(4); ITER(3); ITER(2); ITER(1); ITER(0);

这个插入排序比Daniel Stutzbach的更快,并且在GPU或具有预测功能的计算机上效果特别好,因为ITER可以仅使用3条指令(与SWAP的4条相比)。例如,以下是ARM汇编中的代码“t = d [2]; ITER(1); ITER(0);”:
    MOV    r6, r2
    CMP    r6, r1
    MOVLT  r2, r1
    MOVLT  r1, r6
    CMP    r6, r0
    MOVLT  r1, r0
    MOVLT  r0, r6

对于六个元素,插入排序与排序网络竞争(12个交换次数与15个迭代平衡4个指令/交换次数与3个指令/迭代);冒泡排序当然更慢。但是当大小增长时,情况将不再如此,因为插入排序是O(n^2),而排序网络是O(n log n)。


1
更或者少相关:我向GCC提交了一份报告,以便它可以直接在编译器中实现优化。不确定是否会完成,但至少您可以跟踪其发展。 - Morwenn

11

我将测试套件移植到了一台无法识别的 PPC 架构的机器上(无需修改代码,只需增加测试迭代次数,使用8个测试用例避免使用模块污染结果,并替换x86特定的rdtsc):

直接调用 qsort 库函数:101

朴素实现(插入排序):299

插入排序(Daniel Stutzbach):108

展开的插入排序:51

排序网络(Daniel Stutzbach):26

排序网络(Paul R):85

具有快速交换的12位排序网络:117

重新排列交换的12位排序网络:116

排名排序:56


1
非常有趣。在 PPC 上,无分支交换似乎是一个不好的想法。这也可能是与编译器相关的影响。使用了哪一个? - kriss
这是gcc编译器的一个分支 - 最小值、最大值逻辑可能不是无分支的 - 我会检查反汇编并让您知道,但除非编译器足够聪明,包括像x < y这样的内容而不需要if语句仍然会成为一个分支 - 在x86 / x64上,CMOV指令可能会避免这种情况,但在PPC上没有固定点值的此类指令,只有浮点数。我可能会在明天尝试一下并让您知道 - 我记得Winamp AVS源代码中有一个更简单的无分支min / max,但如果我没记错,它只适用于浮点数 - 但可能是朝着真正无分支方法的良好开端。 - jheriko
4
下面是一个用于PPC无符号输入的无分支min/max函数:subfc r5,r4,r3; subfe r6,r6,r6; andc r6,r5,r6; add r4,r6,r4; subf r3,r6,r3。其中r3/r4是输入寄存器,r5/r6是临时寄存器,输出结果中r3得到最小值,r4得到最大值。这个函数可以手动优化调度。我利用GNU超级优化工具从4条指令的min/max序列开始寻找可组合的两个序列。对于有符号输入,你可以在开始时将所有元素加上0x80000000,在结束时再减去,然后按照无符号数进行处理。 - Paolo Bonzini

7

在您的交换函数中,使用XOR交换可能会很有用。

void xorSwap (int *x, int *y) {
     if (*x != *y) {
         *x ^= *y;
         *y ^= *x;
         *x ^= *y;
     }
 }

如果你保证所有的整数都是唯一的,那么 if 可能会导致你的代码出现太多分歧,但这样做可能会很方便。

1
异或交换也适用于相等的值... x^=y将x设置为0,y^=x将y保持为y(==x),x^=y将x设置为y。 - jheriko
11
xy指向相同位置时,它就不起作用 - hobbs
无论如何,当与排序网络一起使用时,我们永远不会同时调用指向相同位置的x和y。仍然需要找到一种方法来避免测试哪一个更大以获得无分支交换的效果。我有一个实现这一点的想法。 - kriss

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接