prime's diary

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

オセロの着手可能位置生成 【ビット演算テクニック Advent Calendar 2016 23日目】

はじめに

この記事はビット演算テクニック Advent Calendar 2016の23日目の記事です。

www.adventar.org

オセロとビット演算

オセロは盤面サイズが8x8=64マスなので、64bit変数を2つ用意し、1つめを自分の石(あるいは黒石)の位置を、2つめを相手の石(あるいは白石)を表すために使うと合計128ビットで自然に表現できます。*1
自然に表現できるだけでなく、盤面の回転や反転などもDelta swapとその応用 【ビット演算テクニック Advent Calendar 2016 3日目】 - prime's diaryで紹介したような方法で効率的に行なえます。

着手可能位置の生成はタテ・ヨコ・ナナメの一直線上をpextなどで取ってきて表引きをすることにより、ある程度の速度で可能ですが、(38=6561なのでL1キャッシュに乗る程度の大きさとはいえ)メモリアクセスが必要になってしまいそれなりの時間がかかります。

以下、ビットボードには以下のようにマス目とビット列の対応が付いているとします(0が最下位ビット)。

A B C D E F G H
1 0 1 2 3 4 5 6 7
2 8 9 10 11 12 13 14 15
3 16 17 18 19 20 21 22 23
4 24 25 26 27 28 29 30 31
5 32 33 34 35 36 37 38 39
6 40 41 42 43 44 45 46 47
7 48 49 50 51 52 53 54 55
8 56 57 58 59 60 61 62 63

ひっくり返る石の計算

参考文献の方法を使っています。 参考文献のOpenCL風のコードがわかりやすいのでそれを元に動作を解説します。

ulong flip(int pos, ulong P, ulong O)
{
  ulong4 flipped, OM, outflank, mask;

  OM.x = O;
  OM.yzw = O & 0x7e7e7e7e7e7e7e7eUL; // 1
  mask = (ulong4) (0x0080808080808080UL, 0x7f00000000000000UL, 0x0102040810204000UL, 0x0040201008040201UL) >> (63 - pos); // 2
  outflank = (0x8000000000000000UL >> clz(~OM & mask)) & P; // 3
  flipped  = (-outflank * 2) & mask; // 4
  mask = (ulong4) (0x0101010101010100UL, 0x00000000000000feUL, 0x0002040810204080UL, 0x8040201008040200UL) << pos;
  outflank = mask & ((OM | ~mask) + 1) & P; // 5
  flipped |= (outflank - (outflank != 0)) & mask;
  return flipped.x | flipped.y | flipped.z | flipped.w;
}
  1. 盤面の端を超えて石が繋がっていても無視する
  2. posから見て上、左、右上、左上のマス目をマスクするビット列を計算する
  3. posからみて相手の石が途切れたところをclzで求め、自分の石の位置とandを取ることでひっくり返せるかどうかがわかる。outflankの各要素は高々1ビットが立つ。
  4. -outflank * 2でoutflankより上位のビットが全部1になる。これとマスクを取ることでひっくり返る石の位置がわかる。
  5. 後半部分は下位ビットから見ていって相手の石が途切れたところを探します。これは関係ない場所を1埋めしてから+1すればよいです。あとは単に前半部分の逆です。

実際のAVX2をゴリゴリ使った実装はissen/move_generator.cpp at 72f450256878094ffe90b75f8674599e6869c238 · primenumber/issen · GitHubにあります。

このコードでは各方向に対するひっくり返る石をビット演算で求めることができ、AVX2を使うと256bit一気に操作でき高速です。ただ、clzの並列版はないので、立っている一番上のビット位置を頑張って求めます。

Bit scan reverseをAVXで並列化する - prime's diaryと違い、1<<bsrを求める関係で浮動小数点数への変換より下位ビット埋め+vpshufb+最下位ビット抽出のほうが速いです。

立っている一番上のビットを求める実装例: issen/move_generator.cpp at 72f450256878094ffe90b75f8674599e6869c238 · primenumber/issen · GitHub

着手可能位置の計算

着手可能位置(mobility)の計算をビットボードで行う方法は

d.hatena.ne.jp

で言及されているビットシフトを繰り返す方法が知られているようです。今回はそれより高速な方法を解説します。

この方法では、石を置くと左に向かって相手の石をひっくり返せるような位置を一気に列挙します。
これを盤面を回転・反転させて8方向に対して行うことで着手可能位置を計算することができます。 1列分の着手可能位置を求めるコードは以下の通り。

// p: 自分の石, o: 相手の石
uint8_t mobility_line(uint8_t p, uint8_t o) {
  uint8_t p1 = p << 1;
  return ~(p1 | o) & (p1 + o);
}

このプログラムの動作を具体例を元に解説します。 例えば-oop-opo(p: 自分の石, o: 相手の石, -: 空きマス)のとき、

7 6 5 4 3 2 1 0
p 0 0 0 1 0 0 1 0
o 0 1 1 0 0 1 0 1

pを1ビット左シフトして、

7 6 5 4 3 2 1 0
p1 0 0 1 0 0 1 0 0
o 0 1 1 0 0 1 0 1

2つを足し算すると、自分の石より上位にある相手の石が途切れる位置(つまり着手可能位置) または 孤立した相手の石 または 孤立した自分の石の一つ上位 のビットが立ちます。

7 6 5 4 3 2 1 0
p1+o 1 0 0 0 1 0 0 1

~(p1|o)でマスクを取って、

7 6 5 4 3 2 1 0
res 1 0 0 0 1 0 0 0

となり着手可能位置が計算できました。 これを64ビット*2の盤面全体に対して一気にやろうとすると、8ビット単位のシフトや足し算が必要になります。

うまくマスクを取って64ビット演算に帰着してもいいですが、SSE/AVX2には8ビット単位の加算が用意されているので、それを使うと良いでしょう。1ビットの左シフトはp+pで実現できます。 SSEなら2並列、AVX2なら4並列で着手可能位置を求められるので、8方向に対する計算をかなり高速化することができます。

全体をSSE/AVX2で書いたものがこちらになります。

issen/movable_generator.cpp at master · primenumber/issen · GitHub

参考文献

ひっくり返る石の計算方法 Reversi - bitboard and GPU computing

すでに知られている着手可能位置の計算方法 リバーシプログラム - bitboard による合法手生成 - プログラミング練習帳&備忘録

*1:1つ目で両方の石の位置を、2つ目で相手の石の位置を表す流儀もある

GPGPUで爆速でオセロを解く 【KMC Advent Calendar 2016 19日目】

はじめに

この記事はKMC Advent Calendar 2016の19日目の記事です。

www.adventar.org

GPGPUとは

GPGPU(General-purpose computing on graphics processing units; GPUによる汎用計算)とは、GPUの演算資源を画像処理以外の目的に応用する技術のことである。

(出典: GPGPU - Wikipedia )

GPUは並列演算に特化しており、大量のスレッドをハードウェアで処理でき、汎用CPUに比べて演算処理性能が高いという特徴があるので、これを用いて演算を高速に行おうというのがGPGPUのアイデアです。

ただ、良いことばかりではなく、汎用CPUに比べてプログラミング上の制約が大きいなど、性能を出し切るのはかなり難しいです。

GPGPUでオセロを解く

今回はGPGPUを用いてオセロの終盤の完全読み*1を高速化します。

オセロを解く際に基本となるのがアルファベータ法です。

詳しい解説は他に譲るとして、このアルゴリズムの肝はすでに計算した部分の結果を利用して枝刈りを行うことです。そのため、並列化が難しいアルゴリズムになっています。 今回はアルゴリズムレベルでの並列化は後回しにして、たくさんの局面を並列で完全読みすることを考えます。

各スレッドに一つの局面を割り当てて完全読みを行うのが素直な実装ですね。

アルファベータ法は再帰を用いると自然に書けるので再帰で実装してみます。 今回はCUDAで実装します(したがって、この記事では以後GPUNVIDIAGPUと仮定します)。CUDAで再帰するときは必要なスタック数が静的に見積もれないので必要に応じて動的にスタックサイズを指定してやります(今回は4096バイト/スレッドにしました)。*2 高速化のためビット演算を用いてひっくり返る石の位置を計算しています。詳しくはReversi - bitboard and GPU computingを参照してください。

CUDAカーネル部分だけ示します(CPU側のコードはスタックサイズの指定を除いて後で示す改良バージョンと同じです)。

constexpr int threadsPerBlock = 128;

using ull = unsigned long long;

struct Board {
  ull player;
  ull opponent;
  __host__ __device__ Board() = default;
  __host__ __device__ Board(ull p, ull o) : player(p), opponent(o) {}
  __host__ __device__ Board(const Board &) = default;
  __host__ __device__ Board(Board &&) = default;
  __host__ __device__ Board &operator=(const Board &) = default;
  __host__ __device__ Board &operator=(Board &&) = default;
};

__shared__ int count[threadsPerBlock];

__device__ char score(const Board &bd) {
  char pnum = __popcll(bd.player);
  char onum = __popcll(bd.opponent);
  if (pnum == onum) return 0;
  if (pnum > onum) return 64 - 2*onum;
  else return 2*pnum - 64;
}
__constant__ ull mask1[4] = {
  0x0080808080808080ULL,
  0x7f00000000000000ULL,
  0x0102040810204000ULL,
  0x0040201008040201ULL
};
__constant__ ull mask2[4] = {
  0x0101010101010100ULL,
  0x00000000000000feULL,
  0x0002040810204080ULL,
  0x8040201008040200ULL
};

__device__ ull flip(const Board &bd, int pos, int index) {
  ull om = bd.opponent;
  if (index) om &= 0x7E7E7E7E7E7E7E7EULL;
  ull mask = mask1[index] >> (63 - pos);
  ull outflank = (0x8000000000000000ULL >> __clzll(~om & mask)) & bd.player;
  ull flipped = (-outflank * 2) & mask;
  mask = mask2[index] << pos;
  outflank = mask & ((om | ~mask) + 1) & bd.player;
  flipped |= (outflank - (outflank != 0)) & mask;
  return flipped;
}

__device__ ull flip_all(const Board &bd, int pos) {
  return flip(bd, pos, 0) | flip(bd, pos, 1) | flip(bd, pos, 2) | flip(bd, pos, 3);
}

__device__ Board move_pass(const Board &bd) {
  return Board(bd.opponent, bd.player);
}

__device__ char alpha_beta(const Board &bd, char alpha, char beta, bool passed) {
  ++count[threadIdx.x];
  ull puttable = ~(bd.player | bd.opponent);
  if (puttable == 0) return score(bd);
  bool pass = true;
  for (; puttable;) {
    ull bit = puttable & -puttable;
    puttable ^= bit;
    int pos = __popcll(bit-1);
    ull flipped = flip_all(bd, pos);
    if (flipped) {
      pass = false;
      Board next(bd.opponent ^ flipped, (bd.player ^ flipped) | bit);
      alpha = max(alpha, -alpha_beta(next, -beta, -alpha, false));
      if (alpha >= beta) return alpha;
    }
  }
  if (pass) {
    if (passed) {
      return score(bd);
    } else {
      return -alpha_beta(move_pass(bd), -beta, -alpha, true);
    }
  }
  return alpha;
}

__global__ void search_noordering(const Board *bd_ary, int *res_ary, int *nodes_count, const size_t size) {
  size_t index = threadIdx.x + blockIdx.x * blockDim.x;
  while(index < size) {
    count[threadIdx.x] = 0; 
    res_ary[index] = alpha_beta(bd_ary[index], -64, 64, false);
    nodes_count[index] = count[threadIdx.x];
    index += blockDim.x * gridDim.x;
  }
}

GTX1080(ZOTAC Geforce GTX 1080 AMP Edition, Base 1683MHz/Boost 1822MHz)でGGSの棋譜から生成した10石空きの問題1091780局面の完全読みを実行してみると5分48秒かかりました。

一方、CPUで各種高速化を施したソルバー*3ではCore i7-6700K上で同じ処理を行うのに1分13秒しかかかりませんでした。これではGPGPUにした意味がほとんどありません。

高速化する

GPUはSM(Streaming Multiprocessor)という単位でワープを管理していて、GTX 1080(や他のCompute Capability 6.1のGPU)では同時に4つのワープを実行できます。しかし、実際にはSMは4つよりも多くのワープを管理でき、実行可能なワープを選んで実行することで、命令やメモリの読み書きのレイテンシを隠蔽しています。逆に言うと、同時にたくさんのワープを管理できなければレイテンシを隠蔽できず、多くの時間演算器は遊んでしまうことになります。

今回の場合、再帰レジスタをかなりたくさん消費してしまっているので、一つのSMで持てるワープは数個しかありません。

そこで、Shared Memoryに再帰のスタックの内容を移し、再帰をループに直すことにします。 こうした場合今度はShared Memoryの使用量がSMの持てるワープ数を制限してしまうので、4スレッドで1局面を読むことで使用量を減らし、その代わり4スレッドをSIMD的に使うことでひっくり返る石の場所の計算を高速化します。

また、一つの局面を読み終わったらすぐ次の局面を読めるように変更しました。これにより、ワープの中で最も遅いスレッドが律速になるのを軽減できます。

その他の細かい改良を含めた結果、実行時間は最終的に9秒になりました。CPUの8倍高速です。

ソースコード及びベンチデータ

ソースコードは以下の通り。 問題ファイルの読み込みはBase81という1文字で4マス分表すフォーマットです。詳しくはGitHub - primenumber/issen: オセロAIの高速な実装を参照。

othello_solver.cu · GitHub

ベンチマークに使用した入力ファイルは次のリンクからダウンロードできます

prob10.b81 - Google ドライブ

今後の課題

  • 枝刈りの成功率を上げるため、次の手を読む順番を良さそうな順番になるように並べ替える。 並べ替えのためにメモリをかなり消費するので、グローバルメモリ上に置かないと実行効率が上がらないと予想される。

  • 多数の局面を並列に読むのではなく、アルゴリズムを工夫して一つの局面を並列に読めるようにする。 あわよくばCPUで読んだときより速くならないだろうか…

参考文献

Reversi - bitboard and GPU computing 五石空きをminimax法で完全読みするのをGPUにやらせてみたがあまり速くなかったという結果。 ビット演算部分はここを参考にした。

GPU Acceleration for Board Games | NVIDIA Developer この記事と似ていて、並列にオセロのAIを走らせているらしい。Windowsを持ってないので動作を確かめられない。

*1:その局面から両者が最善を尽くしたときに、勝敗及び石差がどうなるかを判定すること。

*2:http://primenumber.hatenadiary.jp/entry/2016/11/27/031429

*3:GitHub - primenumber/issen: オセロAIの高速な実装

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*231byte = 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 化(その2)

を参照してください。

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が合計231回実行されるように、8ビット幅なら226回、64ビット幅なら229回ループを回します。

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 化(その2)

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

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を使っていては遅いのでintelx86の拡張命令セットである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に相当するので)

{\displaystyle 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多項式の自乗を考えてみます。

{\displaystyle (\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それぞれ230回(≒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

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で任意の並べ替えを実現できます。(2nビットなら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);
}

他の例はchessprogramming - Flipping Mirroring and Rotatingにあります。

参考文献

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

はじめに

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

www.adventar.org

の1日目の記事です。

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

表現方法

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

Tc = {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行目で2nビットごとの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としたとき2N回でループを抜けることができます。

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) {
  // ここに処理を書く
}

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

高速ゼータ変換

Sの部分集合を引数に取る関数fがあったとして、 {\displaystyle{ g(U) = \sum_{T \supset U} f(T)}} をSの各部分集合Uに対して求めたいとします。 各Uについて個別に求めれば二項定理により全体の計算量はO(3N)となりますが、Tの取りうる数が2N通りしかないので無駄が多いです。 これを効率よく計算する高速ゼータ変換というアルゴリズムがあって、 あらかじめdp[T]にf(T)の値を入れておくと、

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)];
    }
  }
}

というコードでO(N2N)でできます。

高速メビウス変換・逆変換

Sの部分集合を引数に取る関数fがあったとして、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} f(T)}}] をメビウス変換(ゼータ変換と包含の向きが逆)、 [tex:{\displaystyle{ g(U) = \sum{T \subset U} (-1)^{|U-T|}f(T)}}] をメビウス逆変換といいます。(集合の引き算は差集合。) その名の通りメビウス変換して逆変換すると元に戻ります。

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

逆変換は差集合の要素数に依存する符号があって複雑そうですが、実はこちらも似たようなコードでできます。(符号が逆になっただけ)

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

畳み込み

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

fm = ranked_mobius(f);
gm = ranked_mobius(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 = ranked_mobius_inv(conv);

ランクkのメビウス変換は{\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) {
        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) {
        dp[k][T] -= dp[k][T ^ (1 << i)];
      }
    }
  }
}

でdp[popcount(T)][T]に結果が代入されるように計算できる。(実は畳み込みしなくてももともと結果が入っている)

参考文献

さいごに

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

www.adventar.org

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

CUDAで再帰するとき

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

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

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

cudaDeviceSetLimit(cudaLimitStackSize, 4096);

としてやればよい。