8x8 bitboardへの対称な変換の全列挙にグレイコードを応用する

以前の記事(↓)の続きです。
eukaryote.hateblo.jp
符号なし64bit整数を8x8 bitboardとみなして(1)横方向に鏡映 (2)縦方向に鏡映 (3)行列転置 という3つの処理をこの順番でやってもやらなくてもよいとしたとき、結果得られる最大8通りのbitboardのうち整数としての値が最小のものを得たいという問題を考えます。ベースラインとなるナイーブな実装は以下のようになります:

//horizontal_mirror, vertical_mirror, transposeの実装は省略
uint64_t get_unique_naive(const uint64_t b) {
	uint64_t answer = b;
	for (uint32_t i = 1; i <= 7; ++i) {
		uint64_t candidate = b;
		if (i & 1)candidate = horizontal_mirror(candidate);
		if (i & 2)candidate = vertical_mirror(candidate);
		if (i & 4)candidate = transpose(candidate);
		answer = std::min(answer, candidate);
	}
	return answer;
}

3種類の変換処理をやるかどうかの選択を「変数iを二進数として見たときの下位3桁のビットが立っているかどうか」と対応付けて、iをfor文で回して全通り作っています。

実際の使用例としては、Edaxというオセロソフトでは上記のような形で実装されています。
edax-reversi/src/board.c at 54578eefbb3fe5d734e02c0e6f6ffed11d1bad9e · abulmo/edax-reversi · GitHub

グレイコードの考え方を応用した実装

一方でEgaroucidという別の最強オセロソフトでは異なる実装になっています。
Egaroucid/src/engine/book_accuracy.hpp at 03f7c7df97ae01613f2c205758b76c9fbaa7156a · Nyanyan/Egaroucid · GitHub
この実装を模式的に書くと以下のような感じになります。
(厳密には違っていて、リンク先のEgaroucidの実装では行列転置(=対角線を軸とする鏡映)を2回行っています。3種類の変換処理のうち行列転置が一番重いはずなので、行列転置を1回しかやらない以下のような順序付けのほうが少しだけ高速かもしれません)

uint64_t get_unique_graycode(const uint64_t b) {
	uint64_t answer = b;
	const uint64_t code001 = horizontal_mirror(b);
	answer = std::min(answer, code001);
	const uint64_t code011 = vertical_mirror(code001);
	answer = std::min(answer, code011);
	const uint64_t code010 = horizontal_mirror(code011);
	answer = std::min(answer, code010);
	const uint64_t code110 = transpose(code010);
	answer = std::min(answer, code110);
	const uint64_t code111 = horizontal_mirror(code110);
	answer = std::min(answer, code111);
	const uint64_t code101 = vertical_mirror(code111);
	answer = std::min(answer, code101);
	const uint64_t code100 = horizontal_mirror(code101);
	answer = std::min(answer, code100);
	return answer;
}

3種類の変換処理をやるかどうかの選択を3桁の二進数の各ビットの状態に対応付けるという発想はナイーブ実装と同じです。そのうえでこの実装では「直前のbitboardに対してどれか1種類の変換処理だけを行い次のbitboardを作る」を7回行うことで8通りのbitboardを網羅できています。このように「1ビットだけ逆転させる処理によって任意桁の二進数を巡回できること」やその順序自体については"グレイコード"という名前で広く知られており、様々な分野で応用されています。

ちなみに同様の実装は2048(ゲーム)を強化学習するソフトの実装でも見られます。
TDL2048/board.h at 665d470cd8790c88519352362f91f4209e0986a1 · moporgic/TDL2048 · GitHub
2048の場合は、1マスを4bitとみなして16マスの盤面全体を64bitとみなすことになります。その場合には転置操作をpext命令4回でゴリ押したりできて、縦鏡映や横鏡映と比べたときにどれが最も重い処理なのかはCPUアーキテクチャによって色々と話が変わってきます。

グレイコードの考え方を離れ、クリティカルパスを短くする

上記のグレイコードに基づくアルゴリズムは美しいですが、直前の変換が完了しないと次の変換を計算できない、すなわち計算グラフのクリティカルパスが長いのが欠点です。「変換処理を7回やるだけで8通りを網羅できる」「行列転置は1回しか行わない」という要件でクリティカルパスを短くしようとすると、例えば以下のように書くことができます。

uint64_t get_unique_notgraycode(const uint64_t b) {
	uint64_t answer = b;
	const uint64_t code100 = transpose(b);
	answer = std::min(answer, code100);
	const uint64_t code001 = horizontal_mirror(b);
	answer = std::min(answer, code001);
	const uint64_t code010 = vertical_mirror(b);
	answer = std::min(answer, code010);
	const uint64_t code101 = horizontal_mirror(code100);
	answer = std::min(answer, code101);
	const uint64_t code011 = horizontal_mirror(code010);
	answer = std::min(answer, code011);
	const uint64_t code110 = vertical_mirror(code100);
	answer = std::min(answer, code110);
	const uint64_t code111 = vertical_mirror(code101);
	answer = std::min(answer, code111);
	return answer;
}

この実装を含めた上記3つの実装はすべて同じ結果を返しますが、レイテンシ・スループットの面では一番下のコードが(私の環境では)最も高速でした。特にレイテンシに関しては前回の記事で紹介したAVX2版のものよりも高速でした。

オセロの最大分岐数は33

tl;dr

棋譜は f5f6e6f4g7c6g3e7d6f3e3d3b7d7c2g2g1c3b2b3b4f7g5c4c7c8e2 です。

オセロの局面。33マス空き、白の手番、33マスすべてに着手可能。WZebraで画像作成

経緯(2020)

オセロには相手の石を1つ以上裏返せる場所にしか石を置けないというルールがあります。
わたしはこのブログで2020年に、オセロの最大分岐数(石を置ける場所の数が最大になるような局面)を求めようとした記事をいくつか書いていました。
成果として、初期局面から合法的に到達可能かを考えなければ34箇所が最大であり、35箇所に石を置ける局面は存在しないことを証明しました。また実際に34箇所に置ける局面を構成できました。
オセロの最大分岐数は34以下 - ubiquitinブログ
初期局面から合法的に到達可能な局面に限定する場合に関しては、34箇所が上界となることがこの結果から言えます。

そこで、合法手生成と盤面更新をビット演算と簡単な四則演算で書き下し、初期局面から合法的に到達可能な最大分岐数とその局面をz3に求めさせるコードを書きました。
到達可能局面をSATソルバーに調べさせる - ubiquitinブログ
当時わたしはSMTソルバのz3 v4.6.0 を使っていましたが、そのPython APIでは、与えた制約と充足可能性において等価なCNF(連言標準形)の式を出力する(というか、内部表現を雑に全部書き出させる)ことができました。そこで、その書き出された式をDIMACS CNF形式に変換して、当時のSAT competition 2020でランキング上位だったCaDiCaLというSATソルバに与えてみました。しかし少なくとも数日では解が求まりませんでした。当時使っていたCPUがi7-4770で、2020年時点でも相当古かったせいかもしれませんが、とにかく当時はそこで一旦諦めました。

経緯(2023)

先日それを思い出して、2020年当時に作ったCNF形式のファイルを最新のSAT competition 2022でランキング上位だったkissatのsc2022-bulkyというrelease versionに投げてみたところ、33箇所がSATISFIABLEで(!)34箇所がUNSATISFIABLEだと証明できました。CPUはRyzen 5950Xで、30時間くらいかかりました。
Release SAT Competition 2022 Bulky Release · arminbiere/kissat · GitHub
しかし、その出力(すなわち、制約を満足するような真偽値の割り付け)から、オセロで33箇所に置ける局面への具体的な再構成手順が存在しませんでした。というのも、わたしは2020年時点では33箇所もUNSATなのではと予想していたので、z3の内部表現を標準出力に全部書き出させるという方法を取ったときも、その内部変数に対して真偽値を強制的に渡す方法は調べずにいました。実際、簡単な方法は見当たりませんでした。

そこで、当時z3で書いたコードをC++で書き直しました。CNF形式のファイルだけでなく、論理式上での変数番号とそのオセロソフトにおける意味との対応関係も別ファイルとして書き出すような実装をしました。kissat sc2022-bulkyに投げたところ、8時間くらいで計算できました。出力をPythonで適当に処理すると、結果として冒頭の棋譜が得られ、WZebraで画像を作ることができました。

余談: WZebraの機能で、棋譜の上で何ターン目に置かれたかの番号を石の上に書き出すことができる(Export game as PNG 等)のですが、この棋譜を書き出そうとしたところ石の色が勝手に変わるバグを踏みました。

コードなど

github.com

SVE2のBGRP命令

news.mynavi.jp

この記事で、Scalable Vector Extension 2(SVE2)について語られています。出典は以下のページの右上の"download presentation"ボタンから得られます。

connect.linaro.org

このプレゼンの26ページ目(=一番上の記事のphoto 02)によると、SVE2にはBDEP, BEXT, BGRPという3つのbit shuffle命令が追加されるそうです。(genomics向け?)BDEPとBEXTはBMI2のpdep,pextと等価なんですが、BGRPは一歩進んでいます。

BGRPの仕様は以下のページに書かれています。

https://developer.arm.com/documentation/ddi0602/2021-09/SVE-Instructions/BGRP--Group-bits-to-right-or-left-as-selected-by-bitmask-?lang=en

C++でBMI2で書くとこんな感じになります。第一引数xのビットのうち、第二引数maskで1が立っている場所のビットを下位側に寄せて、残りのビットを上位側に寄せるというものです。

uint64_t BGRP(uint64_t x, uint64_t mask) {
	return pext(x, mask)
		| (pext(x, ~mask) << (_mm_popcnt_u64(mask)));
}

このブログの以前の記事でハッカーのたのしみ7章について扱ったやつがあるのですが、そこで出てきたsag関数に非常に近いです。(完全に同じではありません)

eukaryote.hateblo.jp

(2^n) bit整数のビットを任意の順序にシャッフルする処理は、sag関数(≒BGRP命令)たかだかn回でできます。ちょうどn回でやる方法がハッカーのたのしみ7章に書かれていて、実は場合によってはn回未満でできるということが上のブログ記事で説明されています。

8x8 bitboardへの対称な変換を全部やるのにAVX2を使う

符号なし64bit整数を8x8 bitboardとみなして、(1)横方向に鏡映 (2)縦方向に鏡映 (3)行列転置 という3つの処理をこの順番でやってもやらなくてもよいとしたとき、結果得られる最大8通りのbitboardのうち、整数としての値が最小のものを得たいという問題を考えます。

C++で基本命令だけで書くと以下のようになります。get_unique_naive関数が上記の問題に答える関数です。

uint64_t vertical_mirror(uint64_t b) {
	b = ((b >> 8) & 0x00FF00FF00FF00FFULL) | ((b << 8) & 0xFF00FF00FF00FF00ULL);
	b = ((b >> 16) & 0x0000FFFF0000FFFFULL) | ((b << 16) & 0xFFFF0000FFFF0000ULL);
	b = ((b >> 32) & 0x00000000FFFFFFFFULL) | ((b << 32) & 0xFFFFFFFF00000000ULL);
	return b;
}
uint64_t horizontal_mirror(uint64_t b) {
	b = ((b >> 1) & 0x5555555555555555ULL) | ((b << 1) & 0xAAAAAAAAAAAAAAAAULL);
	b = ((b >> 2) & 0x3333333333333333ULL) | ((b << 2) & 0xCCCCCCCCCCCCCCCCULL);
	b = ((b >> 4) & 0x0F0F0F0F0F0F0F0FULL) | ((b << 4) & 0xF0F0F0F0F0F0F0F0ULL);
	return b;
}
uint64_t transpose(uint64_t b) {
	uint64_t t;
	t = (b ^ (b >> 7)) & 0x00AA00AA00AA00AAULL;
	b = b ^ t ^ (t << 7);
	t = (b ^ (b >> 14)) & 0x0000CCCC0000CCCCULL;
	b = b ^ t ^ (t << 14);
	t = (b ^ (b >> 28)) & 0x00000000F0F0F0F0ULL;
	b = b ^ t ^ (t << 28);
	return b;
}
uint64_t symmetry_naive(const uint32_t s, uint64_t b) {
	uint64_t answer = b;
	if (s & 1)answer = horizontal_mirror(answer);
	if (s & 2)answer = vertical_mirror(answer);
	if (s & 4)answer = transpose(answer);
	return answer;
}
uint64_t get_unique_naive(const uint64_t b) {
	uint64_t answer = b;
	for (uint32_t i = 1; i <= 7; ++i) {
		const uint64_t new_code = symmetry_naive(i, b);
		if (answer > new_code) {
			answer = new_code;
		}
	}
	return answer;
}

これをAVX2で書くことを考えます。まず、vertical_mirror関数とhorizontal_mirror関数は、命令列が完全に同じでオペランドの定数が異なるだけなので、SIMDならば同時に計算できます。例えば以下の関数のように書けます。これの返り値は、(64bit整数4要素の配列とみなすと)[1]が縦方向への鏡映、[2]が横方向への鏡映、[0]と[3]が引数そのままとなります。

__m256i func1(const uint64_t b) {
	constexpr uint64_t FULL = 0xFFFF'FFFF'FFFF'FFFFULL;

	const __m256i bb0 = _mm256_set1_epi64x(int64_t(b));

	const __m256i tt1lo = _mm256_and_si256(_mm256_srlv_epi64(bb0, _mm256_set_epi64x(0, 8, 1, 0)), _mm256_set_epi64x(FULL, 0x00FF00FF00FF00FFLL, 0x5555555555555555LL, FULL));
	const __m256i tt1hi = _mm256_and_si256(_mm256_sllv_epi64(bb0, _mm256_set_epi64x(0, 8, 1, 0)), _mm256_set_epi64x(FULL, 0xFF00FF00FF00FF00LL, 0xAAAAAAAAAAAAAAAALL, FULL));
	const __m256i tt1 = _mm256_or_si256(tt1lo, tt1hi);

	const __m256i tt2lo = _mm256_and_si256(_mm256_srlv_epi64(tt1, _mm256_set_epi64x(0, 16, 2, 0)), _mm256_set_epi64x(FULL, 0x0000FFFF0000FFFFLL, 0x3333333333333333LL, FULL));
	const __m256i tt2hi = _mm256_and_si256(_mm256_sllv_epi64(tt1, _mm256_set_epi64x(0, 16, 2, 0)), _mm256_set_epi64x(FULL, 0xFFFF0000FFFF0000LL, 0xCCCCCCCCCCCCCCCCLL, FULL));
	const __m256i tt2 = _mm256_or_si256(tt2lo, tt2hi);

	const __m256i tt3lo = _mm256_and_si256(_mm256_srlv_epi64(tt2, _mm256_set_epi64x(0, 32, 4, 0)), _mm256_set_epi64x(FULL, 0x00000000FFFFFFFFLL, 0x0F0F0F0F0F0F0F0FLL, FULL));
	const __m256i tt3hi = _mm256_and_si256(_mm256_sllv_epi64(tt2, _mm256_set_epi64x(0, 32, 4, 0)), _mm256_set_epi64x(FULL, 0xFFFFFFFF00000000LL, 0xF0F0F0F0F0F0F0F0LL, FULL));
	const __m256i tt3 = _mm256_or_si256(tt3lo, tt3hi);
	return tt3;
}

ところで、縦方向と横方向の両方への鏡映変換を行うのは、ビットの配置を逆転させることと等価です。64bit整数のビットの配置を逆転させる処理をAVX2で以下のように書けます。

__m256i func2(const uint64_t b) {
	constexpr auto f = [](const uint8_t i) {
		return uint8_t(((i & 1) << 3) + ((i & 2) << 1) + ((i & 4) >> 1) + ((i & 8) >> 3));
	};

	const __m128i r1 = _mm_set_epi8(f(15), f(14), f(13), f(12), f(11), f(10), f(9), f(8), f(7), f(6), f(5), f(4), f(3), f(2), f(1), f(0));

	const __m128i a1 = _mm_set_epi64x((b >> 4) & 0x0F0F'0F0F'0F0F'0F0FULL, b & 0x0F0F'0F0F'0F0F'0F0FULL);
	const __m128i a2 = _mm_shuffle_epi8(r1, a1);
	const __m128i a3 = _mm_shuffle_epi8(a2, _mm_set_epi8(8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 4, 5, 6, 7));
	const __m128i a4 = _mm_shuffle_epi32(a3, 0b00001110);
	const __m128i a5 = _mm_add_epi32(a4, _mm_slli_epi64(a3, 4));
	return _mm256_zextsi128_si256(a5);
}

__m128i型を8bit整数16個とみなして、そこに引数を4bitずつ刻んで入れます。すると、pshufb命令で「4bitを添字とする8bitの配列への表引き」を全部一気に行えますから、各4bitについてビットの配置を逆転させられます。あとは、8bit単位で配置を逆転させて、刻んだのをもとに戻せばOKです。返り値の最下位64bitに欲しい値が入っていることになります。

このfunc1とfunc2は依存関係なく実行できるのでレイテンシを隠蔽できます。返り値2つをblendすることで、「(縦・横)方向への鏡映を(やる・やらない)」の4通りの結果を__m256i型変数1個に格納できます。

あとは、最初に書いた行列転置の関数を単純にSIMDで書き下して入力すれば、所望の8通りの結果を__m256型変数2個に格納できます。最小値を得る処理を含めると、最終的な関数は以下のようになります。

__m256i get_unique_avx2(const uint64_t b) {
	constexpr uint64_t FULL = 0xFFFF'FFFF'FFFF'FFFFULL;

	const __m256i b1 = func1(b);
	const __m256i b2 = func2(b);

	const __m256i tt = _mm256_blend_epi32(b1, b2, 0b00000011);

	const __m256i x1 = _mm256_and_si256(_mm256_xor_si256(tt, _mm256_srli_epi64(tt, 7)), _mm256_set1_epi64x(0x00AA00AA00AA00AALL));
	const __m256i y1 = _mm256_xor_si256(tt, _mm256_xor_si256(x1, _mm256_slli_epi64(x1, 7)));
	const __m256i x2 = _mm256_and_si256(_mm256_xor_si256(y1, _mm256_srli_epi64(y1, 14)), _mm256_set1_epi64x(0x0000CCCC0000CCCCLL));
	const __m256i y2 = _mm256_xor_si256(y1, _mm256_xor_si256(x2, _mm256_slli_epi64(x2, 14)));
	const __m256i x3 = _mm256_and_si256(_mm256_xor_si256(y2, _mm256_srli_epi64(y2, 28)), _mm256_set1_epi64x(0x00000000F0F0F0F0LL));
	const __m256i zz = _mm256_xor_si256(y2, _mm256_xor_si256(x3, _mm256_slli_epi64(x3, 28)));

	const __m256i a1 = _mm256_sub_epi64(zz, tt);
	const __m256i a2 = _mm256_srai_epi32(a1, 32);
	const __m256i a3 = _mm256_shuffle_epi32(a2, 0b11110101);
	const __m256i a4 = _mm256_blendv_epi8(tt, zz, a3);

	alignas(32) uint64_t result[4] = {};
	_mm256_storeu_si256((__m256i*)result, a4);

	result[0] = std::min(result[0], result[1]);
	result[0] = std::min(result[0], result[2]);
	result[0] = std::min(result[0], result[3]);
	return result[0];
}

単体テストベンチマークテストのコードは以下のページにあげてあります。ベンチマーク結果は省略しますが多少速くなりました。

github.com

magic bitboardのmagic numberをSMTソルバで求めようとした

pext命令を軸としてmagic bitboardについて振り返ってみます。

計算速度が重要なソフトで、そのなかである固定されたマスク値でpext命令を頻繁に計算する必要があるがpext命令を使いたくないとします。具体的な問題設定としては、

  • マスク値で立っているビットの数が少ない(16個以下とか)
  • pext命令の返り値を添字として、事前計算された配列にアクセスしたい
  • しかし、肝心のpext命令が使えない(ないしは使いものにならない)

という感じです。

このとき第一に思いつく案はpext命令と等価な関数を基本RISC命令とかで実装することです。ハッカーのたのしみ7章7.4ではcompress関数(=pext命令)の実装方法が議論されていました。わたしも4bitずつ表引きで頑張る方法とかを考えて試してみたことがあります。以下の記事とGitHubにあげてあります。

pdep、pdep抜きで - ubiquitinブログ
pdep-senpai/pext at master · eukaryo/pdep-senpai · GitHub

別の方法として、都合の良いマジックナンバーを見つければ、引数をビットマスクしてからマジックナンバーをかけたときの上位ビットがpext命令の結果と(全)単射になるというのがあります。このことを詳しく説明するために、以下の関数を考えます。

bool is_good_magic_number(
	const uint64_t mask,
	const uint64_t magic_number,
	const uint64_t compromise) {

	std::set<uint64_t> s;
	const uint64_t p = _mm_popcnt_u64(mask);
	assert(p + compromise <= 64);

	const uint64_t siz = 1ULL << p;

	for (uint64_t i = 0; i < siz; ++i) {
		const uint64_t pattern = _pdep_u64(i, mask);
		const uint64_t x = pattern * magic_number;
		const uint64_t result = x >> (64 - p - compromise);
		if (s.count(result) != 0)return false;
		s.insert(result);
	}
	return true;
}

まず、mask値で立っているビットの数がpのとき、任意の整数をmask値でビットマスクした結果の値は2^p通りのどれかになります。上の関数内ではpatternがその全通りを表しています。その各patternにmagic_numberをかけて、論理右シフトで上位ビットを取り出します。結果、「任意の相異なる2つのpatternに関して、取り出される値が必ず異なる」のならば、「patternにmagic_numberをかけて上位ビットを取り出す」操作をpext関数の代用として後段の添字アクセスのために使えるといえます。そのようなmagic_numberであるとき、かつそのときに限り上記の関数はtrueを返します。
第3引数のcompromiseは、うまいマジックナンバーが見つからなかった場合に妥協する度合いを表す値です。妥協しないならcompromise=0で、そのときは取り出された上位ビットがpext命令の結果と全単射の関係になります。compromiseが正の整数の場合、取り出された上位ビットはpext命令の結果と単射の関係になります。後段の配列アクセスにおいて、配列長を(2^compromise)倍の大きさにしなければいけないペナルティが発生しますが、うまいマジックナンバーは容易に見つかるでしょう。

maskが与えられたもとで都合の良いマジックナンバーを見つける既知の方法は、ランダムに生成した数がマジックナンバーとして成立するかを調べまくることです。最初はcompromise=0で調べまくって、見つからなければcompromiseの値を増やして調べなおすわけです。

実例ですが、最強クラスの将棋ソフトのひとつである「やねうら王」は、走り駒のbitboardの処理に関してpextベースの手法とmagic bitboardの両方に対応しています。そのmagic bitboardを使う場合において、飛車が2九または6九の地点にいる場合のmagic numberに関して、compromise=1相当の妥協をしています。(cf: 下記のGitHubのコードの160行目と168行目)

YaneuraOu/bitboard.cpp at f94720b9b72aaa992b02e45914590c63b3d114b2 · yaneurao/YaneuraOu · GitHub

SMTソルバ

ランダムに調べるのではなくSMTソルバに調べさせることもできるのではないかと思い、試してみました。satならばmagic numberが見つかりますし、unsatならば妥協するほかないことを証明できたことになり、すっきりします。コードは以下にあげてあります。

pdep-senpai/magic_bitboard_smt at master · eukaryo/pdep-senpai · GitHub

まずz3で試してみました:

def solve(mask, compromise):

    shiftlen = 64 - popcount(mask) - compromise

    bb = z3.BitVec("bb", 64)
    masked_bb = [z3.BitVec(f"masked_bb_{i}", 64) for i in range(2 ** popcount(mask))]
    magic_number = z3.BitVec("magic_number", 64)
    imul_tmp = [z3.BitVec(f"imul_tmp_{i}", 64) for i in range(2 ** popcount(mask))]
    result_index_short = [z3.BitVec(f"result_index_{i}", popcount(mask) + compromise) for i in range(2 ** popcount(mask))]
    
    s = z3.Solver()
    for i in range(2 ** popcount(mask)):
        s.add(masked_bb[i] == pdep(i, mask))
        s.add(imul_tmp[i] == masked_bb[i] * magic_number)
        s.add(result_index_short[i] == z3.Extract(63, 63 - popcount(mask) - compromise + 1, imul_tmp[i]))
    s.add(z3.Distinct(result_index_short))
    # for i in range(2 ** popcount(mask)):
    #     for j in range(i + 1, 2 ** popcount(mask)):
    #         s.add(result_index_short[i] != result_index_short[j])
    s.add(magic_number >= 1)
    time_start = time.time()
    result = s.check()
    time_end = time.time()

    print(f"mask = {hex64(mask)}")
    print(f"compromise = {compromise}")
    print(f"result = {result}")
    print(f"elapsed time = {int(time_end - time_start)} second")

    if result == z3.unsat:
        return False

    print(f"magic_number = {hex64(s.model()[magic_number].as_long())}")

    return True

maskの立っているビットの数が少ない(9以下とか)ならばうまくいきますが、ビット数13で試したところ、8192要素に対するz3.Distinctがメモリを大量に消費してしまいました。AWSのメモリ384GBのインスタンスで走らせてみたところ、数分でメモリ消費量が50GBを超えて、少ししてsegmentation faultしました。(巨大な配列をなめるときにsize_t型ではなくint型を不用意に使ったときにありがちなことのような気がしますが…)おそらく、z3.Distinctは内部的には(直後にコメントアウトしてあるコードと同様に)愚直に2重にfor文を回しているのではという気がしました。

そこで、集合論理を扱えるcvc5というSMTソルバでも試してみました。

def solve(mask):

    slv = pycvc5.Solver()

    slv.setLogic("QF_ALL")
    slv.setOption("produce-models", "true")
    slv.setOption("output-language", "smt2")


    bitvector64 = slv.mkBitVectorSort(64)
    bitvector_ext = slv.mkOp(pycvc5.kinds.BVExtract, 63, 63 - popcount(mask) + 1)
    bitvector_short = slv.mkBitVectorSort(popcount(mask))
    set_ = slv.mkSetSort(slv.mkBitVectorSort(popcount(mask)))

    shift32 = slv.mkBitVector(64, 32)

    # masked_bb = [slv.mkBitVector(64, pdep(i, mask)) for i in range(2 ** popcount(mask))]
    masked_bb_u = [slv.mkTerm(pycvc5.kinds.BVShl, slv.mkBitVector(64, pdep(i, mask) // (2 ** 32)), shift32) for i in range(2 ** popcount(mask))]
    masked_bb_l = [slv.mkBitVector(64, pdep(i, mask) % (2 ** 32)) for i in range(2 ** popcount(mask))]
    masked_bb = [slv.mkTerm(pycvc5.kinds.BVAdd, masked_bb_u[i], masked_bb_l[i])  for i in range(2 ** popcount(mask))]

    magic_number = slv.mkConst(bitvector64, "magic_number")
    imul_tmp = [slv.mkTerm(pycvc5.kinds.BVMult, masked_bb[i], magic_number) for i in range(2 ** popcount(mask))]
    result_index_short = [slv.mkTerm(pycvc5.kinds.Singleton, slv.mkTerm(bitvector_ext, imul_tmp[i])) for i in range(2 ** popcount(mask))]

    target = [slv.mkTerm(pycvc5.kinds.Singleton, slv.mkBitVector(popcount(mask), i)) for i in range(2 ** popcount(mask))]

    union1 = [slv.mkEmptySet(set_)]
    union2 = [slv.mkEmptySet(set_)]
    for i in range(2 ** popcount(mask)):
        union1.append(slv.mkTerm(pycvc5.kinds.Union, union1[i - 1], result_index_short[i]))
        union2.append(slv.mkTerm(pycvc5.kinds.Union, union2[i - 1], target[i]))

    magic = slv.mkTerm(pycvc5.kinds.Equal, union1[-1], union2[-1])

    print("solve start")
    # print(f"{str(magic)}")

    result = slv.checkSatAssuming(magic)

    print(f"cvc5 reports: magic is {result}")

    if result:
        print(f"For instance, {slv.getValue(magic_number)} is a magic_number.")

この実装だと、マスクのビット数13でもメモリ消費量は5GBくらいのまま動きませんでしたが、そもそも解が求まるまでにz3版よりもかなり長い時間がかかります。数時間待っても求まらないので諦めました。集合論理が使えるからといって早く求まるわけではないのか、それとも私のこの書き方に改善の余地があるのか、cvc5が弱くて他のソルバならもっと速いのか、色々気になります。

pext 6回で8x8 bitboardを転置できる

前回の記事で、符号なし64bit整数を8x8 bitboardとみなして転置する方法として、基本RISC命令のみでやる方法と、AVX2のmovemask命令などを使う方法があるという話について書きました。今回はそれの続きで、ハッカーのたのしみ第7章7.5を改良した方法も加えてベンチマークを取ってみました。

前回の記事:

eukaryote.hateblo.jp

ハッカーのたのしみ第7章7.5の改良についての記事:

eukaryote.hateblo.jp

提案手法

ハッカーのたのしみ7章7.5では、「ワードの一般的順列演算」と筆者らが呼んでいる演算の方法が紹介されていました。その詳細と改良点については上記の記事に書いてあります。それを使うと行列転置はsag関数3回で書けます。平たく書くと以下のようになります。

uint64_t transpose_bitboard_pext(uint64_t x) {
	//引数が8x8 bitboardだとして、転置して返す。

	constexpr uint64_t p = 0xAAAA'AAAA'AAAA'AAAAULL;
	x = (_pext_u64(x, p) << 32) | _pext_u64(x, ~p);
	x = (_pext_u64(x, p) << 32) | _pext_u64(x, ~p);
	x = (_pext_u64(x, p) << 32) | _pext_u64(x, ~p);

	return x;
}

実験結果

ベンチマークでは、xorshift64で生成した乱数を転置させる処理を2^30回行いました。このとき、転置させた結果の値を次の乱数生成への入力にする・しないの両方を試してみました。「しない」場合でも全ての転置結果をxorした値を最後に出力します。CPUは Xeon Gold 6136 (2017年発売、SkyLake世代)で、コンパイラgcc-8.2.0でした。

basic AVX2 pext
する (ms) 6119 4582 6368
しない (ms) 3199 1987 1901

ベンチマークテストのコードは以下のGitHubリポジトリにあげてあります。

github.com

考察

pextをL1T0.5とかにしてほしい(わがまま)