置換表と評価関数の小型化

edaxの置換表は1エントリ192バイト 24バイト(192bit)で、パラメータh(ただし 10 \leq h \leq 30 )を受け取って、2.25  2.0625 * 2^{h} 個のエントリをメモリに確保します。(2.25 2.0625と中途半端な値なのは、3種類のテーブルを使い分けているためです)デフォルトはh=21です。最小値のh=10にした場合、大きさは単純計算で約443 51KBになります。この大きさでは探索がめちゃくちゃ非効率になるんだろうなと思いこんでいましたが、実際試したところそうでもありませんでした。

以下は実験結果です。

  • i5-7400
  • 1年前くらいにビルドしたedax
  • fforum-60-79完全読み
  • h=10とh=21を同時にではなく片方ずつ走らせた

h=10の場合 ( -n 1 l 60 h 10 -solve fforum-60-79.obf )

 # | depth|score|       time   |  nodes (N)  |   N/s    | principal variation
---+------+-----+--------------+-------------+----------+---------------------
  1|   24   +20        0:02.123      72731214   34258697 c2 B4 a4 B5 b8 B6 b3 
  2|   25   -14        0:03.991     136638662   34236698 G1 h6 H5 h3 B7 a8 D7 
  3|   27   +28        1:06.431    2337281357   35183594 E8 h5 H2 f7 A3 a2 D7 
  4|   27   -02        0:33.084    1400723145   42338385 F2 f1 C8 b8 D2 d1 G7 
  5|   27   +20        2:16.815    5159947303   37714778 B4 h6 H7 f7 B5 a6    
  6|   28   +10        4:14.072    9680960313   38103216 g1 H1 h2 H3 h4 H7 a5 
  7|   28   +30        3:01.408    6584206879   36295019 h3 E1 g4 G5 f1       
  8|   28   +22        4:39.972   10119623233   36145126 h3 E8 a6 F7 b6 A4 c7 
  9|   30   +28       21:18.760   41480608536   32438150 e8 D8 a5 H5 a4 A6 b3 
 10|   30   +00        7:43.130   16761152666   36191032 h3 H5 a2 A3 a5 A1 g2 
 11|   30   -24        1:17.396    2497576443   32270097 e3 F3 e8 G8 g3 G4 h4 
 12|   31   +20        9:23.139   17123691973   30407576 D2 a5 E2 h4 B3       
 13|   31   +24       53:33.858  106508195403   33140293 E1 f1 H5 h3 G1 f7    
 14|   31   -04       20:59.928   45885561042   36419193 G4 h6 E8 h5 G7 h8 G3 
 15|   31   -30     1:36:56.862  183881154781   31611744 F1 d1 C1 d8 B6 e8 C8 
 16|   32   +14     3:22:16.458  452762950808   37306020 d2 H8 d1 C1 c2       
 17|   32   +32     7:58:23.752  922599899127   32142136 a3 A2 b3             
 18|   34   +34     1:17:03.995  148981322525   32219179 b7 A5 a7             
 19|   34   +08     1:08:25.921  139248536500   33914081 f1 G1 a2 F3          
 20|   36   +64        0:06.199     187928918   30316005 d7 H2 h8 G1 g8 G7 d8 
---+------+-----+--------------+-------------+----------+---------------------
fforum-60-79.obf: 2113410690828 nodes in 17:33:27.294 (33436184 nodes/s).

h=21の場合 ( -n 1 l 60 h 21 -solve fforum-60-79.obf )

 # | depth|score|       time   |  nodes (N)  |   N/s    | principal variation
---+------+-----+--------------+-------------+----------+---------------------
  1|   24   +20        0:01.625      48623514   29922162 c2 B4 a4 B5 b8 B6 b3 
  2|   25   -14        0:03.227      96479978   29897731 G1 h6 H5 h3 B7 a8 D7 
  3|   27   +28        0:37.071    1087819867   29344228 E8 h5 H2 f7 A3 a2 D7 
  4|   27   -02        0:25.841     977664575   37833852 F2 f1 C8 b8 D2 h3 B7 
  5|   27   +20        1:13.044    2382018714   32610738 B4 h6 H7 f7 B5 a5 D2 
  6|   28   +10        5:07.360   10236493800   33304574 g1 G7 c8 H1 h2 D8 a5 
  7|   28   +30        1:27.713    2737168216   31205958 h3 E1 g4 G5 f1 G1 f2 
  8|   28   +22        3:05.550    6084624924   32792374 h3 F8 a6 A5 f7 G6 b6 
  9|   30   +28       16:07.081   25880362470   26761318 e8 D8 a5 H5 a4 H6 b3 
 10|   30   +00        2:36.389    4780894606   30570530 h3 H5 a2 A3 a5 A1 g2 
 11|   30   -24        0:46.389    1237608124   26678914 e3 F3 e8 G8 g2 F2 g4 
 12|   31   +20        4:20.382    6822196039   26200721 D2 a5 E2 h4 B3 b6 C8 
 13|   31   +24       32:29.660   52344624795   26848079 E1 f1 H5 h3 G1 f7 H2 
 14|   31   -04        5:51.877   11251037955   31974349 G4 h6 E8 h5 G7 h8 G3 
 15|   31   -30       48:27.632   73904989901   25417587 F1 d1 C1 d8 B6 e8 C8 
 16|   32   +14       59:05.968  120241764530   33909433 d2 D1 h5 H8 h2 B5 d8 
 17|   32   +32     1:58:09.752  179160564075   25270357 a3 A2 b3 B6 e1 H4 e2 
 18|   34   +34       55:48.837   87152722897   26024773 b7 A5 a7 E8 e7 A8 b2 
 19|   34   +08       37:16.241   64330548902   28767270 f1 G1 a2 F3 c1 F5 a7 
 20|   36   +64        0:01.345      29085426   21624852 d7 D8 h8 H4 g7 F2 g8 
---+------+-----+--------------+-------------+----------+---------------------
fforum-60-79.obf: 650787293308 nodes in  6:33:02.984 (27595630 nodes/s).

h=10だとnodes(読み切るまでに探索した局面の数)は2倍~4倍とかになっていますが、N/s(1秒間あたりの探索局面数)が数割増えています。置換表全体がキャッシュに乗りうるので妥当でしょう。

評価関数のパラメータを小さくする方法

edaxの評価関数自体の説明はこれとかが詳しいです。

もしオセロの初心者がEdaxの評価関数を読んだら - Qiita

評価関数のパラメータは226315*61*2要素のint16_t型配列に格納されていて、そのままだと約56MBです。最後の"*2"は黒側と白側のプレイヤーから見た場合ごとの値でして、符号が反転されてるだけみたいな感じなので、省略しようと思えば等価な結果を維持しつつ省略できます。2番目の"*61"は空白マスの数(=ゲームの進行度)ごとに異なるパラメータにしてあるものです。近い進行度のパラメータは平均してしまえばパラメータ数を減らせます(が、元の評価関数と等価ではなくなります)。

パラメータの最小値と最大値はそれぞれ[-2736, 2613]です。

全パラメータ値を2べき数になるように丸めて、かつ絶対値32未満の値は全部ゼロに丸めるとかしても、完全読みまでの探索局面数はほぼ変わりませんでした。以下は実験結果です。edaxを私が適当にいじったコードにfforum40-59を途中まで解かせたものです。

h=10, パラメータに手を加えない場合

 depth|score|       time   |  nodes (N)  |   N/s    | principal variation
------+-----+--------------+-------------+----------+---------------------
   20   +38       00:00.466      12156948   26087871 a2 B1 c1 PS b6 C7 a7
   22     0       00:01.101      31135531   28279319 h4 A3 a2 G6 g5 H5 f8
   22   +06       00:01.399      38898038   27804172 g2 H1 c2 G1 f1 D2
   23   -12       00:04.035     103307640   25602884 g3 H6 h5 C7 g7 H8
   23   -14       00:01.272      30148768   23701861 b8 G5 d2 A3 b7 A8 a7
   24   +06       00:15.954     460326682   28853371 b2 C1 g5 H4 h6 G4 b1
   24   -08       00:03.147      89810965   28538597 b3 C1 b1 A3 b2 H3 a5
   25   +04       00:00.998      28776680   28834348 g2 B8 b7 A2 a5 B2 g3 
   25   +28       00:08.715     212344400   24365393 f6 G5 g3 F2 h4 H5 h3
   26   +16       00:18.538     494738501   26687803 e1 H4 g6 H6 h5 G4 b2

h=10, 最近傍2べき数へ丸めた場合

 depth|score|       time   |  nodes (N)  |   N/s    | principal variation
------+-----+--------------+-------------+----------+---------------------
   20   +38       00:00.427      11913534   27900548 a2 B1 c1 PS b6 C7 a7
   22     0       00:01.164      34469110   29612637 h4 A3 a2 G6 g5 H5 f8
   22   +06       00:01.320      39203367   29699520 g2 H1 c2 G1 f1 D2
   23   -12       00:04.159     103926798   24988410 g3 H6 h5 C7 g7 H8
   23   -14       00:01.492      32933003   22073058 b8 G5 d2 A3 b7 A8 a7
   24   +06       00:16.458     469465971   28525092 b2 C1 g5 H6 g4 H3 b1
   24   -08       00:03.420      96862525   28322375 b3 C1 b1 A3 b7 G7 a5
   25   +04       00:01.289      39126801   30354384 g2 B8 b7 A2 a5 B2 g3
   25   +28       00:08.912     217675881   24425031 f6 G5 g3 F2 h4 H5 h3
   26   +16       00:19.661     528606175   26886026 e1 H4 g6 H6 h5 G4 b2 

h=10, 最近傍2べき数へ丸め、絶対値32未満の値は全てゼロにした場合

 depth|score|       time   |  nodes (N)  |   N/s    | principal variation
------+-----+--------------+-------------+----------+---------------------
   20   +38       00:00.385      11739796   30492976 a2 B1 c1 PS b6 C7 a7
   22     0       00:01.034      33581599   32477368 h4 A3 a2 G6 g5 H5 f8
   22   +06       00:01.199      37593608   31354135 g2 H1 c2 G1 f1 D2
   23   -12       00:01.540      40451669   26267317 c7 H4 h5 B8 b7 A8
   23   -14       00:01.816      44839684   24691455 b8 G5 d2 A3 b7 A8 a7
   24   +06       00:15.279     453722125   29695799 b2 C1 g5 H6 g4 H3 b1
   24   -08       00:03.494     106354476   30439174 b3 C1 b1 A3 b2 H3 a5
   25   +04       00:01.054      33844793   32110809 g2 B8 b7 A2 a5 B2 g3
   25   +28       00:08.460     226893910   26819611 f6 G5 g3 F2 h4 H5 h3
   26   +16       00:15.944     455651690   28578254 e1 H4 g6 H6 h5 G4 b2

h=10, 最近傍2べき数へ丸め、絶対値32未満の値は全てゼロにして、空きマス数4つごとに区切って平均化した場合

 depth|score|       time   |  nodes (N)  |   N/s    | principal variation
------+-----+--------------+-------------+----------+---------------------
   20   +38       00:00.413      12939393   31330249 a2 B1 c1 PS b6 C7 a7
   22     0       00:01.004      32424064   32294884 h4 A3 a2 G6 g5 G7 f8
   22   +06       00:01.209      38093236   31508052 g2 H1 c2 G1 f1 D2
   23   -12       00:03.830     104595798   27309607 g3 H6 h5 C7 g7 H8
   23   -14       00:01.355      34276650   25296420 b8 G5 d2 A3 b7 A8 a7
   24   +06       00:17.429     506933268   29085619 b2 C1 g5 H4 h6 G4
   24   -08       00:03.711     109704623   29562011 b3 C1 b1 A3 b7 G7 a5
   25   +04       00:01.214      39935299   32895633 g2 B8 b7 A2 a5 B2 g3
   25   +28       00:09.355     243263550   26003586 f6 G5 g3 F2 h4 H5 h3
   26   +16       00:16.862     467434470   27721176 e1 H4 g6 G4 h5 H6 b2

32~2048の2べき数は7通りなので、符号付きでも4bitで表現できます。頑張れば評価関数と置換表を全てキャッシュに収められるかもしれません。

分岐はないほうがいい?

前回はedaxというオセロソフトのboard_score_1関数を分岐不要に書き換える話をしました。cmov命令がそれに使えるのがポイントでした。私はcmov命令の存在を最近知ったのですが、記事を投稿したあとで"cmov"でググったら達人の間では昔から当たり前だったみたいでした。

https://www.slideshare.net/herumi/cmovmaxps

分岐予測器がうまく機能するならば分岐で書くべきだというのはその通りだと思います。分岐予測器の正解率をいつどうやって予測するのかが鍵のような気がします。状況のほうは

  1. 問題の性質によって予測できるケース
    • コーディング時にちょっと考えれば予想できるケース
    • コンパイル時にprofile-guided optimizationすればわかるケース
  2. 入力クエリの性質によって変わるケース
  3. 様々なユーザーが各々のデバイス(異なる分岐予測器を持っている)の上で走らせるケース

に分かれると思います。

以降はしょうもないネタを書いていきます。

board_solve関数について

edaxでは、両者に合法手がない終局局面(空白マスが残っている場合もある)を入力として最終スコアを返す関数が以下のように実装されています。

edax-reversi/endgame.c at 8e5da6d206b750119fe7157e43ea6985b468e8d0 · abulmo/edax-reversi · GitHub

static int board_solve(const Board *board, const int n_empties)
{
	int score = bit_count(board->player) * 2 - SCORE_MAX;	// in case of opponents win
	int diff = score + n_empties;		// = n_discs_p - (64 - n_empties - n_discs_p)

	SEARCH_STATS(++statistics.n_search_solve);

	if (diff >= 0)
		score = diff;
	if (diff > 0)
		score += n_empties;
	return score;
}

オセロのルールでは、空白マスを残した状態で終局した場合、盤上の石が多い方(=勝者)が空白マスを"総取り"します。上記のコードはちゃんと「先手の勝ち・引き分け・後手の勝ち」の3パターンを考慮したコードになっています。

しかし、空白マスの数(=n_empties)が奇数であることが保証されていれば、引き分けがありえないことを利用してちょっとシンプルにできます。例えば以下のように書けます。

static int board_solve_odd(const Board *board, const int n_empties)
{
	assert(n_empties & 1);
	int score = bit_count(board->player) * 2 - SCORE_MAX;	// in case of opponents win
	int diff = score + n_empties;		// = n_discs_p - (64 - n_empties - n_discs_p)

	SEARCH_STATS(++statistics.n_search_solve);

	if (diff > 0)
		score = diff + n_empties;
	return score;
}

edaxでは残り4手以降は個別の関数を用意していて、かつboard_solveは全部インライン化されるはずなので、残り3手の場所でこの関数を使えばちょっと高速化させられるでしょう。

move_evaluateについて

edaxのmove_evaluateという関数のなかに、以下の記述があります。

https://github.com/abulmo/edax-reversi/blob/master/src/move.c

static void move_evaluate(Move *move, Search *search, const HashData *hash_data, const int sort_alpha, const int sort_depth)
{
(前略)
	const int w_low_parity = 1 << 3;
	const int w_mid_parity = 1 << 2;
	const int w_high_parity = 1 << 1;
(中略)
		if (search->n_empties < 12 && search->parity & QUADRANT_ID[move->x]) move->score += w_low_parity; // parity
		else if (search->n_empties < 21 && search->parity & QUADRANT_ID[move->x]) move->score += w_mid_parity; // parity
		else if (search->n_empties < 30 && search->parity & QUADRANT_ID[move->x]) move->score += w_high_parity; // parity
(後略)

QUADRANT_IDは以下のようなテーブルです。

https://github.com/abulmo/edax-reversi/blob/8e5da6d206b750119fe7157e43ea6985b468e8d0/src/search.c

const unsigned char QUADRANT_ID[] = {
		1, 1, 1, 1, 2, 2, 2, 2,
		1, 1, 1, 1, 2, 2, 2, 2,
		1, 1, 1, 1, 2, 2, 2, 2,
		1, 1, 1, 1, 2, 2, 2, 2,
		4, 4, 4, 4, 8, 8, 8, 8,
		4, 4, 4, 4, 8, 8, 8, 8,
		4, 4, 4, 4, 8, 8, 8, 8,
		4, 4, 4, 4, 8, 8, 8, 8,
		0, 0
	};

この処理の意味は、序盤中盤終盤でparity based move sortingの重要度を変えたいということなのでしょう。

(おさらい)parity based move sortingとは、まずオセロの盤面を4つの象限に分けて、各象限に存在する空白マスの個数の偶奇を記録しておき、奇数個空きの象限のマスを優先して探索するというものです。なぜなら、オセロの終盤においては孤立した空白マスに打つほうが一般に有利だからです。「奇数個空き象限」は「偶数個空き象限」より空きマス数の期待値が小さいので、前者を優先するのが良いといえます。

(i7-4700, Visual Studio 2017の)プロファイリングによると、ここのif文の連続が重いようだったので、(分岐予測器が外しまくっているのではと予想して)分岐不要でcmovと表引きに頼るように書き換えたところ、私の環境では少し早くなりました。以下のような感じでやりました。

(前略)
	const int w_low_parity = 1 << 3;
	const int w_mid_parity = 1 << 2;
	const int w_high_parity = 1 << 1;
(中略)
const int32_t parity_table[64] = {
	w_low_parity,w_low_parity,w_low_parity,w_low_parity,w_low_parity,w_low_parity,w_low_parity,w_low_parity,
	w_low_parity,w_low_parity,w_low_parity,w_low_parity,

	w_mid_parity,w_mid_parity,w_mid_parity,w_mid_parity,w_mid_parity,w_mid_parity,w_mid_parity,w_mid_parity,
	w_mid_parity,

	w_high_parity,w_high_parity,w_high_parity,w_high_parity,w_high_parity,w_high_parity,w_high_parity,w_high_parity,
	w_high_parity,

	0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,
	0,0,
};
(中略)
		//if (search->n_empties < 12 && search->parity & QUADRANT_ID[move->x]) move->score += w_low_parity; // parity
		//else if (search->n_empties < 21 && search->parity & QUADRANT_ID[move->x]) move->score += w_mid_parity; // parity
		//else if (search->n_empties < 30 && search->parity & QUADRANT_ID[move->x]) move->score += w_high_parity; // parity
		const bool f1 = (search->parity & QUADRANT_ID[move->x]) != 0;
		const int32_t p1 = parity_table[search->n_empties];
		move->score += f1 ? p1 : 0;

parity based move sortingについて

edaxで、残り4マス空き時点でのparity based move sortingは以下のように記述されています。

https://github.com/abulmo/edax-reversi/blob/8e5da6d206b750119fe7157e43ea6985b468e8d0/src/endgame.c

	// parity based move sorting.
	// The following hole sizes are possible:
	//    4 - 1 3 - 2 2 - 1 1 2 - 1 1 1 1
	// Only the 1 1 2 case needs move sorting.
	if (!(search->parity & QUADRANT_ID[x1])) {
		if (search->parity & QUADRANT_ID[x2]) {
			if (search->parity & QUADRANT_ID[x3]) { // case 1(x2) 1(x3) 2(x1 x4)
				int tmp = x1; x1 = x2; x2 = x3; x3 = tmp;
			} else { // case 1(x2) 1(x4) 2(x1 x3)
				int tmp = x1; x1 = x2; x2 = x4; x4 = x3; x3 = tmp;
			}
		} else if (search->parity & QUADRANT_ID[x3]) { // case 1(x3) 1(x4) 2(x1 x2)
			int tmp = x1; x1 = x3; x3 = tmp; tmp = x2; x2 = x4; x4 = tmp;
		}
	} else {
		if (!(search->parity & QUADRANT_ID[x2])) {
			if (search->parity & QUADRANT_ID[x3]) { // case 1(x1) 1(x3) 2(x2 x4)
				int tmp = x2; x2 = x3; x3 = tmp;
			} else { // case 1(x1) 1(x4) 2(x2 x3)
				int tmp = x2; x2 = x4; x4 = x3; x3 = tmp;
			}
		}
	}

分岐だらけですね。

これは分岐なしで以下のように書けます。

static const uint16_t QUADRANT_ID_2[] = {
	1, 1, 1, 1, 8, 8, 8, 8,
	1, 1, 1, 1, 8, 8, 8, 8,
	1, 1, 1, 1, 8, 8, 8, 8,
	1, 1, 1, 1, 8, 8, 8, 8,
	64, 64, 64, 64, 512, 512, 512, 512,
	64, 64, 64, 64, 512, 512, 512, 512,
	64, 64, 64, 64, 512, 512, 512, 512,
	64, 64, 64, 64, 512, 512, 512, 512
};
static const uint16_t QUADRANT_ID_3[] = {
	6, 6, 6, 6, 48, 48, 48, 48,
	6, 6, 6, 6, 48, 48, 48, 48,
	6, 6, 6, 6, 48, 48, 48, 48,
	6, 6, 6, 6, 48, 48, 48, 48,
	384, 384, 384, 384, 3072, 3072, 3072, 3072,
	384, 384, 384, 384, 3072, 3072, 3072, 3072,
	384, 384, 384, 384, 3072, 3072, 3072, 3072,
	384, 384, 384, 384, 3072, 3072, 3072, 3072
};
static const int8_t answer_order[32][4] = {
	{0,3,1,2}, {0,3,1,2}, {0,1,2,3}, {0,1,2,3},
	{0,1,2,3}, {0,2,1,3}, {1,2,0,3}, {0,1,2,3},
	{0,1,2,3}, {0,1,2,3}, {0,1,2,3}, {0,1,2,3},
	{2,3,0,1}, {0,1,2,3}, {0,1,2,3}, {0,1,2,3},
	{0,2,1,3}, {0,1,2,3}, {1,3,0,2}, {0,1,2,3},
	{0,1,2,3}, {0,1,2,3}, {0,1,2,3}, {0,1,2,3},
	{3,0,1,2}, {0,1,2,3}, {0,1,2,3}, {0,1,2,3},
	{0,1,2,3}, {0,1,2,3}, {0,1,2,3}, {0,1,2,3},
};
inline bool SortWithParity(const int32_t empties[4], int32_t * const x0, int32_t * const x1, int32_t * const x2, int32_t * const x3) {

	const uint64_t p =
		QUADRANT_ID_2[empties[0]] +
		QUADRANT_ID_2[empties[1]] +
		QUADRANT_ID_2[empties[2]] +
		QUADRANT_ID_2[empties[3]];

	const bool b0 = (p & QUADRANT_ID_3[empties[0]]) == 0;//empties[0]の象限に自分しか居なければtrue
	const bool b1 = (p & QUADRANT_ID_3[empties[1]]) == 0;//empties[1]の象限に自分しか居なければtrue
	const bool b2 = (p & QUADRANT_ID_3[empties[2]]) == 0;//empties[2]の象限に自分しか居なければtrue
	const bool b3 = QUADRANT_ID_2[empties[0]] == QUADRANT_ID_2[empties[1]];//empties[0]の象限が、empties[1]の象限と同じならtrue
	const bool b4 = QUADRANT_ID_2[empties[0]] == QUADRANT_ID_2[empties[2]];//empties[0]の象限が、empties[2]の象限と同じならtrue

	const uint64_t i0 = b0 ? 1 : 0;
	const uint64_t i1 = b1 ? 2 : 0;
	const uint64_t i2 = b2 ? 4 : 0;
	const uint64_t i3 = b3 ? 8 : 0;
	const uint64_t i4 = b4 ? 16 : 0;
	const uint64_t parity_index = i0 + i1 + i2 + i3 + i4;

	*x0 = empties[answer_order[parity_index][0]];
	*x1 = empties[answer_order[parity_index][1]];
	*x2 = empties[answer_order[parity_index][2]];
	*x3 = empties[answer_order[parity_index][3]];

	return (((1ULL << 0b00000) + (1ULL << 0b01000) + (1ULL << 0b10000)) & (1ULL << parity_index)) != 0;//全部2ならtrue
}

引数emptiesは空きマス4つの座標で、x0~x4はソート結果とします。b0~b2は、empties[0]~[2]がその象限に1つなのかどうかを示します。b3とb4は、2つの空きマスが同じ象限かどうかを示します。その5bitの情報から、どのマスを優先すべきかを表引きで求めて返します。parity_indexの5bitの情報だけでは、各象限への偏りのパターンを完全に特定はできませんが、安定ソートでなくてよければ必ずソートできます。(例えばparity_index == 0b01100のとき、(x0,x1,x2,x3) が (3,3,1,3) か (2,2,1,1) のどちらかであることまでしかわかりませんが、x2→x3→x0→x1 とソートすればどちらにも対応できます。

さらに、この方法なら(元のedaxと違い)残り3手時点で再びソートする必要をなくすことができます。理由は以下の2点です。

  • 3-1で偏っているケースを正しくソートできる。
  • 2-2で偏っているケースも"ソート"しておけるので、残り4手時点で打った手と同じ象限に属す手を優先させられる。

これがより速いどうかは元のedaxの処理を分岐予測器がうまく捌けているかなどの諸々によるでしょう。

分岐はないほうがいい

edaxというオセロソフトでは、探索処理の中でも最後の1手を打つ部分を計算する関数を個別に用意しています。board_score_1という関数です。

edax-reversi/endgame.c at master · abulmo/edax-reversi · GitHub

inline int board_score_1(const Board *board, const int beta, const int x)
{
	assert(_mm_popcnt_u64(board->player | board->opponent) == 63);
	assert(0 <= x && x < 64);
	assert(((board->player | board->opponent) & (1ULL << x)) == 0);

	register int score, n_flips;

	score = 2 * bit_count(board->opponent) - SCORE_MAX;

	if ((n_flips = count_last_flip(x, board->player)) != 0) {
		score -= n_flips;
	}
	else {
		if (score >= 0) {
			score += 2;
			if (score < beta) { // lazy cut-off
				n_flips = count_last_flip(x, board->opponent);
				score += n_flips;
			}
		}
		else {
			if (score < beta) { // lazy cut-off
				if ((n_flips = count_last_flip(x, board->opponent)) != 0) {
					score += n_flips + 2;
				}
			}
		}
	}

	return score;
}

(頭に3つあるassertは、わかりやすさのため私が書き加えました。)この関数は引数boardで示される盤面が1マス空きであることを仮定しています。また、引数xが空いているマスの座標を示していることを仮定しています。そのうえで、完全読みのスコアを求めて返します。(相手方から見たスコアです。)引数betaはアルファベータ法のbetaです。ゆえに真の完全読みスコアがbeta以上である場合は、「beta以上かつ真のスコア以下」の値が返されることのみが保証されます。

このなかで呼ばれているcount_last_flip関数は、最終局面においてプレイヤーが位置xに石を置いた場合にひっくり返る石の数の2倍を返す関数です。(2倍である理由は単に速さのためです。表引きで値を計算するので、表の値をあらかじめ2倍しておけば手っ取り早いわけです)これはAVX2で高速化した実装が存在していて、以下のように書かれています。

https://github.com/okuhara/edax-reversi-AVX/blob/master/src/count_last_flip_sse.c

int last_flip(int pos, unsigned long long P)
{
	unsigned char	n_flips;
	unsigned int	t;
	const unsigned char *COUNT_FLIP_X = COUNT_FLIP[pos & 7];
	const unsigned char *COUNT_FLIP_Y = COUNT_FLIP[pos >> 3];

	__m256i	MP = _mm256_and_si256(_mm256_broadcastq_epi64(_mm_cvtsi64_si128(P)), mask_dvhd[pos].v4);

	n_flips  = COUNT_FLIP_X[(unsigned char) (P >> (pos & 0x38))];
	t = _mm256_movemask_epi8(_mm256_sub_epi8(_mm256_setzero_si256(), MP));
	n_flips += COUNT_FLIP_Y[(unsigned char) t];
	t >>= 16;

	n_flips += COUNT_FLIP_Y[t >> 8];
	n_flips += COUNT_FLIP_Y[(unsigned char) t];

	return n_flips;
}

(原コードではプリプロセッサマクロでSSE版と書き分けられていますが、わかりやすさのため適宜削除しました)

board_score_1が分岐だらけの理由と分岐をなくす方法

board_score_1関数はゲーム木の末端で呼ばれるだけあって呼び出し回数が非常に多く、プロファイリングすると実行時間が上位に来る関数です。なので高速化できれば嬉しいです。おそらく分岐が多いせいで時間がかかっているのだろうと考えました。この分岐の多さの原因はbeta cutをやっているからです。beta cutをやりたい理由は、count_last_flip関数が重いのでなるべく呼びたくないからです。board_score_1関数は最初にcount_last_flip関数を1回必ず呼びますが、beta cutできれば2回目は呼ばずに済むわけです。

しかしながら、よく見るとAVX2版count_last_flip関数は命令レベルでかなり逐次的なので、自分と相手のビットボードに対するcount_flipを同時並行で計算してもレイテンシはほぼ変わらないと期待できます。つまり以下のような関数を用意すれば、count_last_flip関数を1回呼ぶのとほぼ同じレイテンシで、2回呼ぶのと等価な結果を得られます。

void count_last_flip_double(uint64_t position, uint64_t player, uint64_t opponent, int32_t *p_flipped, int32_t *o_flipped)
{
	const uint8_t * const COUNT_FLIP_X = COUNT_FLIP[position & 7];
	const uint8_t * const COUNT_FLIP_Y = COUNT_FLIP[position >> 3];

	const __m256i PPPP = _mm256_set1_epi64x(player);
	const __m256i OOOO = _mm256_set1_epi64x(opponent);
	const __m256i mask = _mm256_loadu_si256((__m256i*)mask_dvhd_double[position]);

	const __m256i maskPPPP = _mm256_and_si256(PPPP, mask);
	const __m256i maskOOOO = _mm256_and_si256(OOOO, mask);

	*p_flipped = COUNT_FLIP_X[(player >> (position & 0x38)) & 0xFF];
	*o_flipped = COUNT_FLIP_X[(opponent >> (position & 0x38)) & 0xFF];

	const uint32_t tP = _mm256_movemask_epi8(_mm256_sub_epi8(_mm256_setzero_si256(), maskPPPP));
	const uint32_t tO = _mm256_movemask_epi8(_mm256_sub_epi8(_mm256_setzero_si256(), maskOOOO));

	*p_flipped += COUNT_FLIP_Y[tP & 0xFF];
	*p_flipped += COUNT_FLIP_Y[tP >> 24];
	*p_flipped += COUNT_FLIP_Y[(tP >> 16) & 0xFF];

	*o_flipped += COUNT_FLIP_Y[tO & 0xFF];
	*o_flipped += COUNT_FLIP_Y[tO >> 24];
	*o_flipped += COUNT_FLIP_Y[(tO >> 16) & 0xFF];
}

そのうえで、三項演算子はconditional mov命令で分岐なしで処理されうることを踏まえつつ、board_score_1関数を以下のように書き換えられます。

int32_t board_score_1_branchfree(const Board &board, const uint64_t x) {

	const int32_t score = 2 * bit_count(board.opponent) - SCORE_MAX;

	int32_t p_flip, o_flip;
	count_last_flip_double(x, board.player, board.opponent, &p_flip, &o_flip);

	const int32_t x1 = score - p_flip;
	const int32_t x2 = score + 2;
	const int32_t x3 = x2 + o_flip;
	const bool b1 = p_flip != 0;
	const bool b2 = score >= 0;
	const bool b3 = o_flip != 0;
	const bool b4 = b2 | b3;
	const int32_t ax = b4 ? x3 : score;
	return b1 ? x1 : ax;
}

この関数は元のboard_score_1関数と完全に等価ではなくて、beta cutを行わずに常に真の完全読みスコアを返します。そのため、元の関数より高速でありつつ、探索ノード数もちょっと削減されます。

df-pn亜種でオセロの終盤完全読みを試みた

はじめに一般の詰将棋などの問題設定について説明します。
局面をノードとし、指し手を有向エッジとする木(ゲーム木)があるとします。プレイヤーは「攻め」と「受け」の2人とします。葉ノード(詰将棋の場合は、即詰みor王手をかけられない局面のこと)にはbool値が書かれていて、攻めの勝ちならtrue、受けの勝ちならfalseだとします。内部ノードの値は葉から根に向かって以下のように再帰的に求められます:

  • 攻めの手番の場合、子ノードのどれかがtrueならそのノードはtrueである。すべての子ノードがfalseならそのノードはfalseである。
  • 受けの手番の場合、子ノードのどれかがfalseならそのノードはfalseである。すべての子ノードがtrueならそのノードはtrueである。

前者をORノード、後者をANDノードと呼ぶこともあります。ミニマックス法におけるmaxノード/minノードの特殊例と見ることができます。

何手詰めなのか事前にわからないとしますと、反復深化深さ優先探索を行うのが普通です。つまり、最大深さDを設定して深さ優先探索し、根からの距離Dの内部ノードにたどり着いた場合は0.5みたいな値を返すことにして、(true=1,false=0でミニマックス探索して)求まった根ノードの値が0.5だったらDを増やして再探索するわけです。df-pnやその亜種では最大深さの代わりに証明数・反証数というのを閾値として用い、内部で再帰的に反復深化深さ優先探索をします。

文献

アルゴリズムの詳細な説明をここに書くのはだるいので、自分への備忘録も兼ねて文献を紹介するにとどめておきます。

ゲーム計算メカニズム (コンピュータ数学シリーズ 7)

ゲーム計算メカニズム (コンピュータ数学シリーズ 7)

ゲーム計算メカニズム (コンピュータ数学シリーズ 7)

良い本です。下記のわたしの実装もこれに載っていた疑似コードをもとにしています。

GAME-TREE SEARCH USING PROOF NUMBERS THE FIRST TWENTY YEARS

https://dke.maastrichtuniversity.nl/m.winands/documents/ICGA2012PNS.pdf
本当に良いレビュー論文です。2012年という公開時期も良い感じです。

関連トピック

その1: positionとstateの違い

将棋には千日手というルールがあります。同一局面が4回出現したとき、片方が連続して王手をかけ続けているなら王手をかけている側が負けで、さもなくば引き分けになるやつです。このルール設定は「どの局面がどの順番で何回出たか」をプレイヤーが全部覚えていること(perfect recall)を暗に仮定しています。そのため、局面という単語に曖昧さが生じています。

  • 《局面 / position》:駒の配置と手番のこと。スナップショット的なイメージ
  • 棋譜 / history》:開始局面から現在までのpositionの移り変わりの履歴
  • 《局面 / state》:position + history

みたいな混同がされています。「同一『局面』(←position)が4回出現した『とき』(←state)」みたいな。千日手が絡む詰将棋を正しく解くためには、positionが同じでもhistoryが異なるようなstateは本来は区別して扱う必要があります。positionとstateの区別が曖昧なせいでサイクル問題とかGHI(Graph History Interaction)問題とか呼ばれる問題が出てきます。でも区別すると計算が重くなるじゃんとか色々ありますが、オセロに千日手は無いので今回は無視します。

その2: 2重カウント

オセロには千日手はありませんが、positionが同じでhistoryが異なるstateは無数にあります。transposition tableでこれらを同一視することで探索量を著しく削減できます。しかしこのせいで証明数の二重カウント問題という問題が重く発生します。これに対応するため、df-pnの式をいじったWPNS(weak proof number search)という手法が提案されています。上記のレビュー論文で詳しく言及されています。

その3: 振動問題

これも上記のレビュー論文で詳しく言及されていますが、同じくらい有望な手が複数あったとき、「その手の先の展開をちょっと読みに行って戻ってきてもう片方の手に行く」というのを繰り返し、行きつ戻りつのオーバーヘッドが大きくなる問題です。1+\epsilon trickという対策があります。

その4: データセットの偏り

90年代からゼロ年代にかけての詰将棋ソフト界隈では、「人間が作った詰将棋問題の全てを計算機で解きたい」というのをグランドチャレンジと見ていたようで、江戸時代のスゴい問題集が全部解けたとかそういうのが偉業として語られています。しかしそれはある種の偏りを帯びたデータセットです。そもそも長手数の詰将棋問題を人間が手で作れたという点からして、それはゲーム木として見たときに"メチャクチャ細長い"はずです。なぜならそうでないと手動で検算できないからです。(df-pnをはじめとする)あの界隈の凝った手法たちは、そういうメチャクチャ細長いケースにoverfitしている可能性があるという点は心に留めておくべきです。実際の対局中に出現する長手数の詰み手順が効率的に発見できなくてもがっかりしないようにしましょう。

実装

「完全読みの結果がある値を超えるか否か」という問いにtrue/falseで答えさせるようにして、解の二分探索をする感じで実装しました。長いのでGitHubにあげました。とりあえず単一のcppファイルだけでコンパイルできて動くように書いてあります。

https://github.com/eukaryo/algorithm-study/blob/master/DFPN_single_file.cpp

実装されていない重要テクニックは2つで、終盤は普通のNull Window Searchで読むやつと、ハッシュテーブルを使い切ったときに賢くガベコレするやつです。

枝刈り手法について

edaxが使っている枝刈り手法を説明します。

  • stability cutoff: オセロの終盤では「今後絶対にひっくり返らない石」が生じます。これをstable discと呼びます。四隅が最も有名ですが、それ以外にも例えば上下左右斜めの8方向に空白が一切存在しない石は自明にstableです。もちろんその条件に当てはまらない非自明なstable discも存在しえます。 edaxでは「手番側の石のうち自明にstableな石を高速に列挙する関数」を搭載していて、これで手番側の最終スコアの下界を計算できます。その下界がβ値よりも低ければ枝刈りできます。
  • (enhanced) transposition cutoff: 置換表にスコアが書かれていればそれで枝刈りできます。また、任意の指し手に対して、それで1手進めた先でtransposition cutoffが生じるかを先に判定してしまうことも可能で、これはenhanced transposition cutoffと呼ばれます。これはオーバーヘッドが少し大きいので、中盤の限られた場面でのみ試みます。
  • probcut: 浅く読んだ結果が深く読んだ結果の予測値になっていると考えられますが、edaxではその誤差の標準偏差をモデル化してあります。そのモデルによって一定以上の確率でαβカットされることが予想される場合に探索を打ち切ることができます。これをprobcutと呼びます。上2つと違い完全性が保証されないので完全読みの本筋では使えませんが、move orderingのための読みの中では使えます。

edaxには探索ルーチンがいくつかあります。中盤と終盤とかで色々分けてあるのですが、どれも基本的には以下の手順で探索します。

  1. stability cutoffなどの枝刈り手法を試し、枝刈りできたらそこで終了。
  2. 全合法手とそれでひっくり返る石を全列挙する。
  3. 列挙した合法手各々の有望度を計算する。
  4. 合法手を有望度で降順にソートする。
  5. 有望な順にその先を探索する。

2~4番目に改良の余地があります。いくつかの論点に分けて以下で説明します。

まず、edaxはwipeout(=それを打つとすべての石が自分の石になり即勝利する手)を最も有望と評価します。1手詰みで+64点が得られるのですから当然ですね。そして置換表の手が次に有望とみなされます。edaxの指し手評価ではそれらの手は他の任意の手よりも必ず高く評価されます。この設計自体は妥当なのですが、それならば2番目の処理中にwipeoutの手が存在するとわかった時点で残りのmove orderingを全て放棄してβカットしてよくて、そうすると高速化できます。イメージとしては「wipeout cutoff」みたいなのを2番目の処理と同時に計算する感じです。

2番目の処理を終えた時点でwipeoutが見つからなかった場合、3番目の有望度計算をする前に置換表の手の探索を行ってよいです。それでβカットされれば3・4番目の処理を省略できるからです。

wipeoutの手が存在する局面では、置換表に手が書かれているならばそれは必ずwipeoutの手です。ゆえに、置換表の手を探索するのは*2番目の*処理の前とすべきです。まとめると、「置換表の手があればそれを探索する→2番目の処理をやる(置換表の手がなかった場合はwipeoutかどうかのチェックも行う)→3番目に進む」となります。

4番目のソートですが、すべての手を探索する前にβカットが発生しうるので、最初に全部ソートするのではなく「未探索の手のうち有望度最大の手を得る関数」を都度呼ぶべきです。選択ソートを1ステップずつやるイメージです。あるいはヒープソートのように、最初にヒープをO(N)で構築すればsecond best以降の手を得るのはO(logN)で済むみたいな手法を使う手も考えられます。これはエレガントですが、実際にはオセロの合法手は少ないのでヒープソートとかは定数倍遅いと思います。

言い換えると、2~4番目の処理は一つのコルーチンのように書くべきだということです。「有望度最大の手を特定できたらその手をyield returnして、探索結果がβカットでない場合のみ続きを計算する」わけです。ちなみに、これらの点は将棋ソフトでは以前から考慮されていました。将棋では駒を取る手だけを高速に生成するなどの実装が可能なので、タダ取りとか駒得できる手を素早く得ることができます。また、将棋は合法手の数が比較的多い(歩以外の持ち駒があるとすぐ100通りを超えます)ため、ソートのコストも無視できません。

オセロ終盤の探索をちょっと高速化する手法

edaxというオセロソフトのソースコードをガチャガチャいじって遊んでいたのですが、いろいろわかって満足したので忘れる前に書いておきます。

探索アルゴリズムの基本的な説明

最初なので基本的というか、将棋ソフトとかと共通していてedaxでも使われているキーワードを書いておきます。

ミニマックス探索

ゲームの局面をノードとして合法手を有向エッジとする木をゲーム木と呼びます。ゲーム木の葉ノードでは最終スコアが得られます。勝ち負けだけに興味があるなら「先手の勝ちは+1、引き分けは0、後手の勝ちは-1」とか、適当にスコアに換算します。ゼロサムゲームなので、先手はスコアを高くしたいのと同様に後手はスコアを低くしたいと考えるでしょう。

子ノードの最終スコアが全て判明しているような内部ノードにおいて、それが先手の手番ならばスコア最大の子に行くエッジを選びます。後手は最小のエッジを選びます。すると深さ優先探索によって任意のノードの最終スコアが得られます。これをミニマックス探索と呼びます。

アルファベータ法

深さ優先探索において、n手目(n≧2)を読むとき(i.e.内部ノードを通りがけるとき)には1~(n-1)手目の厳密解がわかっているので、それらの最大スコアを超えられないことが判明した時点でn手目の探索を打ち切ることができます。これをアルファベータ法と呼びます。最良計算量はミニマックス探索の平方根です。ただしそれを達成するためには最善手を最初に探索する必要があります。(これは無茶振りでして、最善手が既知ならそもそも探索する必要はありません)逆に最悪計算量となるのはスコアが小さい順に探索した場合で、そのオーダーはミニマックス探索のと同じです。

scout法

アルファベータ法のアルファとベータの値の差を1にして(α+1=β)探索すると、最終スコアがα以下かβ以上のどちらなのかを高速に判定できます。これをNull Window Search(NWS)と呼びます。2手目以降はまずNWSして、n-1手目までの最大スコアを超えるとわかった場合のみ通常の探索をするわけです。これをscout法とかPrincipal Variation Search(PVS)とか呼びます。

aspiration search

アルファとベータの差を適切に小さい値にすると、その範囲外の超悪手とかは高速に枝刈りされます。しかも、最終結果の値が(α、β)の区間内だった場合は厳密解であることが保証されます。この探索手法をaspiration searchと呼びます。区間外だった場合は区間を広げて再探索します。

transposition table

transpositionとは、原義では「異なる棋譜によって同じ局面に到達すること」を意味します。つまりゲーム木が木ではなくDAGになっている状態のことです。その場合、深さ優先探索でとある内部ノードの厳密解が求まったとき、帰りがけにその値をハッシュテーブルに書き込んでおいて、またその局面に到達したときにテーブルの値を引けば再探索を省略できるわけです。このハッシュテーブルをtransposition tableと呼びます。日本語では置換表と呼ばれています。

iterative deepening

置換表にスコアだけでなく最善手も書き込んでおくこともできます。すると、深く探索する前にあらかじめ浅く探索しておいて有望な手を保存し、深く探索するときにその手を優先するというテクニックが可能になります。これはアルファベータ法の「最善手を最初に探索せよ」という無茶振りへの実用的な答えになっています。iterative deepeningとか反復深化深さ優先探索とか呼ばれています。ちなみに、これをやりつつ再探索を省略するやつもやる場合、完全読みの完全性を保証するためには、スコアの根拠となる探索深さも書き込んでおく必要があります。

評価関数

浅く読むということは、ゲームの途中の局面で探索を打ち切るということです。このときは最終スコアが得られないので、代わりとして局面のスナップショットを入力すると完全読み結果の予測値を返す関数が必要になります。これは(静的)評価関数と呼ばれています。edaxでは石の局所的な配置パターンを素性ベクトルとして重み付き線形和を計算しています。

ちょっと面白い関数

search_solve_3は、3箇所が空いている局面を受け取って厳密解を返す関数です。オセロ特有のテクニックが2つ使われています。

https://github.com/abulmo/edax-reversi/blob/master/src/endgame.c

static int search_solve_3(Search *search, const int alpha)
{
	Board *board = search->board;
	Board next[1];
	SquareList *empty = search->empties->next;
	int x1 = empty->x;
	int x2 = (empty = empty->next)->x;
	int x3 = empty->next->x;
	register int score, bestscore;
	const int beta = alpha + 1;

	SEARCH_STATS(++statistics.n_search_solve_3);
	SEARCH_UPDATE_INTERNAL_NODES();

	// parity based move sorting
	if (!(search->parity & QUADRANT_ID[x1])) {
		if (search->parity & QUADRANT_ID[x2]) { // case 1(x2) 2(x1 x3)
			int tmp = x1; x1 = x2; x2 = tmp;
		} else { // case 1(x3) 2(x1 x2)
			int tmp = x1; x1 = x3; x3 = x2; x2 = tmp;
		}
	}

	// best move alphabeta search
	if ((NEIGHBOUR[x1] & board->opponent) && board_next(board, x1, next)) {
		bestscore = -board_solve_2(next, -beta, x2, x3, search);
		if (bestscore >= beta) return bestscore;
	} else bestscore = -SCORE_INF;

	if ((NEIGHBOUR[x2] & board->opponent) && board_next(board, x2, next)) {
		score = -board_solve_2(next, -beta, x1, x3, search);
		if (score >= beta) return score;
		else if (score > bestscore) bestscore = score;
	}

	if ((NEIGHBOUR[x3] & board->opponent) && board_next(board, x3, next)) {
		score = -board_solve_2(next, -beta, x1, x2, search);
		if (score > bestscore) bestscore = score;
	}

	// pass ?
	if (bestscore == -SCORE_INF) {
		// best move alphabeta search
		if ((NEIGHBOUR[x1] & board->player) && board_pass_next(board, x1, next)) {
			bestscore = board_solve_2(next, alpha, x2, x3, search);
			if (bestscore <= alpha) return bestscore;
		} else bestscore = SCORE_INF;

		if ((NEIGHBOUR[x2] & board->player) && board_pass_next(board, x2, next)) {
			score = board_solve_2(next, alpha, x1, x3, search);
			if (score <= alpha) return score;
			else if (score < bestscore) bestscore = score;
		}

		if ((NEIGHBOUR[x3] & board->player) && board_pass_next(board, x3, next)) {
			score = board_solve_2(next, alpha, x1, x2, search);
			if (score < bestscore) bestscore = score;
		}

		// gameover
		if (bestscore == SCORE_INF) bestscore = board_solve(board, 3);
	}

 	assert(SCORE_MIN <= bestscore && bestscore <= SCORE_MAX);
	return bestscore;
}

ひとつはparity based move sortingです。まず前提として、オセロの終盤では離れて孤立した空白マスに打つのは一般に有望な手です。なぜなら、それでひっくり返した石を相手が取り戻すのは多くの場合不可能だからです。そこで、edaxではオセロ盤面を4つの"象限"に分割して、各象限の空きマスの数の偶奇を常に保持しています。偶奇だけ保持することで計算を軽くできます。そして最終盤では、「空きマス数が奇数であるような象限の中の空きマス」を先に探索します。(偶数の象限を後回しにするとも言えます)なぜなら、合計3マス空きのとき偶数の象限には(空きがあるならば)必ず2箇所の空きがあり、密な可能性が高いからです。

ふたつめはNEIGHBOURというビットマスクです。手を探索する前にそれが本当に合法手なのか(=そこに打ったときにひっくり返る石があるのか)確かめる必要がありますが、しかし合法手生成やflipは比較的重いのが問題です。なので、絶対にひっくり返せないケースを高速に検出してスクリーニングできれば有用です。そのためにedaxでは、座標xの8近傍のビットだけが立っているビットボードの配列 NEIGHBOUR を用意していて、NEIGHBOUR[x] と相手の石のビットボードとのandが非ゼロかどうかでスクリーニングしています。空きマスの8近傍に相手の石がなければ、明らかに何もひっくり返せないので非合法手だと言えるわけです。

https://github.com/abulmo/edax-reversi/blob/master/src/bit.c

/** Conversion array: neighbour bits */
const unsigned long long NEIGHBOUR[] = {
	0x0000000000000302ULL, 0x0000000000000705ULL, 0x0000000000000e0aULL, 0x0000000000001c14ULL,
	0x0000000000003828ULL, 0x0000000000007050ULL, 0x000000000000e0a0ULL, 0x000000000000c040ULL,
	0x0000000000030203ULL, 0x0000000000070507ULL, 0x00000000000e0a0eULL, 0x00000000001c141cULL,
	0x0000000000382838ULL, 0x0000000000705070ULL, 0x0000000000e0a0e0ULL, 0x0000000000c040c0ULL,
	0x0000000003020300ULL, 0x0000000007050700ULL, 0x000000000e0a0e00ULL, 0x000000001c141c00ULL,
	0x0000000038283800ULL, 0x0000000070507000ULL, 0x00000000e0a0e000ULL, 0x00000000c040c000ULL,
	0x0000000302030000ULL, 0x0000000705070000ULL, 0x0000000e0a0e0000ULL, 0x0000001c141c0000ULL,
	0x0000003828380000ULL, 0x0000007050700000ULL, 0x000000e0a0e00000ULL, 0x000000c040c00000ULL,
	0x0000030203000000ULL, 0x0000070507000000ULL, 0x00000e0a0e000000ULL, 0x00001c141c000000ULL,
	0x0000382838000000ULL, 0x0000705070000000ULL, 0x0000e0a0e0000000ULL, 0x0000c040c0000000ULL,
	0x0003020300000000ULL, 0x0007050700000000ULL, 0x000e0a0e00000000ULL, 0x001c141c00000000ULL,
	0x0038283800000000ULL, 0x0070507000000000ULL, 0x00e0a0e000000000ULL, 0x00c040c000000000ULL,
	0x0302030000000000ULL, 0x0705070000000000ULL, 0x0e0a0e0000000000ULL, 0x1c141c0000000000ULL,
	0x3828380000000000ULL, 0x7050700000000000ULL, 0xe0a0e00000000000ULL, 0xc040c00000000000ULL,
	0x0203000000000000ULL, 0x0507000000000000ULL, 0x0a0e000000000000ULL, 0x141c000000000000ULL,
	0x2838000000000000ULL, 0x5070000000000000ULL, 0xa0e0000000000000ULL, 0x40c0000000000000ULL,
	0, 0 // <- hack for passing move & nomove
};

しかし、よく考えると8近傍に相手の石があってもなお、明らかに何もひっくり返せないケースは存在します。

f:id:eukaryo:20200426024507p:plain

例えばこの図(黒の手番だとします)で、D2の8近傍には3つの白石がありますが、これはひっくり返せません。NEIGHBOURのスクリーニングが偽陽性になる例です。しかし「外周の石は絶対にひっくり返せない」とは限らなくて、H2に打てばH3をひっくり返せます。じゃあ「外周の空白ますに打てば外周の石をひっくり返せる」というわけでも必ずしもなく、H2はH1とG1をひっくり返せません。

結論を直観的に言うと、空白マスからある8近傍のマスを見たときの位置関係が、「相手の真後ろに壁があるように見える」形になっている場合は絶対にひっくり返せません。そういう位置関係をビットマスクから除外すると以下のようになります。edaxのNEIGHBOURをこれに差し替えるとほんの少しだけ速くなります。

const uint64_t NEIGHBOUR_WITHOUT_KABEDON[] = {
	0x0000000000000302ULL, 0x0000000000000604ULL, 0x0000000000000e0aULL, 0x0000000000001c14ULL,
	0x0000000000003828ULL, 0x0000000000007050ULL, 0x0000000000006020ULL, 0x000000000000c040ULL,
	0x0000000000030200ULL, 0x0000000000060400ULL, 0x00000000000e0a00ULL, 0x00000000001c1400ULL,
	0x0000000000382800ULL, 0x0000000000705000ULL, 0x0000000000602000ULL, 0x0000000000c04000ULL,
	0x0000000003020300ULL, 0x0000000006040600ULL, 0x000000000e0a0e00ULL, 0x000000001c141c00ULL,
	0x0000000038283800ULL, 0x0000000070507000ULL, 0x0000000060206000ULL, 0x00000000c040c000ULL,
	0x0000000302030000ULL, 0x0000000604060000ULL, 0x0000000e0a0e0000ULL, 0x0000001c141c0000ULL,
	0x0000003828380000ULL, 0x0000007050700000ULL, 0x0000006020600000ULL, 0x000000c040c00000ULL,
	0x0000030203000000ULL, 0x0000060406000000ULL, 0x00000e0a0e000000ULL, 0x00001c141c000000ULL,
	0x0000382838000000ULL, 0x0000705070000000ULL, 0x0000602060000000ULL, 0x0000c040c0000000ULL,
	0x0003020300000000ULL, 0x0006040600000000ULL, 0x000e0a0e00000000ULL, 0x001c141c00000000ULL,
	0x0038283800000000ULL, 0x0070507000000000ULL, 0x0060206000000000ULL, 0x00c040c000000000ULL,
	0x0002030000000000ULL, 0x0004060000000000ULL, 0x000a0e0000000000ULL, 0x00141c0000000000ULL,
	0x0028380000000000ULL, 0x0050700000000000ULL, 0x0020600000000000ULL, 0x0040c00000000000ULL,
	0x0203000000000000ULL, 0x0406000000000000ULL, 0x0a0e000000000000ULL, 0x141c000000000000ULL,
	0x2838000000000000ULL, 0x5070000000000000ULL, 0x2060000000000000ULL, 0x40c0000000000000ULL,
	0, 0 // <- hack for passing move & nomove
};

面白かったポイントは他にもあるのでたぶん続きます。

到達可能局面をSATソルバーに調べさせる

前回から引き続き、「オセロの初期局面から合法的に到達可能で、かつ合法手の数が34通りあるような局面は存在するのか?」を調べました。前回の最後に言及したように、flip関数をz3で書き下すことで、主問題そのものをソルバーに投げられるようにしました。結論から言うと、連言標準形でリテラル数57万くらい、節数370万くらいのSAT問題に帰着しましたが、最新のソルバーに投げても数時間では答えが出ないようでした。

コードは以下のリポジトリにあります。
retrospective-dfs-reversi/direct-dimacs at master · eukaryo/retrospective-dfs-reversi · GitHub

ループ/表引き不要なflip関数を探す

以前の記事で、edaxとそのAVX版に実装されているflip関数たちを紹介しました。それらは(forループで1マスずつ見るナイーブ版を除き)ほぼ全てが事前計算された表への添字アクセスを伴うものでした。これは実戦では有効ですが、これをz3で書き下すのは面倒そうなので別の実装を探しました。結果としてとあるリポジトリを見つけました。

github.com

この人のflip関数は、生のintel intrinsicsで書き下すと以下のようになります。

inline __m256i flipVertical(__m256i dbd) {
	return __m256i(_mm256_shuffle_epi8(dbd, _mm256_set_epi8(8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 4, 5, 6, 7)));
}
inline __m256i upper_bit(__m256i p) {
	p = _mm256_or_si256(p, _mm256_srli_epi64(p, 1));
	p = _mm256_or_si256(p, _mm256_srli_epi64(p, 2));
	p = _mm256_or_si256(p, _mm256_srli_epi64(p, 4));
	p = _mm256_andnot_si256(_mm256_srli_epi64(p, 1), p);
	p = flipVertical(p);
	p = _mm256_and_si256(p, _mm256_sub_epi64(_mm256_setzero_si256(), p));
	return flipVertical(p);
}
inline uint64_t hor(const __m256i lhs) {
	__m128i lhs_xz_yw = _mm_or_si128(_mm256_castsi256_si128(lhs), _mm256_extractf128_si256(lhs, 1));
	return uint64_t(_mm_cvtsi128_si64(lhs_xz_yw) | _mm_extract_epi64(lhs_xz_yw, 1));
}
uint64_t flip(uint64_t player, uint64_t opponent, int pos) {
	const __m256i pppp = _mm256_set1_epi64x(player);
	const __m256i oooo = _mm256_set1_epi64x(opponent);
	const __m256i OM = _mm256_and_si256(oooo, _mm256_set_epi64x(0xFFFFFFFFFFFFFFFFULL, 0x7E7E7E7E7E7E7E7EULL, 0x7E7E7E7E7E7E7E7EULL, 0x7E7E7E7E7E7E7E7EULL));

	__m256i flipped, outflank, mask;

	mask = _mm256_srli_epi64(_mm256_set_epi64x(0x0080808080808080ULL, 0x7F00000000000000ULL, 0x0102040810204000ULL, 0x0040201008040201ULL), 63 - pos);
	outflank = _mm256_and_si256(upper_bit(_mm256_andnot_si256(OM, mask)), pppp);
	flipped = _mm256_and_si256(_mm256_slli_epi64(_mm256_sub_epi64(_mm256_setzero_si256(), outflank), 1), mask);

	mask = _mm256_slli_epi64(_mm256_set_epi64x(0x0101010101010100ULL, 0x00000000000000FEULL, 0x0002040810204080ULL, 0x8040201008040200ULL), pos);
	outflank = _mm256_and_si256(_mm256_and_si256(mask, _mm256_add_epi64(_mm256_or_si256(OM, _mm256_andnot_si256(mask, _mm256_set1_epi8(0xFF))), _mm256_set1_epi64x(1))), pppp);
	flipped = _mm256_or_si256(flipped, _mm256_and_si256(_mm256_sub_epi64(outflank, _mm256_add_epi64(_mm256_cmpeq_epi64(outflank, _mm256_setzero_si256()), _mm256_set1_epi64x(1))), mask));
	return hor(flipped);
}

驚くべきことに表への添字アクセスが全く使われていません。pshufbは使われていますが、実質エンディアン変換なのでビットマスクとシフトで書き換えられます。これならSMTソルバーの制約の元でも簡単に書き下せます。あるいは余談ですが完全準同型暗号で暗号化された棋譜を受け取って暗号化された末端局面を秘匿計算するときにも役に立ちます。(何の意味もないでしょうが)

ちなみに、以前紹介したAVX2版edaxのflip関数をこれに差し替えてベンチマークを取ったところ、(fforum-40-59.obf 完全読みで)わたしの環境ではこれのほうが1%くらい遅かったです。逆アセンブルを見比べたらこれのほうが命令数が数個多かったので妥当だと思います。

変形版オセロのルールについて

序盤を全探索して確認したところ、盤上の石の数が17個から18個に増える瞬間までに最大2回のパスが発生しうることがわかりました。ゆえに、30→31に増える瞬間までに発生するパスの数はたかだか15回です。ということは、変形版オセロのルールとして

  • パスが2回連続してもそこで終局とはしない。
  • (パスを含み)42回目の着手をしたときに終局とする。
  • (終局時点の盤面を含み)登場した全ての盤面における合法手の数の最大値をその棋譜のスコアとする。
  • それ以外はオセロのルールと同じ。

と定義すると、n=33とn=34について、「スコアがn以上の棋譜が存在する」⇔「オセロの初期局面から合法的に到達可能で、かつ合法手の数がn通りあるような局面が存在する」といえます。なぜなら、(n=33を考えると)合法手が33通りあるためには石の数が31個以下でなければならないので、遅くとも「31個目の石を置いた直後にパスをしたその直後」までにスコア33でなければいけません。ゆえに、(石を置く手が26回と、途中のパスがたかだか15回と、最後のパス1回をあわせて)42回の着手を考慮すれば十分です。(同様にn=34は実は40回で十分なのですがさておきます)

n=33とn=34に注目している理由は、edaxのデータ構造がn<=32を仮定していて、それがバグかどうか知りたかったのがそもそものきっかけだったからです。

書き下す~問題作成まで

棋譜を探索すべき変数として、合法手生成とflipと盤面更新を42回繰り返す処理を書き下しました。(コードは下に貼ってありますしGitHubにもあります。)

https://github.com/eukaryo/retrospective-dfs-reversi/tree/master/direct-dimacs

ビットベクトル変数からなる問題は、z3の機能で連言標準形のSAT問題に変換できます。(本当はSAT以外の問題をSAT問題に変換する方法は色々あり、一分野を形成するほどなようですが、今回はとりあえず適当に変換しました)それをDIMACS CNFフォーマットに変換すれば、様々なSATソルバーに投げることができます。

今回は SAT Competitions のUNSAT部門(?)で優勝していて名前がかっこよかったCaDiCaLというやつに投げてみました。n=6とかn=18とかの簡単なケースでは数秒で解いてくれましたが、n-33とn=34はなかなか解けないようでした。

最後に

前回の記事で1600局面を生成しましたが、それらを目視したところ部分的によく似ていました。考えてみれば当然で、例えば外周28マスが全て空白で着手可能のとき、2,3行目のB~G列の12マスの配置はかなり限られるはずですよね。だとすると、もしもそれらの部分的なパターンがそれほど多くなくて、例えば100万通りとかを丁寧に場合分けすれば結論を出せるという状況であるなら、手作業では無理でもSATソルバーならその場合分けに相当する処理をやってくれるのではないかという淡い期待がありました。

コード

import time
import re
import sys

import z3
# z3をimportする簡単な方法:
# https://qiita.com/SatoshiTerasaki/items/476c9938479a4bfdda52


def solve(threshold):

    MAX_TURN = 42

    s = z3.Goal() # z3.Solver()
    
    bb_player = [z3.BitVec(f"bb_player_{i}", 64) for i in range(MAX_TURN+2)]
    bb_opponent = [z3.BitVec(f"bb_opponent_{i}", 64) for i in range(MAX_TURN+2)]
    pop_move = [z3.BitVec(f"pop_move_{i}", 64) for i in range(MAX_TURN+2)]

    directions1 = [z3.BitVec(f"directions1_{i}", 64) for i in range(4)]
    directions2 = [z3.BitVec(f"directions2_{i}", 64) for i in range(4)]

    flip_const1 = [z3.BitVec(f"flip_const1_{i}", 64) for i in range(4)]
    flip_const2 = [z3.BitVec(f"flip_const2_{i}", 64) for i in range(4)]
    flip_const3 = [z3.BitVec(f"flip_const3_{i}", 64) for i in range(4)]

    num_zero = z3.BitVec("num_zero", 64)
    num_one = z3.BitVec("num_one", 64)

    num_threshold = z3.BitVec("num_threshold", 64)

    s.add(
        bb_player[0] == 0x0000001008000000,
        bb_opponent[0] == 0x0000000810000000,
        pop_move[0] == 0,

        directions1[0] == 1,
        directions1[1] == 7,
        directions1[2] == 8,
        directions1[3] == 9,

        directions2[0] == 2,
        directions2[1] == 14,
        directions2[2] == 16,
        directions2[3] == 18,

        flip_const1[0] == 0xFFFFFFFFFFFFFFFF,
        flip_const1[1] == 0x7E7E7E7E7E7E7E7E,
        flip_const1[2] == 0x7E7E7E7E7E7E7E7E,
        flip_const1[3] == 0x7E7E7E7E7E7E7E7E,

        flip_const2[0] == 0x0080808080808080,
        flip_const2[1] == 0x7F00000000000000,
        flip_const2[2] == 0x0102040810204000,
        flip_const2[3] == 0x0040201008040201,

        flip_const3[0] == 0x0101010101010100,
        flip_const3[1] == 0x00000000000000FE,
        flip_const3[2] == 0x0002040810204080,
        flip_const3[3] == 0x8040201008040200,

        num_zero == 0,
        num_one == 1,
        num_threshold == threshold
    )

    for turn in range(MAX_TURN+1):
        bb_occupied = z3.BitVec(f"bb_occupied_{turn}", 64)
        bb_empty = z3.BitVec(f"bb_empty_{turn}", 64)
        masked_bb_opponent = z3.BitVec(f"masked_bb_opponent_{turn}", 64)
    
        movemask = [z3.BitVec(f"movemask_{turn}_{i}", 64) for i in range(4)]

        flip_l_0 = [z3.BitVec(f"flip_l_0_{turn}_{i}", 64) for i in range(4)]
        flip_r_0 = [z3.BitVec(f"flip_r_0_{turn}_{i}", 64) for i in range(4)]
        flip_l_1 = [z3.BitVec(f"flip_l_1_{turn}_{i}", 64) for i in range(4)]
        flip_r_1 = [z3.BitVec(f"flip_r_1_{turn}_{i}", 64) for i in range(4)]
        flip_l_2 = [z3.BitVec(f"flip_l_2_{turn}_{i}", 64) for i in range(4)]
        flip_r_2 = [z3.BitVec(f"flip_r_2_{turn}_{i}", 64) for i in range(4)]
        flip_l_3 = [z3.BitVec(f"flip_l_3_{turn}_{i}", 64) for i in range(4)]
        flip_r_3 = [z3.BitVec(f"flip_r_3_{turn}_{i}", 64) for i in range(4)]
        mask_l = [z3.BitVec(f"mask_l_{turn}_{i}", 64) for i in range(4)]
        mask_r = [z3.BitVec(f"mask_r_{turn}_{i}", 64) for i in range(4)]

        some_moves = [z3.BitVec(f"some_moves_{turn}_{i}", 64) for i in range(4)]
        all_moves = z3.BitVec(f"all_moves_{turn}", 64)
        pop_now = z3.BitVec(f"pop_now_{turn}", 64)

        popcnt_move_tmp1 = z3.BitVec(f"popcnt_move_tmp1_{turn}", 64)
        popcnt_move_tmp2 = z3.BitVec(f"popcnt_move_tmp2_{turn}", 64)

        s.add(
            bb_occupied == bb_player[turn] | bb_opponent[turn],
            bb_empty == bb_occupied ^ 0xFFFFFFFFFFFFFFFF,
            masked_bb_opponent == bb_opponent[turn] & 0x7E7E7E7E7E7E7E7E,
            movemask[0] == masked_bb_opponent,
            movemask[1] == masked_bb_opponent,
            movemask[2] == bb_opponent[turn],
            movemask[3] == masked_bb_opponent
        )
        s.add([z3.And(flip_l_0[i] == (movemask[i] & (bb_player[turn] << directions1[i]))) for i in range(4)])
        s.add([z3.And(flip_r_0[i] == (movemask[i] & z3.LShR(bb_player[turn], directions1[i]))) for i in range(4)])
        s.add([z3.And(flip_l_1[i] == (flip_l_0[i] | (movemask[i] & (flip_l_0[i] << directions1[i])))) for i in range(4)])
        s.add([z3.And(flip_r_1[i] == (flip_r_0[i] | (movemask[i] & z3.LShR(flip_r_0[i], directions1[i])))) for i in range(4)])
        s.add([z3.And(mask_l[i] == (movemask[i] & (movemask[i] << directions1[i]))) for i in range(4)])
        s.add([z3.And(mask_r[i] == (movemask[i] & z3.LShR(movemask[i], directions1[i]))) for i in range(4)])
        s.add([z3.And(flip_l_2[i] == (flip_l_1[i] | (mask_l[i] & (flip_l_1[i] << directions2[i])))) for i in range(4)])
        s.add([z3.And(flip_r_2[i] == (flip_r_1[i] | (mask_r[i] & z3.LShR(flip_r_1[i], directions2[i])))) for i in range(4)])
        s.add([z3.And(flip_l_3[i] == (flip_l_2[i] | (mask_l[i] & (flip_l_2[i] << directions2[i])))) for i in range(4)])
        s.add([z3.And(flip_r_3[i] == (flip_r_2[i] | (mask_r[i] & z3.LShR(flip_r_2[i], directions2[i])))) for i in range(4)])
        s.add([z3.And(some_moves[i] == ((flip_l_3[i] << directions1[i]) | z3.LShR(flip_r_3[i], directions1[i]))) for i in range(4)])
        s.add(
            all_moves == (some_moves[0] | some_moves[1] | some_moves[2] | some_moves[3]) & bb_empty,
            popcnt_move_tmp1 == all_moves - (z3.LShR(all_moves, 1) & 0x7777777777777777) - (z3.LShR(all_moves, 2) & 0x3333333333333333) - (z3.LShR(all_moves, 3) & 0x1111111111111111),
            popcnt_move_tmp2 == ((popcnt_move_tmp1 + z3.LShR(popcnt_move_tmp1, 4)) & 0x0F0F0F0F0F0F0F0F) * 0x0101010101010101,
            pop_now == z3.LShR(popcnt_move_tmp2, 56),
            pop_move[turn + 1] == z3.If(pop_now > pop_move[turn], pop_now, pop_move[turn])
        )

        move_onebit = z3.BitVec(f"move_onebit_{turn}", 64)
        move_bsf_tmp = [z3.BitVec(f"move_bsf_tmp_{turn}_{i}", 64) for i in range(9)]
        move_pos = z3.BitVec(f"move_pos_{turn}", 64)

        s.add(
            move_onebit & (move_onebit - 1) == 0,
            move_onebit & all_moves == move_onebit,

            move_bsf_tmp[0] == move_onebit - 1,
            move_bsf_tmp[1] == move_bsf_tmp[0] | z3.LShR(move_bsf_tmp[0], 1),
            move_bsf_tmp[2] == move_bsf_tmp[1] | z3.LShR(move_bsf_tmp[1], 2),
            move_bsf_tmp[3] == move_bsf_tmp[2] | z3.LShR(move_bsf_tmp[2], 4),
            move_bsf_tmp[4] == move_bsf_tmp[3] | z3.LShR(move_bsf_tmp[3], 8),
            move_bsf_tmp[5] == move_bsf_tmp[4] | z3.LShR(move_bsf_tmp[4], 16),
            move_bsf_tmp[6] == move_bsf_tmp[5] | z3.LShR(move_bsf_tmp[5], 32),
            move_bsf_tmp[7] == move_bsf_tmp[6] - (z3.LShR(move_bsf_tmp[6], 1) & 0x7777777777777777) - (z3.LShR(move_bsf_tmp[6], 2) & 0x3333333333333333) - (z3.LShR(move_bsf_tmp[6], 3) & 0x1111111111111111),
            move_bsf_tmp[8] == ((move_bsf_tmp[7] + z3.LShR(move_bsf_tmp[7], 4)) & 0x0F0F0F0F0F0F0F0F) * 0x0101010101010101,
            move_pos == z3.LShR(move_bsf_tmp[8], 56)
        )

        OM = [z3.BitVec(f"OM_{turn}_{i}", 64) for i in range(4)]
        mask1 = [z3.BitVec(f"mask1_{turn}_{i}", 64) for i in range(4)]
        upperbit_argument = [z3.BitVec(f"upperbit_argument_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp1 = [z3.BitVec(f"upperbit_tmp1_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp2 = [z3.BitVec(f"upperbit_tmp2_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp3 = [z3.BitVec(f"upperbit_tmp3_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp4 = [z3.BitVec(f"upperbit_tmp4_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp5 = [z3.BitVec(f"upperbit_tmp5_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp6 = [z3.BitVec(f"upperbit_tmp6_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp7 = [z3.BitVec(f"upperbit_tmp7_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp8 = [z3.BitVec(f"upperbit_tmp8_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp9 = [z3.BitVec(f"upperbit_tmp9_{turn}_{i}", 64) for i in range(4)]
        upperbit_tmp10 = [z3.BitVec(f"upperbit_tmp10_{turn}_{i}", 64) for i in range(4)]
        upperbit_result = [z3.BitVec(f"upperbit_result_{turn}_{i}", 64) for i in range(4)]
        outflank1 = [z3.BitVec(f"outflank1_{turn}_{i}", 64) for i in range(4)]
        flipped1 = [z3.BitVec(f"flipped1_{turn}_{i}", 64) for i in range(4)]

        mask2 = [z3.BitVec(f"mask2_{turn}_{i}", 64) for i in range(4)]
        outflank2 = [z3.BitVec(f"outflank2_{turn}_{i}", 64) for i in range(4)]
        nonzero = [z3.BitVec(f"nonzero_{turn}_{i}", 64) for i in range(4)]
        flipped2 = [z3.BitVec(f"flipped2_{turn}_{i}", 64) for i in range(4)]

        bb_flip = z3.BitVec(f"bb_flip_{turn}", 64)
        next_player_tmp = z3.BitVec(f"next_player_tmp_{turn}", 64)
        next_opponent_tmp = z3.BitVec(f"next_opponent_tmp_{turn}", 64)

        s.add([z3.And(OM[i] == (bb_opponent[turn] & flip_const1[i])) for i in range(4)])
        s.add([z3.And(mask1[i] == z3.LShR(flip_const2[i], 63 - move_pos)) for i in range(4)])

        s.add([z3.And(upperbit_argument[i] == (~OM[i]) & mask1[i]) for i in range(4)])
        s.add([z3.And(upperbit_tmp1[i] == upperbit_argument[i] | z3.LShR(upperbit_argument[i], 1)) for i in range(4)])
        s.add([z3.And(upperbit_tmp2[i] == upperbit_tmp1[i] | z3.LShR(upperbit_tmp1[i], 2)) for i in range(4)])
        s.add([z3.And(upperbit_tmp3[i] == upperbit_tmp2[i] | z3.LShR(upperbit_tmp2[i], 4)) for i in range(4)])
        s.add([z3.And(upperbit_tmp4[i] == (~z3.LShR(upperbit_tmp3[i], 1)) & upperbit_tmp3[i]) for i in range(4)])
        s.add([z3.And(upperbit_tmp5[i] == (upperbit_tmp4[i] << 32) | z3.LShR(upperbit_tmp4[i], 32)) for i in range(4)])
        s.add([z3.And(upperbit_tmp6[i] == ((upperbit_tmp5[i] & 0x0000FFFF0000FFFF) << 16) | z3.LShR(upperbit_tmp5[i] & 0xFFFF0000FFFF0000, 16)) for i in range(4)])
        s.add([z3.And(upperbit_tmp7[i] == ((upperbit_tmp6[i] & 0x00FF00FF00FF00FF) << 8) | z3.LShR(upperbit_tmp6[i] & 0xFF00FF00FF00FF00, 8)) for i in range(4)])
        s.add([z3.And(upperbit_tmp8[i] == upperbit_tmp7[i] & (-upperbit_tmp7[i])) for i in range(4)])
        s.add([z3.And(upperbit_tmp9[i] == (upperbit_tmp8[i] << 32) | z3.LShR(upperbit_tmp8[i], 32)) for i in range(4)])
        s.add([z3.And(upperbit_tmp10[i] == ((upperbit_tmp9[i] & 0x0000FFFF0000FFFF) << 16) | z3.LShR(upperbit_tmp9[i] & 0xFFFF0000FFFF0000, 16)) for i in range(4)])
        s.add([z3.And(upperbit_result[i] == ((upperbit_tmp10[i] & 0x00FF00FF00FF00FF) << 8) | z3.LShR(upperbit_tmp10[i] & 0xFF00FF00FF00FF00, 8)) for i in range(4)])

        s.add([z3.And(outflank1[i] == upperbit_result[i] & bb_player[turn]) for i in range(4)])
        s.add([z3.And(flipped1[i] == ((-outflank1[i]) << 1) & mask1[i]) for i in range(4)])

        s.add([z3.And(mask2[i] == flip_const3[i] << move_pos) for i in range(4)])
        s.add([z3.And(outflank2[i] == ((OM[i] | (~mask2[i])) + 1) & mask2[i] & bb_player[turn]) for i in range(4)])
        s.add([z3.And(nonzero[i] == z3.If(outflank2[i] == 0, num_zero, num_one)) for i in range(4)])
        s.add([z3.And(flipped2[i] == (outflank2[i] - nonzero[i]) & mask2[i]) for i in range(4)])
        s.add(
            bb_flip == flipped1[0] | flipped1[1] | flipped1[2] | flipped1[3] | flipped2[0] | flipped2[1] | flipped2[2] | flipped2[3],
            next_player_tmp == bb_opponent[turn] ^ bb_flip,
            next_opponent_tmp == bb_player[turn] ^ (bb_flip | move_onebit),
            bb_player[turn + 1] == z3.If(move_onebit == 0, bb_opponent[turn], next_player_tmp),
            bb_opponent[turn + 1] == z3.If(move_onebit == 0, bb_player[turn], next_opponent_tmp)
        )

    s.add(
        pop_move[MAX_TURN+1] >= num_threshold
    )

    t = z3.Then('simplify', 'bit-blast', 'solve-eqs', 'tseitin-cnf')
    subgoal = t(s)
    assert len(subgoal) == 1

    return subgoal[0]

if __name__ == "__main__":

    args = sys.argv
    if len(args) >= 3:
        print("invalid args")
        sys.exit(1)
    threshold = 34
    if len(args) == 2:
        threshold = int(args[1])

    print("threshold = " + str(threshold))
    print("start: solve")
    subgoal = solve(threshold)
    print("finish: solve")
    print("start: write a raw cnf (it may take tens of minutes...)")

    with open("cnf_subgoals_" + str(threshold) + ".txt", "w", encoding='utf-8') as f:

        for x in subgoal:
            s = "".join(str(x).split())
            f.write(s + "\n")

    print("finish: write a raw cnf")

    k_name_to_number = dict()
    clauses = []

    print("start: convert it into DIMACS CNF format")

    with open("cnf_subgoals_" + str(threshold) + ".txt", "r", encoding='utf-8') as f:

        line = f.readline()
        while line:
            answer = line.strip()
            answer = answer.replace(",","),")
            
            x = re.findall(r"k![0-9]+\)", answer)
            for s in x:
                if s not in k_name_to_number:
                    k_name_to_number[s] = str(len(k_name_to_number) + 1)
                answer = answer.replace(s, k_name_to_number[s])
            answer = answer.replace("Or(","")
            answer = answer.replace("Not","-")
            answer = answer.replace("(","")
            answer = answer.replace(")","")
            answer = answer.replace(","," ")
            clauses.append(answer)

            line = f.readline()

    with open("dimacs_cnf_" + str(threshold) + ".txt", "w", encoding='utf-8') as f:
        f.write("p cnf " + str(len(k_name_to_number)) + " " + str(len(clauses)) + "\n")
        for x in clauses:
            f.write(x + " 0\n")

    print("finished")