2016-12-13

Bit scan reverseをAVXで並列化する

Bit scan reverse (bsr)

00000001010100010100000100101101
<-clz->*<---bit scan reverse--->

1になっている最も高位のビットのインデックスを0-indexedで計算します。0に対する結果は未定義（不定）です。似た処理にcount leading zero(clz)があり、こちらは最上位ビット側から数えます。(0に対する結果は整数型のビット数num_bitsになる) したがって、0以外の値に対してはbsr + 1 + clz = num_bitsが成り立ちます。

bsr/lzcnt命令

x86にはズバリbsr命令があり、32ビット/64ビットの値に対してBit scan reverseを計算できます。また、BMI拡張命令セットにはlzcntという命令があり、clzも1命令で計算できます。

しかし、bsr、lzcntともにSSEやAVX拡張命令セットの中に対応する命令がないため、並列に計算することができません。

今回はAVX命令を使って256ビット幅で並列にBit scan reverseを計算する方法を考えます。

また、8/16/32/64ビットそれぞれについて、どのアルゴリズムが最適なのかベンチマークを取ります。

アルゴリズム

スカラーに変換してbsrで処理 (naive)

まずは単純にスカラー値に変換してbsr命令で処理するのが一番自然でしょう。 pextr{b,w,d,q}は意外にコストが高いのでvmovdqaで一括でメモリ上（実際にはL1キャッシュに乗るはず）にコピーして処理することにします。計算結果をYMMレジスタに格納するのは_mm256_set_epi* intrinsicsに任せました（メモリ上に書き戻してvmovdqaで格納するより速かった）

inline __m256i bsr_256_8_naive(__m256i x) {
  alignas(32) uint8_t b[32];
  _mm256_store_si256((__m256i*)b, x);
  return _mm256_setr_epi8(
      __bsrd(b[ 0]), __bsrd(b[ 1]), __bsrd(b[ 2]), __bsrd(b[ 3]),
      __bsrd(b[ 4]), __bsrd(b[ 5]), __bsrd(b[ 6]), __bsrd(b[ 7]),
      __bsrd(b[ 8]), __bsrd(b[ 9]), __bsrd(b[10]), __bsrd(b[11]),
      __bsrd(b[12]), __bsrd(b[13]), __bsrd(b[14]), __bsrd(b[15]),
      __bsrd(b[16]), __bsrd(b[17]), __bsrd(b[18]), __bsrd(b[19]),
      __bsrd(b[20]), __bsrd(b[21]), __bsrd(b[22]), __bsrd(b[23]),
      __bsrd(b[24]), __bsrd(b[25]), __bsrd(b[26]), __bsrd(b[27]),
      __bsrd(b[28]), __bsrd(b[29]), __bsrd(b[30]), __bsrd(b[31]));
}

inline __m256i bsr_256_16_naive(__m256i x) {
  alignas(32) uint16_t b[16];
  _mm256_store_si256((__m256i*)b, x);
  return _mm256_setr_epi16(
      __bsrd(b[ 0]), __bsrd(b[ 1]), __bsrd(b[ 2]), __bsrd(b[ 3]),
      __bsrd(b[ 4]), __bsrd(b[ 5]), __bsrd(b[ 6]), __bsrd(b[ 7]),
      __bsrd(b[ 8]), __bsrd(b[ 9]), __bsrd(b[10]), __bsrd(b[11]),
      __bsrd(b[12]), __bsrd(b[13]), __bsrd(b[14]), __bsrd(b[15]));
}

inline __m256i bsr_256_32_naive(__m256i x) {
  alignas(32) uint32_t b[8];
  _mm256_store_si256((__m256i*)b, x);
  return _mm256_setr_epi32(
      __bsrd(b[0]), __bsrd(b[1]), __bsrd(b[2]), __bsrd(b[3]),
      __bsrd(b[4]), __bsrd(b[5]), __bsrd(b[6]), __bsrd(b[7]));
}

inline __m256i bsr_256_64_naive(__m256i x) {
  alignas(32) uint64_t b[4];
  _mm256_store_si256((__m256i*)b, x);
  return _mm256_setr_epi64x(
      __bsrq(b[0]), __bsrq(b[1]), __bsrq(b[2]), __bsrq(b[3]));
}

popcountに帰着 (popcount)

uint32_t bsr32(uint32_t x) {
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  return popcount32(x>>1);
}

右シフトとORを組み合わせることで1の立っている一番高位のビットより下位のビットを1で埋めることができます。

これをpopcountと組み合わせればbsrが求められます。

inline __m256i popcount_256_8(__m256i x) {
  x = _mm256_sub_epi8(x, _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi8(0x55)));
  x = _mm256_add_epi8(_mm256_and_si256(x, _mm256_set1_epi8(0x33)), _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi8(0x33)));
  return _mm256_and_si256(_mm256_add_epi8(x, _mm256_srli_epi16(x, 4)), _mm256_set1_epi8(0x0F));
}

inline __m256i bsr_256_8_popcnt(__m256i x) {
  x = _mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi8(0x7F)));
  x = _mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi8(0x3F)));
  x = _mm256_and_si256(_mm256_srli_epi16(_mm256_or_si256(x, _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi8(0x0F))), 1), _mm256_set1_epi8(0x7F));
  return popcount_256_8(x);
}

inline __m256i popcount_256_16(__m256i x) {
  x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x5555)), _mm256_and_si256(_mm256_srli_epi16(x, 1), _mm256_set1_epi16(0x5555)));
  x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x3333)), _mm256_and_si256(_mm256_srli_epi16(x, 2), _mm256_set1_epi16(0x3333)));
  x = _mm256_add_epi16(_mm256_and_si256(x, _mm256_set1_epi16(0x0F0F)), _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi16(0x0F0F)));
  return _mm256_maddubs_epi16(x, _mm256_set1_epi16(0x0101));
}

inline __m256i bsr_256_16_popcnt(__m256i x) {
  x = _mm256_or_si256(x, _mm256_srli_epi16(x, 1));
  x = _mm256_or_si256(x, _mm256_srli_epi16(x, 2));
  x = _mm256_or_si256(x, _mm256_srli_epi16(x, 4));
  x = _mm256_srli_epi16(_mm256_or_si256(x, _mm256_srli_epi16(x, 8)), 1);
  return popcount_256_16(x);
}

inline __m256i popcount_256_32(__m256i x) {
  x = _mm256_sub_epi32(x, _mm256_and_si256(_mm256_srli_epi32(x, 1), _mm256_set1_epi8(0x55)));
  x = _mm256_add_epi32(_mm256_and_si256(x, _mm256_set1_epi8(0x33)), _mm256_and_si256(_mm256_srli_epi32(x, 2), _mm256_set1_epi8(0x33)));
  x = _mm256_and_si256(_mm256_add_epi8(x, _mm256_srli_epi32(x, 4)), _mm256_set1_epi8(0x0F));
  return _mm256_srli_epi16(_mm256_madd_epi16(x, _mm256_set1_epi16(0x0101)), 8);
}

inline __m256i bsr_256_32_popcnt(__m256i x) {
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 1));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 2));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 4));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 8));
  x = _mm256_srli_epi32(_mm256_or_si256(x, _mm256_srli_epi32(x, 16)), 1);
  return popcount_256_32(x);
}

inline __m256i bsr_256_64_popcnt(__m256i x) {
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 1));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 2));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 4));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 8));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 16));
  x = _mm256_srli_epi64(_mm256_or_si256(x, _mm256_srli_epi64(x, 32)), 1);
  return popcount_256_64(x);
}

シャッフル命令vpshufbを使った最適化 (rev popcount)

64ビット単位のときのみ、しかも若干効果がある程度ですが、vpshufbを使ってバイト順を反転させ、立っている一番下位のビットを算術演算で求めるテクニックを使うことができます。

inline __m256i bswap_256_64(__m256i x) {
  __m128i swap_table_128 = _mm_setr_epi8(7, 6, 5, 4, 3, 2, 1, 0, 15, 14, 13, 12, 11, 10, 9, 8);
  __m256i swap_table = _mm256_broadcastsi128_si256(swap_table_128);
  return _mm256_shuffle_epi8(x, swap_table);
}

inline __m256i bsr_256_64_rev_popcnt(__m256i x) {
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 1));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 2));
  x = _mm256_or_si256(x, _mm256_srli_epi64(x, 4));
  x = _mm256_andnot_si256(_mm256_srli_epi64(x, 1), x);
  x = bswap_256_64(x);
  x = _mm256_and_si256(x, _mm256_sub_epi64(_mm256_setzero_si256(), x));
  x = bswap_256_64(x);
  return popcount_256_64(_mm256_sub_epi64(x, _mm256_set1_epi64x(1)));
}

ギャザー命令を使ったテーブル参照 (table gather)

32/64ビット単位で表引きができるvpgather{d,q}{d,q}を使ってpopcountを表引きします。32ビット単位なら4byte*2^31byte = 8GBに収まります。表引きは000111...111の形のものしかないので実は32通りしかないので、L1キャッシュに乗る可能性が高いです。しかもmallocでメモリを確保してもアクセスしていない領域には実際にはページ割当がされないので物理メモリ使用量はそこまで大きくなりません。

int32_t *table_32;

void init_table_32() {
  table_32 = (int32_t*)malloc(sizeof(int32_t)*((size_t)1 << 31));
  table_32[0] = 0;
  int64_t j = 1;
  for (int i = 1; i < 31; ++i) {
    table_32[j] = i;
    j = 2*j + 1;
  }
}

inline __m256i bsr_256_32_table_gather(__m256i x) {
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 1));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 2));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 4));
  x = _mm256_or_si256(x, _mm256_srli_epi32(x, 8));
  x = _mm256_srli_epi32(_mm256_or_si256(x, _mm256_srli_epi32(x, 16)), 1);
  return _mm256_i32gather_epi32(table_32, x, 4);
}

vcvtdq2ps命令による浮動小数点数への変換 (cvtfloat)

浮動小数点数のビット表現は

kivantium.hateblo.jp

などを参照してください。

vcvtdq2ps命令で浮動小数点数へ変換して、指数部だけになるようにシフトして、下駄履き(単精度なら127)を引いてやることでbsrが得られます。しかしこれには罠があります。具体的には、

ビットが浮動小数点数の精度以上に連続していると丸めによって間違った結果になる(1大きくなってしまう)。
最上位ビットは符号として扱われるので間違った結果になる(負数扱いになり全く違う値になる)。

これらへの対策については

或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化

或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化（その２）

を参照してください。

16ビット以下では対策は不要になります。また、64ビット版では64ビットを3分割することで問題を回避したバージョンも示します。

inline __m256i bsr_256_32_cvtfloat_impl(__m256i x, int32_t sub) {
  __m256i cvt_fl = _mm256_castps_si256(_mm256_cvtepi32_ps(x));
  __m256i shifted = _mm256_srli_epi32(cvt_fl, 23);
  return _mm256_sub_epi32(shifted, _mm256_set1_epi32(sub));
}

inline __m256i bsr_256_8_cvtfloat(__m256i x) {
  __m256i r0 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi32(0x000000FF)), 127);
  __m256i r1 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi32(x,  8), _mm256_set1_epi32(0x000000FF)), 127);
  __m256i r2 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi32(x, 16), _mm256_set1_epi32(0x000000FF)), 127);
  __m256i r3 = bsr_256_32_cvtfloat_impl(_mm256_srli_epi32(x, 24), 127);
  __m256i r02 = _mm256_blend_epi16(r0, _mm256_slli_epi32(r2, 16), 0xAA);
  __m256i r13 = _mm256_blend_epi16(r1, _mm256_slli_epi32(r3, 16), 0xAA);
  return _mm256_blendv_epi8(r02, _mm256_slli_epi16(r13, 8), _mm256_set1_epi16(0xFF00));
}

inline __m256i bsr_256_16_cvtfloat(__m256i x) {
  __m256i lo = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi32(0x0000FFFF)), 127);
  __m256i hi = bsr_256_32_cvtfloat_impl(_mm256_srli_epi32(x, 16), 127);
  return _mm256_blend_epi16(lo, _mm256_slli_epi32(hi, 16), 0xAA);
}

inline __m256i bsr_256_32_cvtfloat(__m256i x) {
  x = _mm256_andnot_si256(_mm256_srli_epi32(x, 1), x); // 連続するビット対策
  __m256i result = bsr_256_32_cvtfloat_impl(x, 127);
  result = _mm256_or_si256(result, _mm256_srai_epi32(x, 31));
  return _mm256_and_si256(result, _mm256_set1_epi32(0x0000001F));
}

inline __m256i bsr_256_64_cvtfloat(__m256i x) {
  __m256i bsr32 = bsr_256_32_cvtfloat(x);
  __m256i higher = _mm256_add_epi32(_mm256_srli_epi64(bsr32, 32), _mm256_set1_epi64x(0x0000000000000020));
  __m256i mask = _mm256_shuffle_epi32(_mm256_cmpeq_epi32(x, _mm256_setzero_si256()), 0xF5);
  return _mm256_blendv_epi8(higher, _mm256_and_si256(bsr32, _mm256_set1_epi64x(0xFFFFFFFF)), mask);
}

inline __m256i bsr_256_64_cvtfloat2(__m256i x) {
  __m256i bsr0 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(x, _mm256_set1_epi64x(0x3FFFFF)), 127);
  __m256i bsr1 = bsr_256_32_cvtfloat_impl(_mm256_and_si256(_mm256_srli_epi64(x, 22), _mm256_set1_epi64x(0x3FFFFF)), 105);
  __m256i bsr2 = bsr_256_32_cvtfloat_impl(_mm256_srli_epi64(x, 44), 83);
  __m256i result = _mm256_max_epi32(_mm256_max_epi32(bsr0, bsr1), bsr2);
  return _mm256_and_si256(result, _mm256_set1_epi64x(0xFFFFFFFF));
}

vpshufbを使ったテーブル参照 (table)

バイト単位のシャッフル命令vpshufbは、バイト列の並べ替え以外に、4ビット入力、8ビット出力のルックアップテーブルとして用いることができます。これを用いて、4ビット単位でbsrを求めて、それを組み合わせることで8ビット以上のbsrを求めることができます。

inline __m256i bsr_256_8_table(__m256i x) {
  __m128i table128_lo = _mm_setr_epi8(0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3);
  __m128i table128_hi = _mm_setr_epi8(0, 4, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7);
  __m256i table_lo = _mm256_broadcastsi128_si256(table128_lo);
  __m256i table_hi = _mm256_broadcastsi128_si256(table128_hi);
  return _mm256_max_epi8(
      _mm256_shuffle_epi8(table_lo, _mm256_and_si256(x, _mm256_set1_epi8(0x0F))),
      _mm256_shuffle_epi8(table_hi, _mm256_and_si256(_mm256_srli_epi16(x, 4), _mm256_set1_epi8(0x0F))));
}

inline __m256i bsr_256_16_table(__m256i x) {
  __m256i half = bsr_256_8_table(x);
  __m256i mask = _mm256_cmpgt_epi8(x, _mm256_setzero_si256());
  __m256i res8 = _mm256_add_epi8(half, _mm256_and_si256(mask, _mm256_set1_epi16(0x0800)));
  return _mm256_max_epi16(_mm256_and_si256(res8, _mm256_set1_epi16(0x00FF)), _mm256_srli_epi16(res8, 8)); 
}

ベンチマーク

Haswell/Broadwell/Skylakeの各アーキテクチャのCPUで各アルゴリズムを並列に実行した時のの実行時間を計測しました。

bsrが合計2³¹回実行されるように、8ビット幅なら2²⁶回、64ビット幅なら2²⁹回ループを回します。

Parallel bit scan reverse

環境

Haswell:

Google compute engine n1-highmem-2 (2vCPU, 13GB mem)
Xeon 2.3GHz, Turbo Boost時は不明
OS: Ubuntu 16.10
Compiler: GCC 6.1.1

Broadwell:

Google compute engine n1-highmem-2 (2vCPU, 13GB mem)
Xeon 2.2GHz, Turbo Boost時は不明
OS: Ubuntu 16.10
Compiler: GCC 6.1.1

Skylake:

ローカル
Core i7-6700K 4コア8スレッド 4GHz, Turbo Boost 4.2GHz 64GB mem
OS: Arch Linux
Compiler: GCC 6.2.1

コンパイルオプションはいずれも

-O3 -flto -mtune=(アーキテクチャ名) -march=(アーキテクチャ名) -lboost_system -lboost_timer

結果

(異なるアーキテクチャ同士で実行時間を比較するのはTurbo boostなどの影響があるのであまり意味がないと思います。)

8ビット

アルゴリズム	Haswell	Broadwell	Skylake
naive	1.323s	0.993s	0.832s
cvtfloat	0.284s	0.209s	0.165s
popcnt	0.264s	0.201s	0.169s
table	0.079s	0.058s	0.049s

いずれも表引きが一番高速なようです。

16ビット

アルゴリズム	Haswell	Broadwell	Skylake
naive	1.562s	1.192s	0.948s
cvtfloat	0.257s	0.200s	0.146s
popcnt	0.550s	0.423s	0.337s
table	0.262s	0.199s	0.179s

表引きあるいは浮動小数点数への変換が一番高速な結果になりました。

32ビット

アルゴリズム	Haswell	Broadwell	Skylake
naive	1.581s	1.220s	0.993s
cvtfloat	0.416s	0.318s	0.237s
popcnt	1.278s	0.966s	0.701s
table gather	1.202s	0.915s	0.523s

浮動小数点数への変換が圧倒的に速い結果になりました。 gather命令を使った表引きはSkylakeで若干高速になったようです。

64ビット

アルゴリズム	Haswell	Broadwell	Skylake
naive	1.574s	1.220s	0.911s
cvtfloat	1.369s	1.037s	0.804s
cvtfloat2	1.595s	1.225s	0.897s
popcnt	2.521s	1.924s	1.498s
rev popcnt	2.448s	1.898s	1.476s

ここでも浮動小数点数への変換が一番速い結果となりました。3分割してコーナーケースの処理を不要にしたバージョンは思ったほど性能が出ませんでした。 64ビット幅になるとpopcntは命令数が増えてナイーブな実装より遅くなってしまいました。

来るAVX512時代には

少なくとも32ビット幅以上なら並列でclzを求めるvplzcnt{d,q}命令が追加されるのでそれを使えば良さそうです。 8ビット幅の場合vpermi2bを使って7ビットの表引きをして求めるのがいいかもしれません。 16ビット幅の場合はやってみないとわからなさそうです。

参考文献

或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化

或るプログラマの一生 » x86 のビットスキャン命令の SIMD 化（その２）

bsrを並列処理する方法について考察した記事。当時はAVX2が使えるプロセッサはまだ発売されていなかった。

2016-12-11

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

ビット演算

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

www.adventar.org

の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の話」です。お楽しみに！

www.adventar.org

2016-12-03

Delta swapとその応用【ビット演算テクニック Advent Calendar 2016 3日目】

ビット演算

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

www.adventar.org

の三日目の記事です。

ビット列の部分列を入れ替える操作を効率的に行えるDelta swapを紹介します。xorが好きなのでDelta swapも好きなテクニックの一つです。

行える操作

マスクとずらす量を自由に指定できます。

abcdefghijklmnop 元のビット列

0000011000011100 マスク
00000fg0000lmn00 マスクを取ったビット列
00fg0000lmn00000 3ビット左にずらしたビット列

0011000011100000 3ビット左にずらしたマスク
00cd0000ijk00000 マスクを取ったビット列
00000cd0000ijk00 3ビット右にずらしたビット列

ab00e00h000000op 入れ替えない残りのビット列

abfgecdhlmnijkop 結果のビット列

実装

まずは素朴な実装をしてみます。

uint32_t swap_bits_32(uint32_t bits, uint32_t mask, int delta) {
  return ((bits & mask) << delta) | ((bits >> delta) & mask) | (bits & ~(mask | (mask << delta)));
}

それぞれの部分列、そしてそれ以外の場所をそれぞれマスクを取り、シフトしてORで一つにまとめます。(最適化してシフトしてからマスク取ってたりしますがまあいいでしょう)

それがDelta swapでは次のように書けます。

uint32_t delta_swap_32(uint32_t bits, uint32_t mask, int delta) {
  uint32_t x = (bits ^ (bits >> delta)) & mask;
  return bits ^ x ^ (x << delta);
}

だいぶシンプルになりましたね。ビット演算の演算回数も10回（andnotがあれば9回）から6回に減りました。

応用

Delta swapを複数回行ってビット列を任意に並べ替えることができます。

もちろん1ビットずつ入れ替えれば任意に並べ替えることができますが、実際にはより効率的に並べ替えることができます。 32ビットなら次のようにdeltaは固定でマスクをうまく設定することで9回のDelta swapで任意の並べ替えを実現できます。(2ⁿビットなら2n-1回(n>0))

bits = delta_swap_32(bits, mask1,  1);
bits = delta_swap_32(bits, mask2,  2);
bits = delta_swap_32(bits, mask3,  4);
bits = delta_swap_32(bits, mask4,  8);
bits = delta_swap_32(bits, mask5, 16);
bits = delta_swap_32(bits, mask6,  8);
bits = delta_swap_32(bits, mask7,  4);
bits = delta_swap_32(bits, mask8,  2);
bits = delta_swap_32(bits, mask9,  1);

証明は省略しますが、偶数ビット目と奇数ビット目に分けて再帰的に入れ替えを行って最後に偶奇を合わせる感じです。

bit matrix

16bit, 64bitのビット列をそれぞれ4x4、8x8のビットの行列(bit matrix)と思うと二次元のデータを効率的に扱えることがあります。

さて、Delta swapを使うと効率的にbit matrixを反転・回転させることができます。たとえば対角線を軸とした反転は次のように書けます。

uint64_t flipDiagA8H1(uint64_t bits) {
  bits = delta_swap(bits, 0x00000000F0F0F0F0, 28);
  bits = delta_swap(bits, 0x0000CCCC0000CCCC, 14);
  return delta_swap(bits, 0x00AA00AA00AA00AA,  7);
}

他の例はhttps://chessprogramming.wikispaces.com/Flipping+Mirroring+and+Rotatingにあります。

参考文献

https://chessprogramming.wikispaces.com/General+Setwise+Operations#SwappingBits
"The Art of Computer Programming volume 4" chess programming wiki曰く、Delta swapによる任意の入れ替えの証明が載っているそうです。

2016-12-01

ビット列による部分集合表現【ビット演算テクニック Advent Calendar 2016 1日目】

ビット演算競技プログラミング

はじめに

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

www.adventar.org

の1日目の記事です。

この記事ではビット演算を使って有限集合の部分集合を表現、操作するテクニックを紹介します。本文中で説明を省略した部分については、詳しい説明のある参考文献を紹介しておきます。

2021/01/03: 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

表現方法

S={A,B,C,D}という4要素の集合があったとします。各要素を次のように2進数の数字と対応付けます。 A: 0001 B: 0010 C: 0100 D: 1000 こうすると、Sの部分集合は存在する要素のビットごとの論理和を取ることで4ビットのビット列として一意に表すことができ、例えば {A,B,C}: 0111 {B,D}: 1010 のようになります。

以後、部分集合とそのビット列による表現をしばしば同一視します。

基本操作

空集合、全体集合

それぞれ0000, 1111が対応します。

和集合

和集合はビットごとの論理和(|)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∪U = {A,B,C,D}: 1111 = 0111 | 1010

共通部分

共通部分はビットごとの論理積(&)で書くことができます。

T = {A,B,C}: 0111, U = {B,D}: 1010

T∩U = {B}: 0010 = 0111 & 1010

補集合

補集合はビット反転（と必要に応じて&でマスクを取る操作）で書くことができます

T = {A,B,C}: 0111

T^c = {D}: 1000 = ~0111

(上位のビットが無視できない場合は 1000 = ~0111 & 1111)

部分集合全体を渡る

numを全体集合の要素数とすると、

for (int i = 0; i < (1 << num); ++i) {
  //ここに処理を書く
}

(本当はシフト演算子のところの括弧は不要ですが、見やすさの為付けています)

包含関係を判定する

補集合と共通部分を持つかどうかで判定できるので、 TがUに含まれるかどうかは~U&T==0で書けます。

要素数を数える

1になっているビットの数を数える操作(population count, popcount)が相当します。最近のIntelのCPUならずばりPOPCNT命令がありますし、そうでない場合も高速に処理する方法が知られています。 32要素の集合Sの部分集合Tの要素数は、

T = (T & 0x55555555) + ((T >>  1) & 0x55555555);
T = (T & 0x33333333) + ((T >>  2) & 0x33333333);
T = (T & 0x0F0F0F0F) + ((T >>  4) & 0x0F0F0F0F);
T = (T & 0x00FF00FF) + ((T >>  8) & 0x00FF00FF);
T = (T & 0x0000FFFF) + ((T >> 16) & 0x0000FFFF);

で計算できます。ここで、 0x5=0b0101 0x3=0b0011 となっていて、 n行目で2ⁿビットごとのpopcountを計算しています。さらに最適化して、

T -= (T >> 1) & 0x55555555;
T = (T & 0x33333333) + ((T >> 2) & 0x33333333);
T = (T + (T >> 4)) & 0x0F0F0F0F;
T = (T * 0x01010101) >> 24;

と書くこともできます。

応用操作

ある部分集合の部分集合全体を渡る

T⊃UなるUを全部列挙します。全体集合の部分集合全体を渡って、Tに含まれているかを判定すればできますが、要素数が少ない場合は無駄が多いですね。実は、次のようにすることでTの要素数をNとしたとき2^N回でループを抜けることができます。

for (int U = (1 << num) - 1; U >= 0; --U) {
  U &= T;
  // ここに処理を書く
}

UをTでマスクしているのがポイントです。 1を引く(デクリメントする)という操作は1になっている一番下位のビットを0にして、それより下のビットをすべて1にする、という操作なので（0のときは除く）、結果的にもれなくかぶりなく列挙できます。

ある集合を部分集合に持つ集合全体を渡る

T⊂UなるUを全部列挙します。

for (int i = T; i < (1 << num); i=(i+1)|T) {
  // ここに処理を書く
}

仕組みは部分集合列挙とほとんど同じです。

要素数kの部分集合全体を渡る

int x, y;
for (int T = (1 << k) - 1; T < (1 << N); x = T & -T, y = T + x, T = (((T & ~y) / x) >> 1) | y) {
  // ここに処理を書く
}

かなりトリッキーですね。詳しくは参考文献の解説を読んでください。

部分集合と和に関するテクニック

2021/01/03 高速ゼータ変換/メビウス変換/畳み込みに関する部分を全面的に書き換えました。

部分集合について和を取るもの
${\displaystyle{ g(U) = \sum_{T \supset U} f(T)}}$ （Uを包含するTについて和を取る: sum over supersets）
${\displaystyle{ g(U) = \sum_{T \subset U} f(T)}}$ （Uに包含されるTについて和を取る: sum over subsets）
逆変換、要素数ごとに符号が反転する
${\displaystyle{ g(U) = \sum_{T \supset U} (-1)^{|U-T|}f(T)}}$ （Uを包含するTについて和を取る）
${\displaystyle{ g(U) = \sum_{T \subset U} (-1)^{|U-T|}f(T)}}$ （Uに包含されるTについて和を取る）

これらの計算は、上述の「ある部分集合の部分集合全体を渡る」テクニックにより、O(3^N)で計算することができます（計算量の証明は二項定理による）。しかし、さらに効率よく計算するアルゴリズムが存在します。これらはしばしば高速ゼータ変換/高速メビウス変換と呼ばれています。

あらかじめdp[T]にf(T)の値を入れておくことにします。 1は

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 0) {
      dp[T] += dp[T | (1 << i)];
    }
  }
}

2は

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 1) {
      dp[T] += dp[T ^ (1 << i)];
    }
  }
}

3は

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 0) {
      dp[T] -= dp[T ^ (1 << i)];
    }
  }
}

4は

for (int i = 0; i < N; ++i) {
  for (int T = 0; T < (1 << N); ++T) {
    if (((T >> i) & 1) == 1) {
      dp[T] -= dp[T ^ (1 << i)];
    }
  }
}

というコードでO(N2^N)でできます。((T >> i) & 1) が0か1か、Tにビットを付け足すのか取り除くのか、演算が+なのか-なのかが違うだけです。また、ここでは加算演算を利用しましたが、可換で結合的で、単位元がある（つまり、可換モノイド）であれば1と2が、さらに逆元が存在（つまり、可換群（アーベル群））すれば3,4を適用できます。

畳み込み

Sの部分集合を引数に取る関数f,gがあったとして、 ${\displaystyle{ (f \circ g)(U) = \sum_{T \subset U} f(T)g(U-T)}}$ をSの各部分集合Uに対して求めたいとします。 (ランク付きの)sum over subsetsと、その逆変換を用いると計算できます（多項式乗算におけるFFT/逆FFTと似た理屈）

fm = ranked_sum_over_subsets(f);
gm = ranked_sum_over_subsets(g);
for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < (1 << N); ++i) {
    for (int j = 0; j <= k; ++j) {
      conv[k][i] += fm[j][i]*gm[k-j][i];
    }
  }
}
result = inverse_ranked_sum_over_subsets(conv);

ランクkのsum over subsetsは ${\displaystyle{ \hat{f}(k,U) = \sum_{T \subset U, |T| = k} f(T)}}$

実際のコードはこんな感じ(dp[popcount(T)][T]にはf(T)の値を入れておく)

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if (((T >> i) & 1) == 1) {
        dp[k][T] += dp[k][T ^ (1 << i)];
      }
    }
  }
}

逆変換は ${\displaystyle{ f(U) = \sum_{T \subset U} (-1)^{|U-T|} \hat{f}(|U|,T)}}$

dp[k][T]に ${\hat{f}(k,T)}$ の値を入れておくと、

for (int k = 0; k <= N; ++k) {
  for (int i = 0; i < N; ++i) {
    for (int T = 0; T < (1 << N); ++T) {
      if (((T >> i) & 1) == 1) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。また、ここでは加算と乗算を用いましたが、他の可換環でも同様の計算を行うことができます。 (max, +)を用いるトロピカル半環の場合は、元の集合が十分小さい場合にはexp(x)(底は適切に取る)に対して(+,×)の畳み込みを行えばできます。

参考文献

popcount: 明日使えないすごいビット演算
最適化版popcount: Hamming weight - Wikipedia
部分集合の列挙: はてなグループの終了日を2020年1月31日(金)に決定しました - はてなの告知
要素数kの部分集合全体の列挙: プログラミングコンテストチャレンジブック第2版 3-2
高速ゼータ変換: はてなグループの終了日を2020年1月31日(金)に決定しました - はてなの告知
高速ゼータ変換、畳込み: 指数時間アルゴリズム入門
高速メビウス（逆）変換: AOJ 2446 Enumeration + 高速メビウス変換まとめ - simezi_tanの日記
部分集合の和、畳込み: ~~http://web.stanford.edu/~rrwill/presentations/subset-conv.pdf リンク切れ~~ http://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

さいごに

ビット演算テクニック Advent Calendar 2016

www.adventar.org

みんな書いて！おねがい！

2016-11-27

CUDAで再帰するとき

GPGPU

CUDAではDynamic Parallelism（GPU側から新たなカーネルを起動する）以外に普通の（シリアル実行の）再帰が可能です(Compute Capability 2.0以上)。しかし、再帰的に関数を呼び出す場合、必要なスタックサイズが静的に割り出せないことがあり、このときはデフォルトのスタックサイズ（私の手元の環境では1024bytes/thread）が使用されます。当然この状態でスタックサイズの制限を超えてスタックを消費すると不幸な目に合います（何も言わずにカーネルの実行が終了してしまうなど）。

解決方法はカーネル呼び出し前にホスト側でスタックサイズを指定してやること。

例えば1スレッドあたり4096バイトのスタックを割り当てたい場合は、ホスト側でカーネル呼び出し前に

cudaDeviceSetLimit(cudaLimitStackSize, 4096);

としてやればよい。

2015-12-12

コミックマーケット89記念！ANSI C89コンパイラを作る【KMCアドベントカレンダー12日目】

コンパイラ

この記事はKMC Advent Calendar 2015 - Adventar12日目の記事です。

宣伝

サークル「京大マイコンクラブ」は、コミックマーケット89で「木曜日　東地区　"モ"　42b」に配置されました！　コミケWebカタログにてサークル情報公開中です https://t.co/BtKSerlV6R #C89WebCatalog
— 京大マイコンクラブ (@KMC_JP) 2015年10月30日

部誌、自作ゲーム、そして今回の本題であるコンパイラのソースコードを置く予定です。 3日目(12月31日)東地区　"モ"　42bです！よろしくお願いします！！*1　

本題

今度のコミックマーケットは89回目です。コミックマーケット89は通称C89と呼ばれています。 C89といえばC言語の最初の標準規格ANSI-C:1989(通称C89, ANSI C89など)ですね！

こうなったらやることは1つです。C89でC89コンパイラを頒布するのです！

出落ち的アイデアですが、もう二度とC89はやってこないので、このチャンスを逃すわけには行かないと部内でC89コンパイラを作るプロジェクトを立てました。

話し合いの結果コンパイラの名称はkmc89になりました。リポジトリはこちらになります。

github.com

規格書の入手

現行のC言語の標準規格ISO/IEC 9899:2011(通称C11)の規格書はISOのウェブサイトからPDFを購入できます。一つ前のC言語の標準規格ISO/IEC 9899:1999(通称C99)は和訳(JIS X3010:2003)をJISCのウェブサイトから閲覧することが出来ます。

しかし、C89の規格書は上記の方法では購入・閲覧等ができません。そこで、規格書のハードコピーをhttp://infostore.saiglobal.com/store/Details.aspx?productID=434093から注文して入手しました。

これは厳密にはANSI-C:1989ではなくISO/IEC 9899:1990(通称C90)の規格書ですが内容としては同じものです。*2

コンパイラを作る

コンパイラを作ると言っても、あまり時間はありません。 C89のすべてを実装するのは流石に無理があるということで、標準ライブラリ、リンカ、プリプロセッサ等は実装せず（既存のリンカ等を用いる）、狭義のコンパイラに絞って実装することにし、オブジェクトコードとしてLLVM IRを生成することにしました。

実際には規格書を注文してから届くまでにかなり時間がかかった（再送もあったので2ヶ月弱）ので届く前から少しずつ書き始めています。

最終的にセルフホスティングを目指したいのでコンパイラもC89で記述することになり、それにあたってC言語で貧弱な部分をライブラリとして実装する作業が行われています。 C++のstd::vectorやstd::stringに当たるものをC言語で実装しようとしています。KMCアドベントカレンダー9日目の記事template in C - KMC活動ブログも参照してください。

ある程度そろえばあとはflex、bison、LLVM等の便利なライブラリを使うことで12/31のコミックマーケット3日目に間に合うように実装できるはず！！ということで現在鋭意実装中です！応援よろしくお願いします！

このあたりのコンパイラの仕組みの話は今度のKMCの部誌、独習KMC vol.8に少し書いてあるので気になる人はぜひ読んでください。

今後の展望

今回はコンパイラのコア部分に絞って実装することでなんとか間に合わせる作戦ですが、C89の次のコミケはC90ですね。この時までにC89≒C90のコンパイラの周りの部分の実装を進めて完全版を配布したいと考えています。*3

また、C90の部誌ではコンパイラの実装や開発環境等についての苦労話や新たに得た知見なども書きたいと思っています。ご期待ください。

では早速この記事を書き上げたらコンパイラの実装を進めることにしたいと思います。応援よろしくお願いします！

github.com

コミケは3日目(12月31日)東地区　"モ"　42bです！よろしくお願いします！！

絶対C89にC89コンパイラを間に合わせます！間に合わなかったら木の下に埋めてもらって構わないよ！

次回予告

明日のKMC Advent Calendar 2015 - Adventar 13日目の記事はnonamea774(nona7)さんによる「[データ削除済]」です。お楽しみに。

www.adventar.org

*1:残念ながら当日私はいません

*2:C89の文章に章立てを追加したものになっている

*3:C99の時にC99コンパイラを作るかどうかは未定です。

2015-12-04

高速文字列処理ライブラリを作った

ネタ競技プログラミング

この記事はポエムアドベントカレンダー4日目の記事です。

www.adventar.org

大量の文字列データを扱うことの多くなった現代において、文字列処理ライブラリの高速化は重要である。しかしながら、個人レベルで汎用的かつ高速な文字列処理ライブラリを作成することは難しい。今回は汎用性を少し下げることにより圧倒的な高速化をした文字列処理ライブラリ「A」を制作した。

ソースコード

gist.github.com

ライブラリの仕様

文字列の制約

すべての文字がAで構成されていること

制約を満たす文字列の例

AAAAAA

AAAAAAAAAA

制約を満たさない文字列の例

aaa

ABCDE

制約を満たさない文字列を用いた場合、正しくない結果を得る可能性がある。

利用方法

ライブラリ中では専用の効率的なデータ構造により文字列を管理する。そのため、利用するにはstd::stringから変換処理を行う必要がある。

A a(std::string(100, 'A'));

専用のデータ構造からstd::stringに変換するにはto_cppstringメンバ関数を用いる。

実装されているアルゴリズム

検索 `find(A_1, A_2)`

文字列 $A_1$ の部分文字列に正規表現 $A_2$ に一致するものが存在すればその位置（複数あればそのうち最も先頭に近いもの）を返す。計算量: $O(1)$

置換 `replace(A_1, A_2, A_3)`

文字列 $A_1$ の部分文字列に正規表現 $A_2$ に一致するものが存在すればそれをすべて $A_3$ で置き換える。計算量: $O(1)$

編集距離 `edit_distance(A_1, A_2)`

2つの文字列の編集距離を計算する。計算量: $O(1)$

最長共通部分列 `longest_common_subsequence(A_1, A_2)`

2つの文字列の共通部分列で最も長いものを返す。計算量: $O(1)$

最長回文 `longest_palindrome(A)`

回文になっている部分文字列のうち最も長いものとその始点を返す。計算量: $O(1)$

今後の課題

さらなるライブラリの充実
Bのみからなる文字列を高速に扱えるライブラリの開発

Bit scan reverse (bsr)

bsr/lzcnt命令

スカラーに変換してbsrで処理 (naive)

popcountに帰着 (popcount)

シャッフル命令vpshufbを使った最適化 (rev popcount)

ギャザー命令を使ったテーブル参照 (table gather)

vcvtdq2ps命令による浮動小数点数への変換 (cvtfloat)

vpshufbを使ったテーブル参照 (table)

環境

結果

8ビット

16ビット

32ビット

64ビット

来るAVX512時代には

参考文献

ビット単位のunpack

SSE

pdep

pclmulqdq

さいごに

行える操作

実装

応用

bit matrix

参考文献

はじめに

表現方法

基本操作

空集合、全体集合

和集合

共通部分

補集合

部分集合全体を渡る

包含関係を判定する

要素数を数える

応用操作

ある部分集合の部分集合全体を渡る

ある集合を部分集合に持つ集合全体を渡る

要素数kの部分集合全体を渡る

部分集合と和に関するテクニック

畳み込み

参考文献

さいごに

宣伝

本題

規格書の入手

コンパイラを作る

今後の展望

次回予告

ライブラリの仕様

文字列の制約

制約を満たす文字列の例

制約を満たさない文字列の例

利用方法

実装されているアルゴリズム

検索 find(A_1, A_2)

置換 replace(A_1, A_2, A_3)

編集距離 edit_distance(A_1, A_2)

最長共通部分列 longest_common_subsequence(A_1, A_2)

最長回文 longest_palindrome(A)

今後の課題

検索 `find(A_1, A_2)`

置換 `replace(A_1, A_2, A_3)`

編集距離 `edit_distance(A_1, A_2)`

最長共通部分列 `longest_common_subsequence(A_1, A_2)`

最長回文 `longest_palindrome(A)`