ソート済み配列から値を検索する必要がある場合があります。最も単純なアルゴリズムは、値が見つかるか配列が尽きるまで、値を1つずつ順に調べる方法です。このアルゴリズムを線形探索と呼ぶこともあります。C++では、std::find関数でこの効果を得られます。

大規模な配列では、二分探索の方が優れています。二分探索は、ソート済み配列内の目的値を効率的に見つける古典的なアルゴリズムで、探索区間を繰り返し半分に分割します。配列全体から開始し、目的値と中央の要素を比較します。目的値が小さければ上半分を破棄し、大きければ下半分を破棄します。この処理は、目的値が見つかるか区間が空になるまで続けられます。大規模なデータセットでは線形探索よりもはるかに高速です。C++では、値が存在するかどうかを示す真偽値を返すstd::binary_search関数で実装されています。

人気のRoaring Bitmap形式は、サイズが1から4096の16ビット整数の配列を使用します。値が存在するかどうかを確認する必要がある場合があります。そこで二分探索を使用しています。

より高速なアプローチを求めました。そこで2つの洞察を得ました。

  1. 今日のほぼすべてのプロセッサには、データ並列命令（SIMDと呼ばれることもあります）があり、一度に複数の値をチェックできます。64ビットARMおよびx64プロセッサ（Intel/AMD）は、常に単一の命令で8つの16ビット整数と目的値を比較できます。これは、8要素より小さいブロックまで二分探索を進める必要がないことを示唆しています。また、16要素以上を安価に比較することも望ましいでしょう。
  2. 二分探索は一度に1つの値をチェックします。しかし、最近のプロセッサは一度に複数の値をロードしてチェックできます。優れたメモリレベルの並列性を持っています。これは、二分探索の代わりに四分探索を試みることを示唆しています。配列を半分に分割するのではなく、4分の1に分割するのです。その結果、命令数は少し増えるかもしれませんが、命令数がボトルネックになる可能性は低いでしょう。

そこで、私はSIMD Quadアルゴリズムと呼ぶものを作成しました。これは、16ビット符号なし整数のソート済み配列向けの効率的な探索アルゴリズムで、四分補間探索とSIMD（Single Instruction, Multiple Data）を組み合わせています。このアルゴリズムは、配列を16要素の固定サイズブロック（最後のブロックは除く場合あり）に分割し、各ブロックの最後の要素を補間キーとして使用して、探索対象を1つのブロックに素早く絞り込み、その後SIMD命令を使ってそのブロック内の16要素を同時にチェックします。

核となる考え方は、階層的な探索を行うことです。まず、粗いレベル（ブロック境界）で補間探索を使用して、目的値を含む可能性の高いブロックを見つけ、その後ブロック内で細かい並列チェックにSIMDに切り替えます。このハイブリッドアプローチは、アルゴリズム最適化（補間探索は比較回数を対数的に削減）とハードウェアアクセラレーション（SIMDは複数の要素を一度にチェック）の両方の強みを活用します。

  1. 初期チェック: 配列の要素数が16未満の場合は、すべての要素に対して単純な線形探索を実行します。
  2. ブロック分割: 配列を16個の連続した要素のブロックに分割します。サイズがcardinalityの配列の場合、num_blocks = cardinality / 16のフルブロックが存在します。
  3. 四分補間探索: 各ブロックの最後の要素（位置16-1、32-1など）を補間のキーとして使用します。探索は四分（基数4）補間を実行して、目的のposが位置する可能性の高いブロックを見つけます。これには、現在の探索範囲の4分の1ポイントと目的値を比較し、ベースを適宜調整することが含まれます。
  4. ブロック選択: 絞り込み後、補間結果に基づいて適切なブロックインデックスloを選択します。
  5. SIMDチェック: 有効なブロックが見つかった場合、16要素をSIMDレジスタにロードし（ARMではNEON、x64ではSSE2を使用）、目的値と並列に等価比較を実行します。いずれかが一致した場合、trueを返します。
  6. 残余チェック: フルブロックに含まれない要素（残余）については、線形探索を実行します。

性能はどうでしょうか。ベンチマークを作成しました。ベンチマークは次のように動作します。要素数2から4096までの各配列サイズについて、16ビット符号なし整数のソート済み配列を100,000個生成します。各サイズについて、「コールド」モード（各クエリが異なる配列を検索し、キャッシュミスをシミュレート）で1,000万回のメンバーシップクエリと、「ウォーム」モード（クエリが配列ごとにグループ化され、各配列が100回連続して検索され、キャッシュヒットをシミュレート）で1,000万回のクエリを実行します。ベンチマークは、線形探索（std::find）、二分探索（std::binary_search）、および新しいSIMD Quadアルゴリズムの3つのアルゴリズムについて、クエリあたりの平均時間を測定します。

2つのシステムを使用します。Apple LLVM搭載のApple M4と、GCC搭載のIntel Emerald Rapidsプロセッサです。

まず、線形探索と二分探索を比較してみましょう。

Intel/GCC:

Apple/LLVM

結果は明らかです。配列が大きくなると、二分探索は線形探索を上回ります。これは当然のことです。

コールドキャッシュでは、線形探索の方が相対的に悪化します。これは、アクセスするデータが多くなり、キャッシュフォルトが増えるため当然です。

二分探索が線形探索に対して全体的な勝者であることが確認できました。次にSIMD Quadアルゴリズムと比較してみましょう。

Intel/GCC:

Apple/LLVM

結果はIntelプラットフォームとAppleプラットフォームで大きく異なります。Intelプラットフォームでは、SIMD Quadはウォームキャッシュで二分探索の2倍以上高速です。コールドキャッシュでは利点が小さくなります。Appleプラットフォームでは逆で、コールドキャッシュでSIMD Quadが2倍以上高速になる一方、ウォームキャッシュでは利点がより限定的です。

しかし重要な点は、すべてのケースでSIMD Quadが二分探索よりも高速であることです。

アルゴリズムのSIMD部分はかなり単純です。作業を節約する専用命令を使用するため、高速化の理由は容易に理解できます。命令数が少なく、分岐も少なくなります。

しかし「quad」部分はどうでしょうか。重要でしょうか。そこで、同じアルゴリズムのバイナリ版を試してみました。SIMD最適化は同じですが、四分補間探索を削除し、標準的な二分探索に置き換えています。

Intel/GCC:

Apple/LLVM

簡単に言うと、quadアプローチはAppleプラットフォームではほとんど効果がありませんが、Intelプラットフォームでは大規模配列のコールドケースで適切な最適化となります。四分探索は、Intelサーバーのメモリレベルの並列性をより良く活用します。

ソースコードはこちらで入手できます。

結論。私の結果が示唆するのは、教科書的な二分探索はまともなアルゴリズムですが、意味のある方法で改善できるということです。標準アルゴリズムは、非常に多くの並列性を持つコンピュータ向けに設計されていないことが多かったのです。SIMD Quadアルゴリズムは、メモリレベルとデータ並列性の両方を活用しようとしています。さらに、私のアルゴリズムよりもさらに改善できると私は考えています。創造性を発揮しましょう！

参考文献: Faster intersections between sorted arrays with shotgun

付録（ソースコード）

bool simd_quad(const uint16_t *carr, int32_t cardinality, 
            uint16_t pos) {
    constexpr int32_t gap = 16;
    if (cardinality < gap) {
      for (int32_t j = 0; j < cardinality; j++) {
          if (carr[j] == pos) return true;
        }
        return false;
    }
    int32_t num_blocks = cardinality / gap;
    int32_t base = 0;
    int32_t n = num_blocks;
    while (n > 3) {
      int32_t quarter = n >> 2;

      int32_t k1 = carr[(base + quarter + 1) * gap - 1];
      int32_t k2 = carr[(base + 2 * quarter + 1) * gap - 1];
      int32_t k3 = carr[(base + 3 * quarter + 1) * gap - 1];

      int32_t c1 = (k1 < pos);
      int32_t c2 = (k2 < pos);
      int32_t c3 = (k3 < pos);

      base += (c1 + c2 + c3) * quarter;
      n -= 3 * quarter;
    }
    while (n > 1) {
        int32_t half = n >> 1;
        base = (carr[(base + half + 1) * gap - 1] < pos) 
                 ? base + half : base;
        n -= half;
    }
    int32_t lo = (carr[(base + 1) * gap - 1] < pos) 
                ? base + 1 : base;

    if (lo < num_blocks) {
        const uint16_t *blk = carr + lo * gap;
#ifdef __ARM_NEON
        uint16x8_t needle = vdupq_n_u16(pos);
        uint16x8_t v0 = vld1q_u16(blk);
        uint16x8_t v1 = vld1q_u16(blk + 8);
        uint16x8_t hit = vorrq_u16(vceqq_u16(v0, needle), 
                  vceqq_u16(v1, needle));
        return vmaxvq_u16(hit) != 0;
#else
        __m128i needle = _mm_set1_epi16((short)pos);
        __m128i v0 = _mm_loadu_si128((const __m128i *)blk);
        __m128i v1 = _mm_loadu_si128((const __m128i *)(blk + 8));
        __m128i hit = _mm_or_si128(_mm_cmpeq_epi16(v0, needle),
                                   _mm_cmpeq_epi16(v1, needle));
        return _mm_movemask_epi8(hit) != 0;
#endif
    }

    for (int32_t j = num_blocks * gap; j < cardinality; j++) {
        uint16_t v = carr[j];
        if (v >= pos) return (v == pos);
    }
    return false;
}