prime's diary

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

ビット演算マニアのためのAVX512入門 【KMCアドベントカレンダー 22日目】

この記事はKMCアドベントカレンダー22日目の記事です。

adventar.org

この記事では全国一億三千万のビット演算マニアのために、AVX512命令セットからビット演算に使えそうなものを紹介します。

AVX512の基本

AVX512とは

Advanced Vector Extensions(AVX)というx86/x64のベクトル演算拡張命令を512ビット幅に拡張したものです。 普通のCPUではSkylake-SP/X/Wから使えるようになりました。 512bitのデータを一括で処理できます。ビット演算なら512並列ですよ512並列。

AVX512になって変わった主なこと

AVX512の命令セット

AVX512は複数の要素から構成されており、プロセッサによってサポートの状況が異なります。

以下に書くサブセットとその説明を書きます。

  • AVX512-F (Foundation) 基本命令セット。AVX512をサポートするプロセッサは必ずこれをサポートする。
  • AVX512-CD (Conflict Detection) データの衝突検出。Hist処理とかがベクトル化できる。 vpconflictd - Qiita
  • AVX512-DQ (Doubleword and Quadword) 32bit/64bit単位での処理。
  • AVX512-BW (Byte and Word) 8bit/16bit単位での処理。
  • AVX512-VL (Vector Length) 128bit(XMM)/256bit(YMM)レジスタを操作する。

最新のCPU(Skylake-SP/X/W)でもまだ使えない拡張もたくさんあります。

  • AVX512-PF (Pre Fetch) プリフェッチ命令。
  • AVX512-ER (Exponential and Reciprocal) 指数関数や逆数関数。
  • AVX512-IFMA52 (Integer Fused Multiply Add) 整数の積和算を倍精度の積和算演算器を借りて52ビット精度で実行する。
  • AVX512-VBMI (Vector Byte Manipulation Instructions) BWに入っていないbyte操作の命令。
  • AVX512-VNNI (Vector Neural Network Instructions) ニューラルネットワーク処理用の整数演算命令。16bit整数同士の積を32bit整数と足す。
  • AVX512-4VNNIW (Vector Neural Network Instructions Word variable precision) ニューラルネットワーク処理用の整数演算命令。16bit整数同士の積を32bit整数と足すのを4回繰り返す。
  • AVX512-4FMAPS (Fused Multiply Accumulation Packed Single precision) 単精度の内積を4回繰り返す。これもニューラルネットワーク処理用。
  • AVX512-VPOPCNTDQ (Vector POP CouNT) ベクトル化されたpopcount。
  • AVX512-VBMI2 (Vector Byte Manipulation Instructions2) BWにもVBMIにも入っていないbyte操作の命令。
  • AVX512-BITARG (不明) ビット操作関連の命令。

マスク演算

AVX512では比較演算の結果などはマスクレジスタと呼ばれる特殊なレジスタに格納されます。

AVX512では命令にマスクとしてマスクレジスタを指定することができます。

VPANDQ (and演算, 64bit単位のマスク)を例にして疑似コードで表すと、

for j in range(0, 8):
  i = j*64
  if mask[j]:
    dest[i+63:i] = src[i+63:i]
  else:
    dest[i+63:i] = a[i+63:i] AND b[i+63:i]

このようにマスクでビットが立っている場所を演算せずに値を残すことが可能です。また、値を残す代わりに0にすることもできます。

intrinsicでは

dest = _mm512_and_epi64(a, b); // マスク無し
dest = _mm512_mask_and_epi64(src, mask, a, b); // マスクあり(マスクされたらsrcを残す)
dest = _mm512_maskz_and_epi64(mask, a, b); // マスクあり(マスクされたら0)

のように使います。

マスク同士の演算

このマスクレジスタSIMDレジスタ(XMM/YMM/ZMMレジスタ)とも通常のレジスタとも違うので、通常のレジスタとマスクレジスタとの間で値をやりとりするのにもコストがかかります。 そこで、いくつかの演算(and, or, addなど)はマスクレジスタ同士で行えるようになっています。

命令の解説

AVX512-F

VPTERNLOG{D,Q}

ビット単位の任意の3項演算ができます。演算内容は即値で与えます。 A & (B | ~C)なら

A B C A&(B|~C)
1 1 1 1
1 1 0 1
1 0 1 0
1 0 0 1
0 1 1 0
0 1 0 0
0 0 1 0
0 0 0 0

なので即値として0b11010000 = 0xD0を与えてやれば良いです。 intrinsicは_mm512_ternarylogic_epi64(A, B, C, 0xD0)のように使います。 これを使えば他のビット論理演算命令いらないじゃんとなりそうですが、怪しいページによるとスループットが2命令/サイクルとandやorなどの3命令/サイクルに比べると劣るので適材適所で使うことになりそうです。

VP{AND,OR,XOR,ANDNOT}{D,Q}

and, orなどの演算です。マスクが32bit/64bit単位で指定できます。

VPS{L,R}{L,A}{,V}{W,D,Q}

16bit/32bit/64bit単位のビット{論理,算術}シフトです。V付きはシフト量もベクトルで指定することで、各データごとに異なるシフト量を指定できます。

VPRO{L,R}{,V}{D,Q}

待望のローテート命令です。でも32bit/64bitだけ。

AVX512-CD

VPLZCNT{D,Q}

leading zero countです。上位ビットから見て最初に1が立っているビットの位置を調べる。整数→浮動小数の変換で前からほぼ同じことができましたが、キャストとか0の特別扱いとか諸々せずに済むようになりました。 ところでなんでConflict Detectionに入っているんでしょうか。

AVX512-VBMI

VPERMI2B

各バイトの下位7bitをインデックスにして、2つのZMMレジスタをくっつけた128bytesの中から表引きします。 これまでの128bit境界に捕われていたり、即値で指定しないといけなかったりしたシャッフル命令とはおさらばできます。 とはいえスループットやレイテンシが大きかったら従来のシャッフル命令も併用しないといけなさそうですね。

VPMULTISHIFTQB

各64bit要素に対して各バイトの分だけ右シフトしたデータを返します。

for i in range(0, 8):
  q = i*64
  for j in range(0, 8):
    tmp8 = 0
    ctrl = src1[q+j*8+7:q+j*8] & 0x3F
    for l in range(0, 8):
      tmp8[l] = src2[q+((ctrl+l)&0x3F)]
    dst[q+j*8+7:q*j*8] = tmp[7:0]

AVX512-VBMI2

VPSH{L,R}D{,V}{W,D,Q}

2つの値をつなげてシフトします。多倍長シフトとかに使えますね。V付きの命令はシフト量もベクトルで指定することで、各データごとに異なるシフト量を指定できます。

AVX512-POPCNT

VPOPCNT{D,Q}

32bit/64bit単位でのpopulation count(1になっているビットの数を数える)です。

AVX512-BITARG

VPOPCNT{B,W}

8bit/16bit単位でのpopcountです。なぜかVPOPCNT{D,Q}とは別の拡張になっています。

VPSHUFBITQMB

謎命令。各64bit要素に対して、各バイトをインデックスとして64bitデータの指定したビットを取ってきてマスクにします。疑似コードを見るのが早いと思います。

for i in range(0, 8): # qword
  for j in range(0, 8): # byte
    m = src2[i*64+j*8+7:i*64+j*8] & 0x3F
    res[i*8+j] = src1[i*64+m]
return res

AVX512以外の拡張

VPCLMULQDQ

繰り上がりなしの64bit乗算を並列で行います。ビット演算的にはCRCの計算とかpclmulqdqを用いたビット単位のunpack【ビット演算テクニック Advent Calendar 2016 11日目】 - prime's diaryとかに使えそうですね。

便利Intrinsic

_mm512_reduce_{or,and}_epi{32,64}

32bit/64bitデータ16/8個のOR/ANDを取ります。

_mm512_reduce_and_epi64なら、

reduced[63:0] = 0xFFFFFFFFFFFFFFFF
for j in range(0, 8):
    i = j*64
    reduced[63:0] = reduced[63:0] AND a[i+63:i]
return reduced[63:0]

まとめ

AVX512にはどう使うかよくわからない便利そうな命令がたくさん追加されています。

使えるCPUを手に入れた暁にはうまく使ってビット演算ライフを満喫しましょう!

明日のKMCアドベントカレンダーはwass88さんによる「役に立たないこと」の予定です。お楽しみに!

adventar.org

ICPC 2017 Tsukuba Asia Regional 参加記

2017/12/16-18に行われたICPC 2017 Tsukuba Asia RegionalにチームPrimeDragonとしてamano, nikuttoと参加してきました。

コンテスト前日まで

15日

生活リズム的に16日の朝早くから出られそうにないので徹夜することにした。Redbull飲んで徹夜する。

16日

徹夜のおかげで無事新幹線に乗れる。

つくばエクスプレス乗車。(京都大阪間に比べると)値段高いですね…

中学校の知人とばったり出会う。偶然ではなくてICPCに出てるらしい。

JAXA見学。管制システムかっこいい…

リハーサルではコンテスト時の段取り(誰が最初にテンプレートとか書くのかとか)を確認、コンテストシステムに慣れる。 プログラムが実行されるサーバーと手元の速度差を確認。(ほとんど違いはなかった) 印刷やClarを投げる練習もした。

筑波学院大学でWelcome Party。 京大の伝統(?)でチーム紹介は短めに。 料理がだいたい茶色系だったのはなんでだろう…おいしかったけど

歩いて会場から宿まで歩いた。さむい。

ちょうどARCがあったので参加。3完できたのでまあまあ良い結果。

寝る。

コンテスト当日(17日)

緊張してあんまり寝られなかったけどおかげで寝坊せず起きられた。

朝食はつくばカピオの近くのCASAってところで食べた。natsugiriさんとかがいた。

コンテスト本番

記憶が曖昧なところがあるので時系列おかしいかも知れないです。

コンテスト開始。

まず、Xmodmapの設定(Caps LockをCtrlにする)、.vimrcやテンプレートを書く。 その間にAやBを読んでもらう。

書けた時点でAが読めて解法もわかったらしいので教えてもらって書く。AC。 Bが解けるらしいので交代する。

Cを読む。特に難しいところはなさそう。Bを解いている間に実装を詰める。 BがACしたので交代してCの実装をする。サンプルが通らなくて焦るが、インデックスがずれているだけだった。 2個もずれててこんなずれるはずないだろと思っていたが、サンプルは通るので投げる。AC。

この時点で10位ぐらいだった気がする。出だしとしてはいい感じ。

ここからは難易度順に並んでいないので簡単そうな問題をさがす。 Fが始点と終点からDijkstraすれば良さそうなので解けると主張するが、それだけだと経路が長くなるか判定できないと言われてたしかにとなる もうちょっと考えると最短経路上の辺で作った部分グラフ(DAGになる)で橋になっているかがわかれば良いことがわかる。 コードを書いてサンプルも通ったので投げる。WA。コード印刷するが原因わからず。

その間にGで実装軽そうな方針が立ったので書いてもらう。AC。

まだわからないので他の問題を考えていた。Iが解けるらしいので書いてもらう。AC。

ようやくFで橋の判定アルゴリズムが間違っていたことに気づき(DAGという条件で軽く書けると思っていたが、嘘解法だった)、ライブラリの橋検出を写す。AC。

この時点で5位くらいだったかな?かなり良い順位。このまま比較的解いている人の多いEを通せば良さそう。

しかし、EでDP解(これが正解)ではなく貪欲法解をずっと考察していていくつも嘘解法を作ってはだめとなっていた。

結局Eは通せずにコンテスト終了。

コンテスト後

解説&順位発表&表彰。D実装つらそう…こんなんコンテスト中に通せる気がしない… 結果は6完で7位。

夕食&企業ブース見学。料理うまい。前日に比べて彩りもあってよかった。 企業ブースを一通り回る。 レコチョクのパズルを解いて景品(モバブ、アメなど)をもらう。このことが後々生死を分かつことになる。 あとドワンゴのWA投げでトランプをもらった。

企業賞みたいなのでレコチョクからワイヤレスイヤホンをもらう。結構良いやつっぽい。

E通してればWF行けたんでは?とか一応まだWFの可能性ありそう、とかのコメントももらった。まさかWFの可能性がある順位取れると思ってなかったので嬉しいが、余計Eを解けなかったことが辛くなる。

昨日と同じく歩いて宿に帰る。つくばカピオからは近い。

この日はなかなか眠れなかったのでオセロでモンテカルロ木探索を実装した。合法手生成でバグってうまく動かなかった。

ホテルの自販機でDrPepperが売っていたので人生で初めて飲む。謎の味がした。 よく考えたら寝れなくて困ってるのにカフェイン入り飲料を飲むのは間違っていた。

結局朝6時頃に就寝。

コンテスト翌日(18日)

企業見学のためつくばから東京に移動しようとするが、つくばエクスプレスが30分ぐらい遅れて遅刻する。

企業見学前半はLINE。以前ISUCONで行ったことがあるが、当時は渋谷にあったので新しい新宿のオフィスに来るのは初めて。 窓からの眺めが良かった。あと叙々苑の焼肉弁当が美味しかった。

企業見学後半はPFN。ICPC/JOI出身者がすごく多いらしい。ロボットの制御もやっていたのは知らなかった。

解散。このあとせっかく東京に着たのでオフ会して終電で帰った。途中で持ってきていたモバブの充電が切れたが、レコチョクからもらったモバブのおかげで生きながらえた。

感想

E通したかった…

オセロの着手可能位置生成 【ビット演算テクニック 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にした意味がほとんどありません。

高速化する

2018/03/31追記

このプログラムが遅い最も大きな要因は、以下で説明するレイテンシの隠蔽失敗によるものではなく、再帰的な分岐により演算器が遊んでしまうことによるものです。 実際、Shared Memoryを少し節約しつつ、1スレッドで1局面読むようにすると下記の結果の9秒より速く計算することが出来ます。

追記終わり

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の高速な実装を参照。

CUDA Othello Solver · 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にあります。

参考文献