std::sortよりも高速なC++ Quicksort

blqsort - std::sortよりも高速なC++ Quicksort

blqsort は非常に高速なクイックソートアルゴリズムです。ランダムデータに対して、これより高速な汎用比較ベースソートアルゴリズムは他に知りません。小規模なプリミティブデータ型を扱う場合のみ、SIMDベースのアルゴリズムが一般的に高速です。

実装はMITライセンスの下でGitHubで利用可能です。C言語版やマルチスレッド版もあり、それぞれ単一のヘッダーファイルとして提供されています。

blqsort on Github

---

以下のベンチマークは、5000万個のランダムなdoubleをソートする際の実行時間を示しています。測定は、Clang を使用した Apple M1 システムと、GCC を使用した Intel i5-9300H Linux システムで行われ、どちらも -O3 オプションでコンパイルされました。

---

実世界のデータはしばしば部分的にソートされています。ここでは90%ソート済みの入力に対するベンチマーク結果を示します。

---

ランダムなid値を持つこれらのstructを5000万個ソートする際の実行時間です。

struct entry {
    int id;
    int value;
    bool operator<(const entry& other) const { return id < other.id; }
};

公平な比較のために、ここではblqsortのシングルスレッドバージョンが使用されました。M1では、マルチスレッドバージョンはさらに約3倍高速です。実行時間に関して、CバージョンはC++バージョンとほとんど違いがありません。

既にソート済みおよび逆ソート済みの入力は検出され、効率的に処理されます。多くの重複キーを持つ入力も最適化されています。

blqsortの主要な概念は、分岐なしのパーティショニング、小規模パーティションのためのソートネットワーク、および優れたピボット選択です。

基本的なクイックソートからblqsortへ

実行時間は、特に指定がない限り、Apple M1 システムに適用されます。

A small, elegant, recursive quicksort:

#include <cstdio>
#include <cstdlib>
#include <chrono>

constexpr int SIZE = 50000000;
double data[SIZE];

double time() {
    static auto time0 = std::chrono::high_resolution_clock::now();
    auto t = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff = t - time0;
    return diff.count();
}

void sortq(double ar[], int len) {
    if (len < 2) return;
    double piv = ar[0];
    int mid = 0;
    for (int i = 1; i < len; i++) {
        if (ar[i] < piv) {
            mid += 1;
            std::swap(ar[mid], ar[i]);
        }
    }
    std::swap(ar[mid], ar[0]);

    sortq(ar + mid + 1, len - mid - 1);
    sortq(ar, mid);
}

int main() {
    double t0;

    printf("Sorting %d million doubles\n", SIZE / 1000000);
    for (int i = 0; i < SIZE; i++) data[i] = rand() / 1024.0;

    t0 = time();
    sortq(data, SIZE);
    printf("sortq: %.2fs\n", time() - t0);
}

c++ -O1 main.cpp && ./a.out

Sorting 50 million doubles
sortq: 4.61s

ソート済みのデータを再度ソートすると、パーティションのすべての要素がピボット要素以上になります。2番目のパーティションは毎回1つずつ小さくなり、関数が数百万回再帰的に呼び出されるため、数時間後にスタックオーバーフローが発生します。

これをテールコール最適化で回避できます。つまり、最初の再帰呼び出しはそのままにし、2番目の呼び出しをループイテレーションに置き換えます。この場合、小さい方のパーティションが再帰的に呼び出されることを保証する必要があります。したがって、必要に応じてそれらをスワップします。

void sortq(double* ar, int len) {
    while (len >= 2) {
        double piv = ar[0];
        int mid = 0;
        for (int i = 1; i < len; i++) {
            if (ar[i] < piv) {
                mid += 1;
                std::swap(ar[mid], ar[i]);
            }
        }
        std::swap(ar[mid], ar[0]);
        double* ar2 = ar + mid + 1;
        int len2 = len - mid - 1;
        len = mid;
        if (len2 > len) {
            std::swap(len, len2);
            std::swap(ar, ar2);
        }
        sortq(ar2, len2);
    }
}

これでプログラムはクラッシュしなくなりましたが、完了するのに数日かかります。この不都合なパーティショニングのため、O(n²) の最悪ケースシナリオが現実に起こります。ただし、挿入ソートでもこれより速くなることはありません。

中央の要素をピボットとして使用するだけで、ソート済みデータをソートする問題は解消されます。

    .
    while (len >= 2) {
        std::swap(ar[len / 2], ar[0]);
        double piv = ar[0];
    .

残るリスクは、誰かが意図的に特別に細工されたデータを送信してこの動作を引き起こすことです。この脅威を排除するために、ピボットポイントはランダムに選択されるべきです。

この実装は多くの目的に十分であり、また比較的迅速です。ちなみに、ここでは「Lomuto」パーティショニングが使用されています。「Hoare」パーティショニングはより少ない比較とスワップを行いますが、最新のCPUではキャッシュの動作が不利になる可能性があります。

std::sort が存在するのに、なぜ独自の С++ 実装を書くのでしょうか？

.
#include <algorithm>
.
    .
    t0 = time();
    std::sort(data, data + SIZE);
    printf("std::sort: %.2fs\n", time() - t0);
    .

std::sort: 1.31s

これはARM Macでの結果です。x86システムでは、std::sort は約10%しか速くありません。

クイックソートをstd::sortより高速に実行するにはどうすればよいでしょうか？

ソートネットワーク

約20要素までのパーティションでは、挿入ソートはクイックソートより高速です。しかし、小規模パーティションにはさらに高速なものがあります。それはソートネットワークです。各パーティションサイズは独自のプログラムブランチで処理されます。2つの要素のパーティションサイズの場合、2つの要素を比較し、必要に応じてスワップするだけです。しかし、要素数が増えるにつれて、ロジックはますます複雑になります。

inline void sort2(double& a, double& b) {
    if (a < b) std::swap(a, b);
}
inline void sort3(double& a, double& b, double& c) {
    sort2(a, b); sort2(b, c); sort2(a, b);
}
inline void sort4(double& a, double& b, double& c, double& d) {
    sort2(a, b); sort2(c, d); sort2(a, c);
    sort2(b, d); sort2(b, c);
}
inline void sort5(double& a, double& b, double& c, double& d, double& e) {
    sort2(b, c); sort2(d, e); sort2(b, d);
    sort2(a, c); sort2(a, d); sort2(c, e);
    sort2(a, b); sort2(c, d); sort2(b, c);
}

    .
    while (len >= 6) {
        .
        .
    }
    if (len == 5) sort5(ar[0], ar[1], ar[2], ar[3], ar[4]);
    else if (len == 4) sort4(ar[0], ar[1], ar[2], ar[3]);
    else if (len == 3) sort3(ar[0], ar[1], ar[2]);
    else if (len == 2) sort2(ar[0], ar[1]);
}

if が遅くする場合、それを避ける

簡単なタスクを考えてみましょう。配列を反復処理し、500未満のすべての数値を新しい配列にコピーします。C での直接的な実装は次のようになります。

smlen = 0;
for (int i = 0; i < 1000; i++) {
    if (numbers[i] < 500) {
        small_numbers[smlen] = numbers[i];
        smlen += 1;
    }
}

数値がランダムに分布している場合、if 条件の結果は分岐予測にとって予測不可能になります。これにより、予測ミス率が高くなり、プロセッサが常にパイプラインをフラッシュして実行を再開する必要があるため、パフォーマンスに深刻な影響を与えます。

分岐なしのアプローチを使用することで、この問題を回避できます。

smlen = 0;
for (int i = 0; i < 1000; i++) {
    small_numbers[smlen] = numbers[i];
    smlen += (numbers[i] < 500);
}

このバージョンでは、条件が満たされない場合でも無条件に書き込み操作（メモリへの書き込み）が発生しますが、関連するオーバーヘッドは、分岐予測ミスによるパフォーマンス低下よりも一般的に大幅に低くなります。

分岐なしは、ジャンプが全くないことを意味するわけではありません。この用語は、予測が困難な比較のみを指します。明らかに、ループジャンプは残りますが、予測が容易であるため問題はありません。1回だけスキップされます。

パフォーマンス比較

意味のある結果を得るために、100,000個のランダムな数値を使用して両方のルーチンを1000回実行します。

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define SIZE 100000
int numbers[SIZE];
int small_numbers[SIZE];
int smlen;

double ts(void) {
    static double t0;
    struct timeval tv;
    gettimeofday(&tv, NULL);
    double h = t0;
    t0 = tv.tv_sec + tv.tv_usec / 1000000.0;
    return t0 - h;
}
void init(void) {
    srand(1);
    for (int i = 0; i < SIZE; i++) numbers[i] = rand() % 1000;
    ts();
}
void small_if(void) {
    smlen = 0;
    for (int i = 0; i < SIZE; i++) {
        if (numbers[i] < 500) {
            small_numbers[smlen] = numbers[i];
            smlen += 1;
        }
    }
}
void small_bl(void) {
    smlen = 0;
    for (int i = 0; i < SIZE; i++) {
        small_numbers[smlen] = numbers[i];
        smlen += (numbers[i] < 500);
    }
}
int main(void) {
    init();
    for (int i = 0; i < 1000; i++) small_if();
    printf("Time if: %.3f\n", ts());
    init();
    for (int i = 0; i < 1000; i++) small_bl();
    printf("Time bl: %.3f\n", ts());
}

ベンチマーク結果は、Apple M1 で -O1 最適化を使用して取得されました。このベンチマークでは、-O3 は追加の速度向上をもたらさず、生成コードが大きくなるだけでした。

Time if: 0.345
Time bl: 0.036

アセンブリ解析

objdump -d を使用すると、CPUがこれらのループを処理する方法の違いを確認できます。

if バージョンでは、コンパイラは条件付き分岐 b.gt を生成します。

ldr     w13, [x10], #0x4           ; load numbers[i]
cmp     w13, #0x1f3                ; compare with 499 (0x1f3)
b.gt    0x100000678                ; branch if greater (The bottleneck)
str     w13, [x12, w8, sxtw #2]    ; store value in small_numbers[j]
add     w8, w8, #0x1               ; j += 1

分岐なしバージョンでは、コンパイラは cinc（条件付きインクリメント）命令を使用します。

ldr     w13, [x10], #0x4           ; load numbers[i]
str     w13, [x12, w8, uxtw #2]    ; unconditional store to small_numbers[j]
cmp     w13, #0x1f4                ; compare with 500
cinc    w8, w8, lt                 ; conditional increment (No branch)

このバージョンでは、データに基づいたジャンプはありません。制御フローは完全に線形であり、CPUはパイプラインの全幅を使用できます。

コンパイラはなぜ if ケースを自動的に最適化しないのか？

コンパイラはデータの性質を知りません。もし数値のわずかな割合だけが500未満だった場合、予測子はほぼ常に正しいので、分岐バージョンの方が実際に高速になる可能性があります。

分岐なしバージョンは、毎回書き込み操作（str）を実行します。コンパイラは、元のソースコードの書き込み操作が特定の条件でのみ発生するように意図されていた場合（例えば、割り当てられたバッファの末尾を超えて書き込むのを避けるため）、書き込みを許可できると安全に仮定できません。

クイックソートは基本的にデータを大小に分割することによって機能するため、このテクニックを使用して高速化できるはずです。

分岐なしクイックソート（明示的なバージョン）

したがって、前述のように、パーティショニングステップを分岐なしにすることができます。（comp は比較関数です）：

int h = comp(*left, piv);
*lwr = *rwr = *left++;
lwr += h; rwr -= !h;

また、ソートネットワークの基本操作である sort2 を分岐なしにすることもできます。

inline static void sort2(double& a, double& b) {
    double x = a; double y = b;
    bool m = x < y;
    a = m ? x : y; b = m ? y : x;
}

仕組みは（要するに）次のとおりです。

まず、約250個の要素が1つのページから補助バッファにコピーされ、後続の操作のためのスペースが確保されます。次に、ジャンプを使用せずに、ブロックを交互に左と右にコピーします。

ピボット要素も慎重に選択されます。パーティションのさまざまな要素から複数のメディアンが決定され、最終的なメディアンがこれらのメディアンから計算されます。これにより、データの分布がはるかに均一になります。一貫して最悪ケースの動作を引き起こす入力を防ぐために、ピボット選択プロセスにいくつかのランダム性が追加されます。

ただし、分布がまだ非常に不均一な場合、アルゴリズムは予防措置としてこれらのパーティションにヒープソートに切り替わります。不均一な分布は、多くの同一要素がある場合にしばしば発生します。この場合、それらは識別され、グループ化されます。

さらに、コードにはソートネットワークが含まれており、これは最大16のパーティションサイズに使用されます。

そして for (int i = 16; i--;) { ... } ループで、コンパイラはループを展開し、ループ本体の16個のコピーを生成します。

このバージョンは単一のC++ヘッダーファイルで構成されており、簡単にインクルードできます。std::sort と同様に、任意のデータ型をソートでき、比較関数をユーザーが提供できます。

blqs.h

#ifndef BLQS_H
#define BLQS_H

#include <cstddef>
#include <cstdint>
#include <cstring>
#include <functional>
#include <utility>

#ifndef BLQS_RAND
#define BLQS_RAND

#include <random>
static thread_local uint64_t blqs_seed;

static inline unsigned blqs_rand(void) {
    if (!blqs_seed) {
        std::random_device rd;
        blqs_seed =
            ((uintptr_t)&blqs_seed) ^
            (((uint64_t)rd() << 32) | (uint64_t)rd());
        if (!blqs_seed) blqs_seed = 1;
    }

    uint64_t x = blqs_seed;
    x ^= x << 13;
    x ^= x >> 7;
    x ^= x << 17;
    blqs_seed = x;
    return (unsigned)x;
}
#endif

namespace blqs {

constexpr int SMALLPART = 1024;
constexpr int SWSZ = 512;
constexpr int UNROLL = 16;

template<typename T, typename Compare>
static inline void sort2(T& a, T& b, Compare comp) {
    T x = a; T y = b;
    bool m = comp(x, y);
    a = m ? x : y; b = m ? y : x;
}

template<typename T, typename Compare>
static inline void sort3(T& a, T& b, T& c, Compare comp) {
    sort2(a, b, comp); sort2(b, c, comp); sort2(a, b, comp);
}

template<typename T, typename Compare>
static inline void sort4(T& a, T& b, T& c, T& d, Compare comp) {
    sort2(a, b, comp); sort2(c, d, comp); sort2(a, c, comp);
    sort2(b, d, comp); sort2(b, c, comp);
}

template<typename T, typename Compare>
static inline void sort5(T& a, T& b, T& c, T& d, T& e, Compare comp) {
    sort2(b, c, comp); sort2(d, e, comp); sort2(b, d, comp);
    sort2(a, c, comp); sort2(a, d, comp); sort2(c, e, comp);
    sort2(a, b, comp); sort2(c, d, comp); sort2(b, c, comp);
}

template<typename T, typename Compare>
static inline void sort6(T& a, T& b, T& c, T& d, T& e, T& f, Compare comp) {
    sort2(a, b, comp); sort2(c, d, comp); sort2(e, f, comp);
    sort2(a, c, comp); sort2(b, d, comp); sort2(e, f, comp);
    sort2(a, e, comp); sort2(b, f, comp); sort2(c, e, comp);
    sort2(d, f, comp); sort2(b, c, comp); sort2(d, e, comp);
    sort2(c, d, comp);
}

template<typename T, typename Compare>
static inline void sort7(T& a, T& b, T& c, T& d, T& e, T& f, T& g, Compare comp) {
    sort2(a, g, comp); sort2(c, d, comp); sort2(e, f, comp);
    sort2(a, c, comp); sort2(b, e, comp); sort2(d, g, comp);
    sort2(a, b, comp); sort2(c, f, comp); sort2(d, e, comp);
    sort2(b, c, comp); sort2(e, g, comp);
    sort2(c, d, comp); sort2(e, f, comp);
    sort2(b, c, comp); sort2(d, e, comp); sort2(f, g, comp);
}

template<typename T, typename Compare>
static inline void sort8(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, Compare comp) {
    sort2(a,b,comp); sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp);
    sort2(a,c,comp); sort2(b,d,comp); sort2(e,g,comp); sort2(f,h,comp);
    sort2(b,c,comp); sort2(f,g,comp);
    sort2(a,e,comp); sort2(b,f,comp); sort2(c,g,comp); sort2(d,h,comp);
    sort2(c,e,comp); sort2(d,f,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,g,comp);
}

template<typename T, typename Compare>
static inline void sort9(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, Compare comp) {
    sort2(a,d,comp); sort2(b,h,comp); sort2(c,f,comp); sort2(e,i,comp);
    sort2(a,h,comp); sort2(c,e,comp); sort2(d,i,comp); sort2(f,g,comp);
    sort2(a,c,comp); sort2(b,d,comp); sort2(e,f,comp); sort2(h,i,comp);
    sort2(b,e,comp); sort2(d,g,comp); sort2(f,h,comp);
    sort2(a,b,comp); sort2(c,e,comp); sort2(d,f,comp); sort2(g,i,comp);
    sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,g,comp);
}

template<typename T, typename Compare>
static inline void sort10(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, Compare comp) {
    sort2(a,b,comp); sort2(c,f,comp); sort2(d,g,comp); sort2(e,h,comp); sort2(i,j,comp);
    sort2(a,g,comp); sort2(b,i,comp); sort2(c,e,comp); sort2(d,j,comp); sort2(f,h,comp);
    sort2(a,c,comp); sort2(b,d,comp); sort2(e,f,comp); sort2(g,i,comp); sort2(h,j,comp);
    sort2(a,b,comp); sort2(c,h,comp); sort2(d,f,comp); sort2(e,g,comp); sort2(i,j,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp);
    sort2(b,d,comp); sort2(c,e,comp); sort2(f,h,comp); sort2(g,i,comp);
    sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp);
}

template<typename T, typename Compare>
static inline void sort11(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k, Compare comp) {
    sort2(a,j,comp); sort2(b,g,comp); sort2(c,e,comp); sort2(d,h,comp); sort2(f,i,comp);
    sort2(a,b,comp); sort2(d,f,comp); sort2(e,k,comp); sort2(g,j,comp); sort2(h,i,comp);
    sort2(b,d,comp); sort2(c,f,comp); sort2(e,h,comp); sort2(i,k,comp);
    sort2(a,e,comp); sort2(b,c,comp); sort2(d,h,comp); sort2(f,j,comp); sort2(g,i,comp);
    sort2(a,b,comp); sort2(c,g,comp); sort2(e,f,comp); sort2(h,i,comp);
    sort2(j,k,comp);
    sort2(c,e,comp); sort2(d,g,comp); sort2(f,h,comp); sort2(i,j,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp);
    sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp);
}

template<typename T, typename Compare>
static inline void sort12(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k, T& l, Compare comp) {
    sort2(a,i,comp); sort2(b,h,comp); sort2(c,g,comp); sort2(d,l,comp); sort2(e,k,comp); sort2(f,j,comp);
    sort2(a,c,comp); sort2(b,e,comp); sort2(d,f,comp); sort2(g,i,comp); sort2(h,k,comp); sort2(j,l,comp);
    sort2(a,b,comp); sort2(c,j,comp); sort2(e,h,comp); sort2(f,g,comp); sort2(k,l,comp);
    sort2(b,d,comp); sort2(c,h,comp); sort2(e,j,comp); sort2(i,k,comp);
    sort2(a,b,comp); sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp); sort2(i,j,comp); sort2(k,l,comp);
    sort2(b,c,comp); sort2(d,f,comp); sort2(g,i,comp); sort2(j,k,comp);
    sort2(c,e,comp); sort2(d,g,comp); sort2(f,i,comp); sort2(h,j,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp); sort2(j,k,comp);
}
template<typename T, typename Compare>
static inline void sort13(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k, T& l, T& m,
                          Compare comp) {
    sort2(a,m,comp); sort2(b,k,comp); sort2(c,j,comp); sort2(d,h,comp); sort2(f,l,comp); sort2(g,i,comp);
    sort2(b,g,comp); sort2(c,d,comp); sort2(e,l,comp); sort2(h,j,comp); sort2(i,k,comp);
    sort2(a,e,comp); sort2(b,c,comp); sort2(d,g,comp); sort2(h,i,comp); sort2(j,k,comp); sort2(l,m,comp);
    sort2(e,g,comp); sort2(f,j,comp); sort2(i,l,comp); sort2(k,m,comp);
    sort2(a,f,comp); sort2(d,i,comp); sort2(e,h,comp); sort2(g,l,comp); sort2(j,k,comp);
    sort2(a,b,comp); sort2(c,f,comp); sort2(g,j,comp); sort2(h,i,comp); sort2(k,l,comp);
    sort2(b,d,comp); sort2(c,e,comp); sort2(f,g,comp); sort2(j,k,comp);
    sort2(b,c,comp); sort2(d,e,comp); sort2(f,h,comp); sort2(g,i,comp);
    sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp); sort2(i,j,comp);
    sort2(d,e,comp); sort2(f,g,comp);
}
template<typename T, typename Compare>
static inline void sort14(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k, T& l,
                        T& m, T& n, Compare comp) {
    sort2(a,b,comp); sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp); sort2(i,j,comp); sort2(k,l,comp); sort2(m,n,comp);
    sort2(a,c,comp); sort2(b,d,comp); sort2(e,i,comp); sort2(f,j,comp); sort2(k,m,comp); sort2(l,n,comp);
    sort2(a,e,comp); sort2(b,c,comp); sort2(d,h,comp); sort2(f,i,comp); sort2(g,k,comp); sort2(j,n,comp); sort2(l,m,comp);
    sort2(a,g,comp); sort2(b,f,comp); sort2(d,j,comp); sort2(e,k,comp); sort2(h,n,comp); sort2(i,m,comp);
    sort2(c,k,comp); sort2(d,l,comp); sort2(e,g,comp); sort2(h,j,comp);
    sort2(b,d,comp); sort2(c,i,comp); sort2(f,l,comp); sort2(g,h,comp); sort2(k,m,comp);
    sort2(b,e,comp); sort2(c,g,comp); sort2(d,f,comp); sort2(h,l,comp); sort2(i,k,comp); sort2(j,m,comp);
    sort2(c,e,comp); sort2(d,g,comp); sort2(f,i,comp); sort2(h,k,comp); sort2(j,l,comp);
    sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp); sort2(j,k,comp);
    sort2(g,h,comp);
}
template<typename T, typename Compare>
static inline void sort15(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k, T& l,
                        T& m, T& n, T& o, Compare comp) {
    sort2(b,c,comp); sort2(d,k,comp); sort2(e,o,comp); sort2(f,i,comp); sort2(g,n,comp); sort2(h,m,comp); sort2(j,l,comp);
    sort2(a,o,comp); sort2(b,f,comp); sort2(c,i,comp); sort2(d,h,comp); sort2(g,j,comp); sort2(k,m,comp); sort2(l,n,comp);
    sort2(a,h,comp); sort2(b,g,comp); sort2(c,j,comp); sort2(e,k,comp); sort2(f,l,comp); sort2(i,n,comp); sort2(m,o,comp);
    sort2(a,g,comp); sort2(c,e,comp); sort2(d,f,comp); sort2(h,l,comp); sort2(i,k,comp); sort2(j,m,comp); sort2(n,o,comp);
    sort2(a,d,comp); sort2(b,c,comp); sort2(e,h,comp); sort2(f,j,comp); sort2(g,i,comp); sort2(k,l,comp); sort2(m,n,comp);
    sort2(a,b,comp); sort2(c,d,comp); sort2(e,g,comp); sort2(h,j,comp); sort2(k,m,comp); sort2(l,n,comp);
    sort2(b,c,comp); sort2(d,f,comp); sort2(i,k,comp); sort2(l,m,comp);
    sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp); sort2(j,k,comp);
    sort2(c,d,comp); sort2(e,f,comp); sort2(g,h,comp); sort2(i,j,comp); sort2(k,l,comp);
    sort2(f,g,comp); sort2(h,i,comp);
}
template<typename T, typename Compare>
static inline void sort16(T& a, T& b, T& c, T& d, T& e, T& f, T& g, T& h, T& i, T& j, T& k,
                            T& l, T& m, T& n, T& o, T& p, Compare comp) {
    sort2(a,n,comp); sort2(b,m,comp); sort2(c,p,comp); sort2(d,o,comp); sort2(e,i,comp);
    sort2(f,g,comp); sort2(h,l,comp); sort2(j,k,comp);
    sort2(a,f,comp); sort2(b,h,comp); sort2(c,j,comp); sort2(d,e,comp); sort2(g,n,comp);
    sort2(i,o,comp); sort2(k,p,comp); sort2(l,m,comp);
    sort2(a,b,comp); sort2(c,d,comp); sort2(e,f,comp); sort2(g,i,comp); sort2(h,j,comp);
    sort2(k,l,comp); sort2(m,n,comp); sort2(o,p,comp);
    sort2(a,c,comp); sort2(b,d,comp); sort2(e,k,comp); sort2(f,l,comp); sort2(g,h,comp);
    sort2(i,j,comp); sort2(m,o,comp); sort2(n,p,comp);
    sort2(b,c,comp); sort2(d,m,comp); sort2(e,g,comp); sort2(f,h,comp); sort2(i,k,comp);
    sort2(j,l,comp); sort2(n,o,comp);
    sort2(b,e,comp); sort2(c,g,comp); sort2(f,i,comp); sort2(h,k,comp); sort2(j,n,comp);
    sort2(l,o,comp);
    sort2(c,e,comp); sort2(d,g,comp); sort2(j,m,comp); sort2(l,n,comp);
    sort2(d,f,comp); sort2(g,i,comp); sort2(h,j,comp); sort2(k,m,comp);
    sort2(d,e,comp); sort2(f,g,comp); sort2(h,i,comp); sort2(j,k,comp); sort2(l,m,comp);
    sort2(g,h,comp); sort2(i,j,comp);
}

template<typename T, typename Compare>
static inline void sorting_network(T* l, int partsz_min1, Compare comp) {
    switch (partsz_min1) {
    case 0: break;
    case 1:  sort2(l[0], l[1], comp); break;
    case 2:  sort3(l[0], l[1], l[2], comp); break;
    case 3:  sort4(l[0], l[1], l[2], l[3], comp); break;
    case 4:  sort5(l[0], l[1], l[2], l[3], l[4], comp); break;
    case 5:  sort6(l[0], l[1], l[2], l[3], l[4], l[5], comp); break;
    case 6:  sort7(l[0], l[1], l[2], l[3], l[4], l[5], l[6], comp); break;
    case 7:  sort8(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], comp); break;
    case 8:  sort9(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], comp); break;
    case 9:  sort10(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
         l[8], l[9], comp); break;
    case 10: sort11(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], comp); break;
    case 11: sort12(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], l[11], comp); break;
    case 12: sort13(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], l[11], l[12], comp); break;
    case 13: sort14(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], l[11], l[12], l[13], comp); break;
    case 14: sort15(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], l[11], l[12], l[13], l[14], comp); break;
    case 15: sort16(l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7],
        l[8], l[9], l[10], l[11], l[12], l[13], l[14], l[15], comp); break;
    default: break;
    }
}

template<typename T, typename Compare>
static inline void heap_sift_down(T* left, ptrdiff_t n, ptrdiff_t j, T& k, Compare& comp) {
    while (j * 2 + 1 < n) {
        ptrdiff_t child = j * 2 + 1;
        if (child + 1 < n && comp(left[child], left[child + 1])) child++;
        if (!comp(k, left[child])) break;
        left[j] = std::move(left[child]);
        j = child;
    }
    left[j] = std::move(k);
}

template<typename T, typename Compare>
void heap_sort(T* left, T* right, Compare comp) {
    ptrdiff_t n = right - left + 1;
    if (n < 2) return;

    for (ptrdiff_t i = n / 2; ; ) {
        if (i > 0) {
            i--;
            T k = std::move(left[i]);
            heap_sift_down(left, n, i, k, comp);
        }
        else {
            n--;
            if (n == 0) return;
            T k = std::move(left[n]);
            left[n] = std::move(left[0]);
            heap_sift_down(left, n, 0, k, comp);
        }
    }
}
template<typename T, typename Compare>
static inline void med5(T& a, T& b, T& c, T& d, T& e, Compare comp) {
    sort2(a, b, comp); sort2(c, d, comp);
    sort2(a, c, comp); sort2(b, d, comp);
    sort2(b, c, comp); sort2(c, e, comp);
    sort2(b, c, comp);
}

template<typename T, typename Compare>
static inline void med7(T& a, T& b, T& c, T& d, T& e, T& f, T& g, Compare comp) {
    sort2(a, b, comp); sort2(c, d, comp); sort2(a, c, comp);
    sort2(b, d, comp); sort2(b, c, comp); sort2(e, f, comp);
    sort2(e, g, comp); sort2(a, e, comp);
    sort2(b, f, comp); sort2(c, g, comp); sort2(b, e, comp);
    sort2(c, e, comp);
    sort2(d, f, comp); sort2(d, e, comp);
}

template<typename T, typename Compare>
static inline T* partition_small(T* left, T* right, Compare comp) {

    T* outerleft = left;
    T* pivp = left + (right - left) / 2;

    T piv = *pivp;
    T l1 = left[1], l2 = left[2], l3 = left[3];
    T r2 = right[-2], r1 = right[-1], r0 = *right;
    med7(l1, l2, l3, piv, r2, r1, r0, comp);
    left[1] = l1; left[2] = l2; left[3] = l3;
    right[-2] = r2; right[-1] = r1; *right = r0;
    left += 4;
    right -= 3;

    *pivp = *outerleft;
    T swbuf[SMALLPART];
    T *sw = swbuf, *lwr = left;

    while (right - left >= UNROLL) for (int i = UNROLL; i--;) {
        int h = comp(*left, piv);
        *lwr = *sw = *left++;
        lwr += h; sw += !h;
    }
    while (left <= right) {
        int h = comp(*left, piv);
        *lwr = *sw = *left++;
        lwr += h; sw += !h;
    }
    std::memcpy(lwr, swbuf, (sw - swbuf) * sizeof(T));
    lwr -= 1;
    *outerleft = *lwr;
    *lwr = piv;
    return lwr;
}

template<typename T, typename Compare>
static T* partition_large(T* left, T* right, Compare comp) {

    T* outerleft = left;
    T* pivp = left + (right - left) / 2;
    T piv = *pivp;

    int j = blqs_rand() & 7;
    med5(left[1 + j], left[2 + j], left[3 + j], left[4 + j], left[5 + j], comp);
    med5(left[11 + j], left[12 + j], left[13 + j], left[14 + j], left[15 + j], comp);
    med5(left[21 + j], left[22 + j], left[23 + j], left[24 + j], left[25 + j], comp);
    med5(pivp[-2], pivp[-1], piv, pivp[1], pivp[2], comp);
    med5(right[-24 - j], right[-23 - j], right[-22 - j], right[-21 - j], right[-20 - j], comp);
    med5(right[-14 - j], right[-13 - j], right[-12 - j], right[-11 - j], right[-10 - j], comp);
    med5(right[-4 - j], right[-3 - j], right[-2 - j], right[-1 - j], right[ 0 - j], comp);

    med7(left[3 + j], left[13 + j], left[23 + j], piv,  right[-22 - j], right[-12 - j],
        right[-2 - j] , comp);

    left += 1;
    *pivp = *outerleft;

    T swbuf[SWSZ];
    T *lwr = left, *rwr = right, *sw = swbuf;

    while (UNROLL < SWSZ - (sw - swbuf) && left < right - UNROLL) {
        ptrdiff_t avail = std::min(right - left, SWSZ - (sw - swbuf));
        T* endp = right - avail;
        while (right > endp + UNROLL) {
            for (int i = UNROLL; i--;) {
                int h = comp(*right, piv);
                *rwr = *sw = *right--;
                rwr -= !h; sw += h;
            }
        }
    }

    while (right - left >= UNROLL &&
            (rwr - right > UNROLL || left - lwr > UNROLL)) {

        while (rwr - right > UNROLL && right - left >= UNROLL) {
            for (int i = UNROLL; i--;) {
                int h = comp(*left, piv);
                *lwr = *rwr = *left++;
                lwr += h; rwr -= !h;
            }
        }
        while (left - lwr > UNROLL && right - left >= UNROLL) {
            for (int i = UNROLL; i--;) {
                int h = comp(*right, piv);
                *rwr = *lwr = *right--;
                rwr -= !h; lwr += h;
            }
        }
    }
    while ((lwr < left||rwr > right) && left <= right) {
        while (rwr > right && left <= right) {
            T x = *left++;
            if (comp(x, piv)) *lwr++ = x;
            else *rwr-- = x;
        }
        while (lwr < left && left <= right) {
            T x = *right--;
            if (comp(x, piv)) *lwr++ = x;
            else *rwr-- = x;
        }
    }
    while (left <= right && !comp(*right, piv)) {
        right--;
        rwr--;
    }
    std::memcpy(lwr, swbuf, (sw - swbuf) * sizeof(T));
    *outerleft = *rwr;
    *rwr = piv;
    return rwr;
}

template<typename T, typename Compare>
void smallsort(T* left, T* right, Compare comp) {
    while (right - left > 15) {
        T* mid = partition_small(left, right, comp);
        smallsort(left, mid - 1, comp);
        left = mid + 1;
    }
    sorting_network(left, right - left, comp);
}

template<typename T, typename Compare>
void blqsort(T* left, T* right, Compare comp) {
    while (1) {
        ptrdiff_t partszm1 = right - left;
        if (partszm1 <= SMALLPART) break;

        T* mid = partition_large(left, right, comp);

        if (mid - left < partszm1 / 16) {
            if (mid > left) blqsort(left, mid - 1, comp);
            T piv = *mid;
            mid += 1;
            // collect duplicates
            for (T* p = mid; p <= right; p++) {
                if (!comp(piv, *p)) {
                    std::swap(*mid, *p);
                    mid++;
                }
            }
            left = mid;
            if (right - left < SMALLPART) break;
            // second chance before fallback to heapsort
            mid = partition_large(left, right, comp);
            if (mid - left < partszm1 / 16) {
                if (mid > left) heap_sort(left, mid - 1, comp);
                if (mid < right) heap_sort(mid + 1, right, comp);
                return;
            }
        }
        if (mid - left < right - mid) {
            blqsort(left, mid - 1, comp);
            left = mid + 1;
        }
        else {
            blqsort(mid + 1, right, comp);
            right = mid - 1;
        }
    }
    smallsort(left, right, comp);
}

template <typename T, typename Compare = std::less<T>>
void sort(T* first, T* last, Compare comp = Compare()) {

    if (last - first < 2) return;
    blqsort(first, last - 1, comp);
}
}
#endif

test.cpp

#include <cstdio>
#include <cstdlib>
#include <chrono>

#include "blqs.h"

constexpr int SIZE = 50000000;
double data[SIZE];

double time() {
    static auto time0 = std::chrono::high_resolution_clock::now();
    auto t = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double> diff = t - time0;
    return diff.count();
}

int main() {
    double t0;
    printf("Sorting %d million doubles\n", SIZE / 1000000);
    for (int i = 0; i < SIZE; i++) data[i] = rand() / 1024.0;
    t0 = time();
    blqs::sort(data, data + SIZE);
    printf("blqs::sort: %.2fs\n", time() - t0);
}

c++ -O3 test.cpp && ./a.out

Sorting 50 million doubles
blqs::sort: 0.95s

std::sort は1.33秒かかりました。x86システムでは、これはstd::sortの3倍の速さになります。

blqsortに含まれるその他の機能

blqs.h をインクルードするだけで、std::sort と同じくらい簡単に使用できます。

C++ マルチスレッド版（C++ スレッドを使用）では、blqs.h の代わりに blqs_thr.h をインクルードしてください。関数呼び出しは同じです。

実際には、カスタムデータ構造をソートする必要があることがよくあります。ここで、Google Highway のようなSIMDライブラリは、単純な数値には非常に高速ですが、使用が困難になります。

std::sort または blqsort を使用すると、より柔軟性が得られます。

#include "blqs.h"

struct entry {
    int id;
    int value;
    bool operator<(const entry& other) const { return id < other.id; }
};

struct entry data[SIZE];

blqs::sort(data, data + SIZE);

コンパイラとプラットフォームによっては、if バージョンまたは明示的な分岐なしバージョンが使用されます。

コピーコストが高く is_trivially_copyable でない型（文字列など）の場合、バッファベースの分岐なしアプローチは非効率的になります。このような場合、BlockQuicksort（Edelkamp および Weiß）のバリアントが使用されます。これは、要素インデックスのみを分岐なしで処理します。これにはC++17以降が必要です。

c++ -O3 -std=c++17 test_double.cpp && ./a.out

---

christof.kaser@gmail.com - Blog

---