霍尔的反驳和泡泡排序的卷土重来

2020-05-31 02:53:36

编者按:对于这篇博客,我欢迎我的朋友兼同事Gerben Stavenga作为客座作者。

最近,Andrei Alexandrescu发表了一篇关于使用Lomuto分区方案优化QuickSort的有趣文章。这篇文章的本质是,在许多情况下,快速排序的性能完全由分支错误预测主导,并且可以通过编写无分支代码来实现大幅加速。这已经被许多人观察到了,并且已经提出了各种无枝排序例程。Andrei观察到,从两个著名的快速排序分区方案来看,Lomuto很容易实现无分支,这对于对小的原语进行排序确实要好得多。我最近尝试了类似的想法,但把它们带到了一个不同但有趣的方向。我发现,Hoare和Lomuto的混合方案可以带来很大的改善,即使与无分支的Lomuto相比也是如此。最后令人惊讶的是,气泡排序在小数组中夺得桂冠。所有这些成功的关键是利用指令级并行性和减少依赖链。

快速排序是指对所有共享相同轮廓的数组进行排序的一类算法。

void QuickSort(T*Left,T*Right){if(Right-Left>;kCutOff){auto Pivot=ChoosePivotElement(Left,Right);//此处重要但不聚焦auto p=Partition(Pivot,Left,Right);//主工作循环QuickSort(Left,p);QuickSort(p,Right);//尾部调用,理想情况下是最大的子间隔}Else{SortSmallArray(Left,Right);}}。

在kCutOff的选择、用于小阵列的排序算法的选择和枢轴元素的选择方面存在着无数的变化。这些对于性能很重要,但是快速排序执行的主要工作是在Partition函数中完成的。有两种实现划分的规范方案:原始的Hoare方案和Lomuto方案。HoarePartition方案的工作原理是从数组前面将违反分区属性的元素与数组后面的元素交换,从外部向内处理数组,收敛到中间某处的分割点。

T*HoarePartition(T轴,T*Left,T*Right){While(Left<;Right){Left=ScanForward(Pivot,Left,Right);IF(Left==Right)Break;Right=ScanBackward(Pivot,Left,Right);If(Left==Right)Break;SWAP(*Left,*Right);}返回左侧;}。

相反,Lomuto方案从前到后处理数组,为到目前为止在每个步骤中处理的元素维护一个正确分区的数组。

T*LomutoPartition(T Pivot,T*Left,T*Right){T*p=Left;for(auto it=Left;it<;right;it++){IF(*it<;Pivot){std::SWAP(*it,*p);//可以是自交换p++;}}返回p;}

一个非常简单的循环,还有一个额外的属性,就是它对于小于轴心的元素的排序是稳定的。很容易看出,上面的循环可以在没有条件分支的情况下实现。条件交换可以通过条件移动来实现,并且条件指针增加可以作为无条件的p+=(*it;轴心)来实现。Andrei的博客展示了在各种标准库中,这个简单的无分支循环相对于产品质量实现的性能提升。

在这里,我想进一步进行优化,更深入地研究排序算法的性能分析。当条件分支被移除时,代码的性能趋向于变得更加稳定和更容易理解,数据依赖性成为关键。我们将分析一个自解释伪汇编代码的代码项,其中命名的变量应该被认为是显式的CPU寄存器和对内存的加载/存储。在此表示法中,上面的基本循环变成。

循环:val=加载(It)//1 prev_val=加载(P)//2 is_sirer=val<;PIVOT//3valnew=cmov(is_size,prev_val,val)//4 prev_val=cmov(is_size,val,prev_val)//5 Store(it,valnew)//6 Store(p,prev_val)//7 p+=is_size//。右)转到循环。

它将如何在并行的无序CPU上运行?作为第一近似值,我们可以假设CPU是任意并行的,即。只要指令是独立的,它就可以并行执行无限的指令。重要的是要了解携带依赖链的循环。它们决定了循环可能运行的最小延迟。在上面的代码中,您可以看到循环迭代之间的依赖关系由它携带,p。只有第8行和第9行参与循环携带的链,并且这两行都是单周期指令。因此,我们确定循环可能会以每次迭代1个周期的吞吐量运行。

然而,这种依赖分析并不十分正确。还有对内存的加载和存储。如果您将某些内容存储到某个内存地址,然后再从该地址加载,则必须读取前一个存储的值。这意味着加载依赖于存储,或者至少在地址重叠的情况下。这里,它和p是动态值,并且它们肯定可以重叠,如上图中的虚线所示。因此,让我们添加一个事实,即第1行和第2行的加载与第6行和第7行的存储之间存在依赖关系。

这完全改变了游戏,现在有了一个携带数据依赖性的长循环,第6行和第7行依赖于第4行和第5行,这两行都依赖于第3行,而第3行依赖于第1行和第2行的负载,而第1行和第2行可能依赖于前一次迭代第6行和第7行的存储。如果我们计算周期,我们得到5个周期用于加载(加载本身可以并行完成),1个周期用于比较,1个周期用于条件移动,1个周期用于存储,因此该循环将运行~8个周期。我们粗略的讨论表明,这与1个周期的迭代相去甚远。

虽然通常不可能对存储和加载进行重新排序,但这样做对CPU的性能至关重要。让我们来看一个简单的memcpy循环。

如果加载不能与前一个存储一起重新排序,则这是一个5+1=6周期延迟循环。然而,在memcpy中,可以保证加载和存储永远不会重叠。如果CPU重新排序,上面的循环将以每个周期一次迭代的吞吐量执行。它的执行应该是这样的,忽略了控制流所需的指令。

val_0=加载(It_0);it_1=it_0+1;//循环1 val_1=加载(It_1);it_2=it_1+1;//循环2 val_2=加载(It_2);it_3=it_2+1;//循环3 val_3=加载(It_3);it_4=it_3+1;//循环4//周期1的负载值变为可用。//从这一点开始,每个周期执行循环的所有指令。val_4=加载(It_4);it_5=it_4+1;存储(Val_0);//周期5 val_5=加载(It_5);it_6=it_5+1;存储(Val_1);//周期6。

实际上,指令流中加载之前的大多数存储实际上都指向不同的地址,并且可以对存储之前的加载重新排序。因此,CPU在存储之前对加载进行重新排序,这称为推测性加载,但是,只有在确认没有存储使推测性加载无效时,才会提交加载的影响。如果在前存储实际上使加载无效,则执行回滚到加载,并且CPU重新开始。人们可以想象,这是非常昂贵的,非常类似于分支机构的错误预测。虽然许多存储和加载位于不同的地址,但也有大量的存储和加载位于同一地址,例如,考虑到寄存器溢出。因此,CPU使用基于指令地址的预测模型来确定加载是否具有对先前存储的依赖性。一般来说,CPU是相当保守的,错误排序的成本非常高。在上面的代码中,CPU将遇到来自与最近存储相同的地址的加载,并且会犹豫是否进行必要的重新排序。

仔细看,p的加载才是有问题的。实际上,来自it的加载总是来自与前面存储不同的地址,此外,p处的加载还负责TheLoop中的大量额外工作。这是必要的,否则p处的存储将损坏数组中的值。被覆盖的值是先前在扫描中遇到的值,并且是大于透视的那些元素。相反,如果我们将这些值保存在临时缓冲区中,则不需要交换。

//将元素[LEFT,RIGHT]分布到//暂存缓冲区T*DistributeForward(T Pivot,T*Left,T*Right,T*ScratchSize-1;ptrdiff_t Offset=0;for(auto it=Left;it<;right;it++){auto val=*it;bool is_Larger=val>;=Pivot;auto dst=){auto val=*it;bool is_size=val>;;ptrdiff_t offset=0;for(auto it=Left;it<;right;it+){auto val=*it;bool is_size=val>;=voft;auto dst=。SCRATCH_END:it;dst[Offset]=val;Offset-=IS_LAGER;}RETURN RIGHT+OFFSET;}T*ModifiedLomutoPartition(T Pivot,T*Left,T*Right,T*Scratch){auto p=DistributeForward(Pivot,Left,Right,Scratch);//要完成分区,需要将Scratch中的元素复制//到数组末尾。auto size=right-p;memcpy(p,ScratchSize+kScratchSize-size,size*sizeof(T));返回p;}。

这是一个简单得多的循环,每个迭代只有一个加载和一个存储。更重要的是,加载永远不会与以前的存储冲突。这个循环比原来的循环快得多,在我的机器上,它不是每次迭代1个周期,而是2.5个周期。这表明它正在饱和CPU的ILP。不幸的是,上面的代码不再是适当的,它需要O(N)个额外的内存用于暂存缓冲区。

相反,如果我们使用较小的固定大小的临时缓冲区,我们仍然可以使用上面的代码,除非我们需要在固定缓冲区已满时中止。那我们该怎么办?函数返回离开循环的分割点p。此时,[Left,p]在数组前面有比Pivot更小的正确元素。暂存缓冲区充满了大于或等于Pivot的元素,[p,p+kScratchSize]包含我们不再需要的信息。我们的想法是,我们可以做相同的算法,但是向后,我们可以使用[p,p+kScratchSize)作为临时缓冲区。请注意DistributeForward()如何从后到前填充暂存缓冲区;向后版本将从前到后填充暂存缓冲区。因此,使用Interval[p,p+kScratchSize)作为临时变量执行DistributeBackwards()会将遇到的所有较小元素整齐地打包到正确的位置。这将一直持续到临时空间已满,但是现在在数组末尾打开了一个新的临时空间。等等,这看起来像是Hoare的算法,但它与Lomuto启发的分布式函数混合在一起。

T*ModifiedHoarePartition(T Pivot,T*Left,T*Right,T*Scratch){auto pleft=DistributeForward(Left,Right,Pivot,Scratch);if(right-pleft<;=kScratchSize){auto size=right-pleft;std::memcpy(pleft,sccratch+kScratchSize-size,size*sizeof(T));return pleft;}Left=pleft+kScratchSize-size,size*sizeof(T));返回pleft;}Left=pleft+kScratchSize。WHILE(TRUE){RIGHT=DistributeBackward(LEFT,RIGHT,

下一步是短数组的后备,其中递归的开销开始占据主导地位。在文献中,最常推荐的是插入排序和无数的微优化。在这一点上,我发现分区是如此之快,以至于快速排序比插入排序一直到只有几个元素。问题是插入排序有不可预测的分支,基本上每个插入有1个未命中。解决方案是冒泡排序。冒泡排序是一种非常可预测的访问模式,并且交换可以无分支地实现。再优化一点,您会发现您不需要交换。Onecan将最大值保存在寄存器中,并存储最小值。

void BubbleSort(T*arr,size_t n){for(size_t i=n;i>;1;i--){auto max=arr[0];for(size_t j=1;j<;i;j++){auto y=arr[j];arr[j-1]=(max<;=y?max:y);max=(max<;=y?y:max);}arr[i-1]=max;}}。

在上面的内循环中,max是参与最长循环承载的数据链的变量,并在其上进行比较和条件移动。这使得上面的循环在每次迭代中执行2个周期。不过,我们可以做得更好。不是让最大值冒泡,我们是在冒泡两个最大的元素,我们只需要重复冒泡阶段++n/2++次,而不是++n++次。事实证明,使用一个巧妙的实现,可以在同一个2个周期中冒泡两个元素。

void BubbleSort2(T*arr,size_t n){for(size_t i=n;i>;1;i-=2){auto x=arr[0];auto y=arr[1];if(y<;x)std::exchange(x,y);for(size_t j=2;j<;i;j++){auto z=arr[j];BOOL IS_SILEGER=y<;=z;AUTO w=IS_SIMPLETER?Y:Z;Y=IS_SIMPLETER?z:y;是否较小=x<;=z arr[j-2]=(是否较小?x:z);x=IS_SIMPLETER?w:x;}arr[i-2]=x;arr[i-1]=y;}}。

bm_SmallSort<;exp_gerbens::BubbleSort>;/2 2 ns 2 ns 400700000 bm_SmallSort<;exp_gerbens::BubbleSort>;/8 5 ns 5 ns 155200000 BM_SmallSort<;exp_gerbens::BubbleSort>;/32 14 ns 14 ns 52600000 BM_SmallSort<;/32 9 ns 9 ns 78400000 BM_SmallSort<;插入排序&>;/2 4 ns 4 ns 183600000 BM_SmallSort>;插入排序&>;/8 11 ns 11 ns 11 ns 63500000 BM_SmallSort<;插入排序&>;/32 17 ns 17 ns 17 ns 42000000。

事实上,可以对此进行推广,在CPU具有任意ILP的模型中,您可以在每次迭代的顶部++N++非常循环上冒泡。

结果是什么?结果令人叹为观止。下面是一个对100000个整数进行排序的简单基准。时间由元素的数量归一化,所以它是每个元素花费的时间。我在这里使用clang+libc++,因为GCC在发出分支免费代码方面明显更差。

处理器:采用超线程(36核)的英特尔Skylake Xeon DL1:32 KB DL2:1024KB DL3:24 MB基准测试时间(Ns)CPU(Ns)迭代--BM_SORT<;std_SORT>;51.651.6 10000000 BM_SORT<;std_STRATE_SORT>;65.65.6 10000000 BM_SORT<;lib_qsort>;90.490.5 7800000 BM_SORT<;Andrei_SORT>;32.632.6 21500000 BM_SORT<;exp_gerbens::QuickSort>;16.4 16.4 43200000

我们说的是复制自安德烈的实现的2倍的胜利。这些代码可以在my GitHub上找到,尽管它不包含由Andrei发布的代码的基准,因为它不包含许可证。

我们已经看到了解数据依赖关系对于优化代码有多么重要。特别是加载和存储之间的隐藏内存依赖关系会极大地影响工作循环的性能。理解代码的数据依赖图通常是真正的性能收益所在,但在博客圈子里很少有人关注它。我读过很多关于分支预测失误的影响、数据局部性和缓存的重要性的文章,但关于数据依赖性的文章要少得多。我敢打赌,像这样的问题“为什么链表很慢?”在局部性、缓存或不可预测的随机存储器访问方面得到了许多人的回答。至少我经常听到这些理由,就连斯特劳斯特鲁普也这么说。这些赛季可以起到一定作用,但不是主要原因。从根本上说,迭代链表在关键路径上有一个要使用的加载,这使得它比迭代平面数组慢5倍。此外,访问平面阵列允许循环展开,这可以进一步提高ILP。

这就带给我们一个答案,为什么快速排序比其他排序更快,理论复杂度更好,甚至更好。这都是关于数据依赖性的,上面的快速排序分区循环展示了一个独特的特性。它下一步要处理的元素不依赖于前几次迭代中的比较结果。将此与合并排序进行比较。在合并排序中,将比较两个头元素,但是需要比较的下一个元素取决于此比较的结果。实现无合并排序分支非常简单。它看起来就像。

val1=加载(LEFT+k-RIGHT_IDX)//1val2=加载(RIGHT+RIGHT_IDX)is_Smaller=val2<;val1//2tmp=cmov(IS_Smaller,val2,val1)存储(Out+k,tmp);k++;Right_IDX+=IS_Smaller//3。

这大约是每次迭代8个周期来更新RIGHT_IDX,我们在第1行有一个LOAD和非平凡索引(6个周期),第2行和第3行都是1个周期。类似的分析适用于堆排序,恢复堆属性要求比较最大的子级的子树上的两个子级和递归。

LEFT=加载(2*IDX)RIGHT=加载(2*IDX+1)IS_SILER=LEFT<;RIGHT TMP=cmov(IS_SIMPLER,RIGHT,LEFT)存储(IDX,TMP);IDX=2*IDX+IS_SNORER。

同样,这是环路承载的数据链上的8个周期。这就切中了为什么快速排序速度快的核心问题,尽管从理论上讲,与堆和合并排序相比,快速排序的行为较差。通过构造堆和合并排序,在隐含的树结构中非常均匀地划分数据,堆结构是最小深度的二叉树,合并排序执行的递归细分也是如此。这意味着他们进行的比较次数是++n\lg(N)++,紧贴着信息理论下限++\lg(n!)++。

相比之下,快速排序押注于以高概率获得合理平衡的树。从与枢轴的比较中提取的信息比特取决于它在数组中的排名。只有接近中值的枢轴才会导致每次比较获得1位信息。求解具有均匀枢轴选择的快速排序递归方程,给出了++2n\ln{n}++作为平均比较次数。因此,快速排序比堆排序或合并排序平均执行因子++2\ln(2)=1.4++更多的比较,或者相当于内部工作循环中更多的迭代。然而,基本工作回路中的巨大速度差异足以补偿这一信息理论因素。

此外,在对大型阵列进行分区时,花费少量工作改进旋转透视选项的中位数为3,可将此系数降至约1.2。此外,

..