prime's diary

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

「サイゼリヤで1000円あれば最大何kcal摂れるのか」を全探索で解いてみた

概要

qiita.com

この記事や後述の先行研究に触発されて、「サイゼリヤで1000円あれば最大何kcal摂れるのか」を最も単純に全探索で解いてみました。

先行研究

先人たちがいろいろな方法で解決を試みています。

qiita.com

qiita.com

qiita.com

qiita.com

qiita.com

qiita.com

qiita.com

anqou.net

実装

あらかじめ実装したものがこちらになります。

github.com

メニューデータは

github.com

にあるmenu.csvをダウンロードしてきました(make時に自動でダウンロードされます)

実装の注意点としては、探索の途中で合計金額があらかじめ指定した金額(1000円)を超えるような組み合わせは調べないようにしています。*1

実際に計算してみます。

$ ./solve 1000
Loading menu...
Solving...
Solved!
ラージライス: 219 yen, 454 kcal
アーリオ・オーリオ(Wサイズ): 574 yen, 1120 kcal
ポテトのグリル: 199 yen, 366 kcal
Price: 992 yen
Calorie: 1940 kcal

先行研究の結果と一致しているので正しく計算できているようです。

私の手元のPCでは0.4秒弱で計算が終わりました。全探索しても大したことなかったですね。

つまり数々のアルゴリズムはすべて無駄だったわけです。

…というわけではなく、この方法だと上限金額を大きくするとあっという間に計算時間が爆発します。

試しに1500円で計算させたところ、30秒弱かかりました。1000円から1700円程度までの実行時間を元に予測すると、おおよそ指数的に計算時間が増えて、2000円だと30分程度、3000円だと100日程度かかってしまう計算になります。

手法によっては10000円でも現実的な時間で解くことができているので、サイゼリヤで豪遊したい場合はアルゴリズムを学ぶことが必須なことがわかります。

感想

もっと計算時間がかかるかと思っていたので拍子抜けしています。本当はGPGPUで殴って楽しい高速化ライフ!という記事にしたかったのですが…

*1:この工夫を入れないと2115通り調べることになり計算が終わりません、合法な組み合わせすべてを試すので全探索と主張することは許して…

PEZY-SC/SC2を使った話 【KMCアドベントカレンダー 2日目】

この記事はKMC Advent Calendar 2018 - Adventarの2日目の記事です。

adventar.org

概要

理研に設置されているスーパーコンピューター、菖蒲(Shoubu)・菖蒲SystemBでPEZY Computing社のPEZY-SC/SC2を使わせてもらいました。

性能とか使いやすさについて感想とか知見を述べます。

PEZY-SC/SC2について

PEZY Computing社の開発したプロセッサで、メニーコアのMIMD(Multiple Instruction, Multiple Data)プロセッサであることが特徴です。

MIMDなので各コアで全く異なる命令を独立に実行できる点でGPUなどのSIMDプロセッサと異なります。

使えるようになるまで

理研に設置されているスーパーコンピューター、菖蒲(Shoubu)は一般の人でも申請して認められれば無料で使うことができます。 (CUDAまたはOpenCLでの開発経験があったほうが良さそうですが)

Shoubu システム利用募集のご案内 | 理化学研究所情報システム部(所内および所外向け)

私は一昨年12月頃(ちょうどPEZY Computing社の社長(当時)が逮捕された頃です)に、並列オセロ終盤ソルバーの開発の名目で申請しました。

申請して通ったら書類にサインしたのをスキャンして送ったりsshの公開鍵を送ったりといった手続きをすればアカウントが作成されて使えるようになります。

また、申請して使えるようになるのはPEZY-SCを搭載した菖蒲ですが、PEZY-SC2でしか使えない機能(アトミックなど)を使いたい場合はPEZY-SC2を搭載した菖蒲SystemBを使わせてもらえるようです。

使ってみる

基本的な使い方はGPUと同じで、ホストとなるCPU側のコードとデバイスアクセラレータ)側のコードの両方を書きます。

バイス側のプログラムはclangに食わせるので普通のC++14が書けます(標準ライブラリは使えませんが)。 この点は古いOpenCLで書かざるを得ない環境と比べるとだいぶ楽だと思います。

メモリ空間はホストとデバイスで別なので明示的に転送しないといけません(CUDAのUnified Memoryのようなことはいまのところできない)。

今回はサンプルとしてMIMDの恩恵が受けられるN-Queen問題を実装してみます。

問題

N*Nの盤面上にチェスのクイーンを互いに効きがないように置く方法は何通りあるか? チェスのクイーンは上下左右と斜めに動ける(飛車と角行を合わせたような動き)。

詳しくは エイト・クイーン - Wikipedia などを見てください。

この問題は各行についてどこに置けるかを試すバックトラックで解けることが知られています。 また、駒の利きをビット列で管理することで、ビット演算を用いて高速に解けることも知られています。 並列化は数行を全探索し、その結果を各スレッドに分配して計算させることでできます。

実装

#include "pzc_builtin.h"
#include "../solver.hpp"

uint64_t solve(int N, int depth, uint32_t left, uint32_t mid, uint32_t right) {
  if (depth == N) return 1;
  uint64_t sum = 0;
  for (uint32_t pos = (((uint32_t)1 << N) - 1) & ~(left | mid | right);
      pos; pos &= pos-1) {
    uint32_t bit = pos & -pos;
    sum += solve(N, depth+1, (left | bit) << 1, mid | bit, (right | bit) >> 1);
  }
  return sum;
}

void pzc_Solve(const Problem * const probs, uint64_t * const result, size_t count) {
  size_t offset = get_tid() + get_pid() * get_maxtid();
  result[offset] = 0;
  for (size_t index = offset; index < count; index += get_maxpid() * get_maxtid()) {
    const Problem &prob = probs[index];
    result[offset] += solve(prob.N, prob.depth, prob.left, prob.mid, prob.right);
  }
  flush();
}

GPUなどのSIMDアーキテクチャと違って、PEZYプロセッサでは分岐のペナルティーがとても小さいという特徴があるので、このように再帰の中で分岐をするコードを書いても素直かつ高速に動きます。 一方、GPUで高速に動かそうとすると再帰をスタックを用いたループで書くなどコードを変形する必要があり、プログラミングが困難になります。

ホスト側では数行分の全探索をするコードとデバイス側のプログラムを呼び出すコードを書きます。 デバイス側のプログラムを呼び出すのはOpenCLとほぼ同じです。

コード

CPU/GPU版はこちら。

github.com

GPU版は再帰をスタックを用いたループで書き直す最適化をしてあります(solve_nonrec関数)。 詳しい最適化の内容については

www.slideshare.net

を見てください(102ページ目から)。

PEZY-SC/SC2版はこちら。

github.com

結果

マシン 時間(秒)
Core i7-6700K 78.857027
GeForce GTX 1080 OC 8.400348
Tesla V100 PCIe 3.936454
PEZY-SC 4.67289
PEZY-SC2 2.891

分岐ペナルティーの小ささなどからPEZY-SC2が最も高速という結果になりました。

さらなる最適化

GPU/PEZY-SC/SC2版ともに、タスク(計算する盤面)によって計算量が異なることにより、後半は遊んでしまうコアが出てきます。 GPUやPEZY-SC2ではアトミック命令が使えるので、アトミック命令を使うなどして動的にタスクを割り振ることでさらに高速化することができます。 アトミック命令の使い方などはNDAを結ばないと教えてもらえませんが、ソースコードは公開OKなので、私の公開しているオセロソルバーのソースコード

github.com

などから類推することでおおよそはわかるのではないでしょうか。 実際、PEZY-SC2上でアトミック命令を用いて最適化したところ、2.40199秒まで高速化されました。

GitHub - primenumber/PEZY-nqueen at sc2

オセロソルバーの実装

オセロソルバーとはあるオセロの局面が与えられたときにそこから両者が最善を尽くしたときの試合結果(何石差か)を求めるプログラムのことです。 基本的にはαβ法という再帰的なアルゴリズムを用いて解きます。

すでにCPU/GPU版を実装していたので、簡単にPEZY-SC/SC2版も実装できると思っていました。

CPU版: github.com

GPU版: github.com

しかし、N-Queen問題と同じようにオセロソルバーも実装しようとしたのですが、空きマスが多くなると(=再帰の深さが深くなると)エラーが出て正しく動きませんでした。 原因を探ってみると、どうやらスタックオーバーフローしているようです。 実は、PEZY-SCにはローカルメモリが16KBあり、これを1コア内の8スレッドで分割してスタック領域として使います(ちなみに、スタック領域を縮小して余ったローカルメモリをコア内の共有メモリとして使うこともできます)。 そのため、スタック領域は各スレッドにつき2KBしかありません。PEZY-SC2でローカルメモリが20KBになりましたが、それでもスレッドあたりは2.5KBです。 N-Queenより複雑なオセロソルバーではスタック領域が不足してしまっていたのです。 (再帰で書いているのも一因ですが、スタック消費量を気にせずコードを吐くLLVM側も悪い気はします)

そこで、泣く泣くスタック領域を大量に消費する再帰で書くことを諦めて、自分で実装したスタックをグローバルメモリ上に置いてループで解くことにしました。 これで解けるようにはなりましたが、プログラミングが非常に難解になってしまいました。 このあたりのプログラミングの難しさが解消されるともっと良いと思うのですが…(スタック領域は基本グローバルメモリに置かれて、キャッシュで高速化するなど、難しそうではありますが)

何とかプログラミングの困難さを乗り越えて実装した結果が

github.com

です。 根の1段だけFastest-First Heuristic(速さ優先探索)を行っています(複数段に適用しようとするとスタックが足りなかった)。 また、PEZY-SC2向けに最適化を行ったものが

GitHub - primenumber/PEZY-Othello at sc2

にあります。 最適化の内容は

  • 速さ優先探索を複数段に適用
  • ビット演算命令を使用
  • アトミック命令を使用して動的なタスク割り振り
  • 置換表の利用
  • 葉の近くで速さ優先探索を適用していないところで隅を先に読む静的move ordering
  • NegaScout法の適用

になります(まだパラメータ調整をちゃんとしていなくて速度が出ない場面もありますが…)。 これらの改良により、20石空きの局面を219万局読ませて、1ノード(8モジュール)で6時間弱で解けるようになりました(Tesla V100でも測るべきですが、クラウドで8GPUとかぶん回すと結構高いのでやってません、あとGPU版はいくつかの改良を取り入れてないです)。 菖蒲SystemBでは最大1週間のジョブを投げられるので、22~23石空きまでは1ノードで読ませることができそうです(1石空きマスが増えるごとにだいたい3倍ぐらいかかる時間が増える)(ちなみにノード数を増やしてもゲーム木の大きい局面が律速になるのでそんなに高速化しない)。

今後は統計評価関数を利用したmove orderingを実装して更なる高速化を果たして、24石空き程度まで読ませたいと考えています。

感想

PEZY-SC/SC2はメニーコアのMIMDプロセッサという他にはあまりないアクセラレータです。 今回は再帰の中で分岐するようなMIMDに有利なワークロードではPEZY-SC2がTesla V100を超えるパフォーマンスを出せることを示せました。次世代のPEZY-SC3も近いうちに登場するみたいなので楽しみですね。

一方で、スタック領域の制約にはかなり苦しめられました。せっかく再帰的なタスクに有利なので、たくさんスタック領域が使えるとプログラミングがしやすくてよいと思うので、この辺が改善されてほしいです。

最後に

KMCアドベントカレンダー 3日目の記事は id:nojima718 (KMC-ID: nojima) さんの「ぷよぷよのプレイ動画を解析して棋譜を生成する」です。お楽しみに!

adventar.org

ビット演算マニアのための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の高速な実装