この記事はビット演算テクニック Advent Calendar 2016

の11日目の記事です。

ビット単位のunpack

C++の固定長ビット配列を扱うクラスstd::bitsetを例として解説します。

std::bitset<128> unpacklo(const std::bitset<128> &a, const std::bitset<128> &b) {
  std::bitset<128> result;
  for (int i = 0; i < 64; ++i) {
    result[2*i]   = a[i];
    result[2*i+1] = b[i];
  }
  return result;
}

std::bitset<128> unpackhi(const std::bitset<128> &a, const std::bitset<128> &b) {
  std::bitset<128> result;
  for (int i = 0; i < 64; ++i) {
    result[2*i]   = a[i+64];
    result[2*i+1] = b[i+64];
  }
  return result;
}

要は2つのビット列から交互に1ビットずつ取ってきたデータを生成します。今回は簡単のため128ビット固定とします。

SSE

普通にstd::bitsetを使っていては遅いのでintelのx86の拡張命令セットであるSSE(Streaming SIMD Extensions)を利用します。 (追加された時期によりSSE, SSE2, ..., SSE4.2がありますがここではまとめてSSEと呼ぶことにします)

SSEを用いると128ビット単位で算術演算、ビット演算、シャッフルなど様々な操作を行うことができます。 SSEには8ビット単位のunpackを行うpunpacklbw/punpackhbwという命令があるので、まずこれを用いて8ビットごとにunpackし、残りを3日目に解説したDelta Swapを使ってunpackするのがぱっと思いつく実装でしょうか。 intrinsicsを使って実装してみます。

#include <x86intrin.h>

__m128i mm_unpacklo_epb_unpack_dswap(__m128i a, __m128i b) {
  __m128i unpack8 = _mm_unpacklo_epi8(a, b);
  __m128i unpack4 = mm_delta_swap_epi64(unpack8, _mm_set1_epi16(0x00F0), 4);
  __m128i unpack2 = mm_delta_swap_epi64(unpack4, _mm_set1_epi8(0x0C), 2);
  return mm_delta_swap_epi64(unpack2, _mm_set1_epi8(0x22), 1);
}

__m128i mm_unpackhi_epb_unpack_dswap(__m128i a, __m128i b) {
  __m128i unpack8 = _mm_unpackhi_epi8(a, b);
  __m128i unpack4 = mm_delta_swap_epi64(unpack8, _mm_set1_epi16(0x00F0), 4);
  __m128i unpack2 = mm_delta_swap_epi64(unpack4, _mm_set1_epi8(0x0C), 2);
  return mm_delta_swap_epi64(unpack2, _mm_set1_epi8(0x22), 1);
}

pdep

ビット演算に詳しい人ならBMIという拡張命令セットにpdepという命令があってこれが使えるということを思いつくかもしれません。

#include <x86intrin.h>

__m128i mm_unpacklo_epb_pdep(__m128i a, __m128i b) {
  uint64_t alo = _mm_cvtsi128_si64(a);
  uint64_t blo = _mm_cvtsi128_si64(b);
  return _mm_set_epi64x(
      _pdep_u64(alo >> 32, UINT64_C(0x5555555555555555)) | _pdep_u64(blo >> 32, UINT64_C(0xAAAAAAAAAAAAAAAA)),
      _pdep_u64(alo, UINT64_C(0x5555555555555555)) | _pdep_u64(blo, UINT64_C(0xAAAAAAAAAAAAAAAA)));
}

__m128i mm_unpackhi_epb_pdep(__m128i a, __m128i b) {
  uint64_t ahi = _mm_extract_epi64(a, 1);
  uint64_t bhi = _mm_extract_epi64(b, 1);
  return _mm_set_epi64x(
      _pdep_u64(ahi >> 32, UINT64_C(0x5555555555555555)) | _pdep_u64(bhi >> 32, UINT64_C(0xAAAAAAAAAAAAAAAA)),
      _pdep_u64(ahi, UINT64_C(0x5555555555555555)) | _pdep_u64(bhi, UINT64_C(0xAAAAAAAAAAAAAAAA)));
}

pclmulqdq

さて、また別の拡張命令セットCLMULにはpclmulqdqという命令があります。これは繰り上がりを無視しての2進数の掛け算を計算する命令で、CRCの計算やガロア体 ${GF(2^k)}$ 上の乗算を高速に行うために使われます。 pclmulqdqの動作をstd::bitsetを用いた疑似コードで書くとこんな感じでしょうか。

std::bitset<128> pseudo_pclmulqdq(const std::bitset<128> &a, const std::bitset<128> &b, const int index) {
  std::bitset<64> x, y;
  if (index & 0x01) {
    for (int i = 0; i < 64; ++i) {
      x[i] = a[i+64];
    }
  } else {
    for (int i = 0; i < 64; ++i) {
      x[i] = a[i];
    }
  }
  if (index & 0x10) {
    for (int i = 0; i < 64; ++i) {
      y[i] = b[i+64];
    }
  } else {
    for (int i = 0; i < 64; ++i) {
      y[i] = b[i];
    }
  }
  std::bitset<128> result;
  for (int i = 0; i < 64; ++i) {
    for (int j = 0; j < 64; ++j) {
      if (x[i] && y[j]) {
        result.flip(i+j);
      }
    }
  }
  return result;
}

最初のindexにより分岐している部分で128ビットのうちどの64ビットを使って演算するかを選択し、64ビットx64ビットの繰り上がりなしの掛け算を行います。繰り上がりなしの掛け算は有限体 $F_2$ 係数の多項式 $F_2[x]$ の乗算だと思うことができます。(足し算はXOR、掛け算はANDに相当するので)

$pclmulqdq(\sum_{i=0}^na_ix^i, \sum_{j=0}^mb_jx^j) = \sum_{i=0}^n\sum_{j=0}^m a_ib_jx^{i+j} = (\sum_{i=0}^na_ix^i)(\sum_{j=0}^mb_jx^j)$

ここで、 $F_2$ の多項式の自乗を考えてみます。

$(\sum_{i=0}^na_ix^i)^2 = \sum_{i=0}^na_ix^{2i}+2\sum_{0 \leq i < j \leq n}a_ia_jx^{i+j} = \sum_{i=0}^na_ix^{2i}$

これはつまり

for (int i = 0; i < 64; ++i) {
  result[2*i]   = a[i];
}

と同じことになります。これを用いるとかなり簡潔にunpackを実装できます。

__m128i mm_unpacklo_epb_pclmulqdq(__m128i a, __m128i b) {
  return _mm_or_si128(_mm_clmulepi64_si128(a, a, 0x00), _mm_slli_epi32(_mm_clmulepi64_si128(b, b, 0x00), 1));
}

__m128i mm_unpackhi_epb_pclmulqdq(__m128i a, __m128i b) {
  return _mm_or_si128(_mm_clmulepi64_si128(a, a, 0x11), _mm_slli_epi32(_mm_clmulepi64_si128(b, b, 0x11), 1));
}

ベンチマーク

以上3種類のアルゴリズムの実行コストを実際に計測してみました。

unpack.cpp

Hi/Loそれぞれ2³⁰回(≒10億回)unpackを計算しその時間を計測します。

1回の処理が軽いので、関数呼び出しやデータ生成が遅いと差が出る比較ができないので、（もともと最適化したら関数呼び出しは消えますが）関数にインライン化指定をして、データ生成はxorshiftを用いています。

ベンチマーク結果は以下のとおりです。

CPU: Core i7-6700K 4.0GHz Turbo Boost 4.2GHz
Memory: DDR-2133 64GB
OS: Arch Linux (Linux 4.8.10-1-ARCH x86_64)
Compiler: GCC 6.2.1
Compiler Options: -O3 -flto -mtune=skylake -march=skylake -lboost_system -lboost_timer

アルゴリズム	unpacklo	unpackhi
8bitunpack+Delta Swap	3.767s	3.792s
pdep	2.624s	2.938s
pclmulqdq	2.252s	2.265s

Hi/Loともにpclmulqdqを用いた場合が一番高速という結果になりました。ちなみにpdepバージョンがHiとLoで実行時間に開きがあるのはmovq命令とpextrq命令の違いによると思われます。

pclmulqdq命令は使いどころが難しいですがスポッとハマると気持ちいいですね。

さいごに

ビット演算テクニックアドベントカレンダー明日の12日目はocxtalさんによる「たぶんbitDPの話」です。お楽しみに！

prime's diary

そすうの日々を垂れ流しちゃうやつだよ

pclmulqdqを用いたビット単位のunpack【ビット演算テクニック Advent Calendar 2016 11日目】

ビット単位のunpack

SSE

pdep

pclmulqdq

ベンチマーク

さいごに