この記事はBrainf*ck Advent Calendar 2019 3日目の記事です。

2日目はみみねこさんによる

でした。4日目はmatsu7874さんによる「RustでBrainfuck インタプリタを実装した話を書けるだろうか？」の予定です。

今回はBrainf*ckの素直な派生言語で、かつメタプログラミングが可能な自作言語、BFmetaの紹介をします。

Brainf*ck派生言語の惨状

Brainf*ckは言語仕様が単純で、コンパイラやインタプリタが作りやすいという特徴があります。そのため、巷にはBrainf*ckの仕様をちょっといじった派生言語がいっぱいあります。

特に多いのが8種類の命令を別の文字列に置き換えただけの派生言語で、Ook!やジョジョ言語などがあります。しかし、これらの言語は単なる文字列の置き換えに過ぎないため、ソースコードの見栄えが異なる以上の違いはなく、プログラミングする上ではあまり面白くないです。

ほかにも、命令の追加でプログラミングの難易度を下げた言語もあり、ビット演算命令を追加するなどしたBrainCrashが代表的な例です。しかし、普通に命令を追加しただけでは多少プログラミングが簡単になるだけであまり面白くありません。また、命令数を増やした分Brainf*ckのシンプルさから来る美しさが失われてしまいます。

かくして、ほとんどのBrainf*ck派生言語はプログラミング言語としてはちっともおもしろくないわけで、その惨状にひどくうなされた私は、夜しか眠れない日々が続きました。

なんとかBrainf*ckの美しさを継承しつつ、本質的な差異をもつ言語が作れないか考えた結果、ついに生み出されたのが、Brainf*ckでメタプログラミング *1を可能にする派生言語、BFmetaです。

BFmeta言語仕様

命令数はBrainf*ckと同じ8種類です。（プログラムの停止を表すヌル文字を9種類目と見ることもできるでしょう）しかし、Brainf*ckとBFmetaの最も大きな違いはプログラム空間とデータ空間が統合されていることです。これにより柔軟なメタプログラミングが可能になります。

まず、8bit(1Byte)*2符号なし整数*3の双方向に無限長*4の配列が用意され*5、0番地からプログラムのソースコードが順に格納され、配列のそれ以外の領域は0で初期化される。データポインタとプログラムポインタという二つのポインタがあり、どちらも最初は0番地を指している。

命令の意味は次の通り、オリジナルのBrainf*ckとほとんど変更されていません。プログラムポインタは毎ステップの最後に1進められ、プログラムポインタの指す値が次の9種類の文字のいずれかのASCIIコードに等しい場合、それに対応する処理が行われます。それ以外の場合、何も行わず（単にプログラムポインタを1進めて）そのステップを終了します。

+ データポインタの指す値を1増やす
- データポインタの指す値を1減らす
> データポインタの値を1増やす
< データポインタの値を1減らす
. データポインタの指す値を標準出力に出力する
, 標準入力から1Byte読み、データポインタの指す先に格納する
[ データポインタの指す値が0なら、プログラムポインタをこの[に対応する]の位置まで進める（今のプログラムポインタの位置から正方向に配列を舐めていき、[の数と]の数が等しくなったら停止）
] データポインタの指す値が0以外なら、プログラムポインタをこの]に対応する[の位置まで進める（今のプログラムポインタの位置から負方向に配列を舐めていき、[の数と]の数が等しくなったら停止）
\0（ヌル文字）プログラムを終了する

Brainf*ckの仕様をご存じの方なら、Brainf*ckとほとんど同じ仕様であることがわかると思います。しかしながら、実際のBFmetaはBrainf*ckよりかなり表現力に長けた*6 プログラミング言語になっています。次節ではその表現力の一部をお見せいたします。

BFmetaプログラミング

Hello World

ソースコード中に書いたデータを実行中に読むことができるので、ソースコード中に文字列を埋め込んでそれを表示するだけです。

[>]<<<<<<<<<<<<<[.>]Hello World!

一旦データの末尾まで移動した後、Hello World!の文字数分だけ戻って、そこから1文字ずつ出力しています。ソースコードにヌル文字を許せば*7、次のように更に短く書けます。

[>]>[.>]\0Hello World!

Brainf*ck インタプリタ

標準入力から最初にヌル文字に到達するまでをソースコード、それ以降を標準入力とみなすことにします。

[>],[>,]>

まずBFmetaソースコードの末尾まで移動し、そこからBrainf*ckのソースコードを読み込み、1Byteデータポインタを進めます（これがないとBrainf*ckプログラムの終了時にBrainf*ckプログラムの配列の先頭が0でないときに誤作動してしまう）。最初のソースコードはここまでですが、連続してさっき読み込んだBrainf*ckプログラムが格納してあり、データポインタはここからそれ以降はすべて0で初期化されているため、このままプログラムを実行し続ければ、読み込んだBrainf*ckプログラムをBrainf*ckプログラムとして実行したときの実行結果が得られます。

Brainf*ck自身でBrainf*ckインタプリタを書くと短いものでも https://arxiv.org/html/cs/0311032 これくらいありますが、それより遥かに短く書けました。

BFmetaでBFmeta自身のインタプリタを書くのはもう少し難しいです。（配列が負の方向に伸びないと仮定するなら、Brainf*ckインタプリタとほぼ同じにはなるのですが…）

Quine

空のソースコードをQuineと言い張ればBFmeta以外にもBrainf*ckやRubyでも0ByteのQuineが書けたことになりますが、それでは面白くないので、大抵は正の長さのソースコードであることが求められます。

keisukenakano.hatenablog.com

では404BytesのBrainf*ckのQuineが紹介されています。これでも十分短いのですが、BFmetaではなんと1文字のQuineが存在します。

まあソースコードをデータとして保持しているので、当たり前ではあります。

このソースコードには改行を含んではいけないことに注意してください。そのため、環境によってはうまく表示されないことがあります。

改行を含めたソースコードをQuineしたい場合は次のようにしてください。

[.>]

これでもBrainf*ckに比べてコード長が1/80です。

BFmeta処理系

こんなに表現力が向上したプログラミング言語、書いてみたい、実行してみたいですよね！そう思って予めインタプリタ（及びデバッガ）をRubyで書いたものがこちらになります。

github.com

実行速度があまり速くないので、今度C++かRustで速いインタプリタを書きたいと思っていますが、それはまた別の記事にしたいと思います。そして、その速いインタプリタによって、Brainf*ckに対して本質的な計算量の改善ができる可能性についても、また別の記事にしたいと思います。乞うご期待。

*1:より正確にはリフレクション

*2:2008年になってようやく、IEC 80000-13により1Byte=8bitと定義されました

*3:オーバーフロー時にビットレベルで同じ挙動をするなら、符号付き整数でも良い

*4:初期は正の方向にのみ無限長だったが、メタプログラミングしやすくするためには、負の方向にも進めるようにしたほうが良い

*5:本当に用意する必要はなく、遅延評価をすることで有限のメモリでインタプリタを実現可能

*6:どちらもチューリング完全なので、計算可能なクラスは同一です

*7:POSIXによるテキストファイルの定義を満たさなくなってしまいますが…

2019-05-27

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

ネタ

概要

qiita.com

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

先行研究

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

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:この工夫を入れないと2¹¹⁵通り調べることになり計算が終わりません、合法な組み合わせすべてを試すので全探索と主張することは許して…

2019-04-01

LeapMind株式会社に入社しました

その他

LeapMind株式会社に入社する

2018-12-02

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

コンピューターオセロ GPGPU

この記事は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での開発経験があったほうが良さそうですが）

http://accc.riken.jp/shoubu_info/application/

私は一昨年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関数)。詳しい最適化の内容については

プログラムを高速化する話Ⅱ 〜GPGPU編〜 from 京大マイコンクラブ

www.slideshare.net

を見てください（102ページ目から）。

PEZY-SC/SC2版はこちら。

github.com

結果

(2020/10/02 追記) N=18、CPUで予め展開する行数は最適なものを選択

マシン	時間(秒)
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

2017-12-22

ビット演算マニアのための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になって変わった主なこと

レジスタ幅が最大512bit(ZMMレジスタ)に（それはそう）
論理レジスタ数が16から32に倍増
マスクレジスタ(k0-k7)の追加
各種ビット操作・バイト操作命令の追加

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をインデックスにして、２つの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}

２つの値をつなげてシフトします。多倍長シフトとかに使えますね。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

2017-12-20

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通したかった…

2016-12-26

オセロの着手可能位置生成【ビット演算テクニック 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などで取ってきて表引きをすることにより、ある程度の速度で可能ですが、(3⁸=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;
}

盤面の端を超えて石が繋がっていても無視する
posから見て上、左、右上、左上のマス目をマスクするビット列を計算する
posからみて相手の石が途切れたところをclzで求め、自分の石の位置とandを取ることでひっくり返せるかどうかがわかる。outflankの各要素は高々1ビットが立つ。
-outflank * 2でoutflankより上位のビットが全部1になる。これとマスクを取ることでひっくり返る石の位置がわかる。
後半部分は下位ビットから見ていって相手の石が途切れたところを探します。これは関係ない場所を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つ目で相手の石の位置を表す流儀もある

	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

	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

	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