重ならないビット列を3進数とみなして2進数にする

2進数(binary number)は0と1だけからなりますが、3進数(ternary number)は0と1と2からなります。2進40桁以下の符号なし整数を受け取って、それを『偶然0と1だけでかけるような3進数だった』と解釈して、その値を2進数にして返す関数は以下のように書けます。なぜ40桁以下かというと、返り値が64桁に収まるギリギリの桁数だからです。

uint64_t ternarize(uint64_t x) {
	assert(x <= 0x0000'00FF'FFFF'FFFFULL);
	uint64_t answer = 0;
	for (uint64_t value = 1; x; x /= 2, value *= 3) {
		if (x & 1ULL) {
			answer += value;
		}
	}
	return answer;
}

一般の3進数は0と1と2からなります。言い換えると各桁が2進2桁の数だとみなせます。そのうえで、n桁の3進数を"ビットプレーン分解"してn桁の2進数2つ $(l,u)$ の組に分解することができます。曖昧に"ビットプレーン分解"と書きましたが、例えば $12012211_{ternary}$ は $l=0b10010011$ と $u=0b01001100$ に分解されます。

このように分解された $(l,u)$ を入力として、それを3進数と解釈したうえで2進数にして返す関数は以下のように書けます。さっきと同様、元の3進数は3進40桁以下だとします。（返り値をuint64_tに収めるため）

uint64_t ternarize_naive1(uint64_t u, uint64_t l) {

	assert((u | l) <= 0x0000'00FF'FFFF'FFFFULL);
	assert((u & l) == 0);

	return ternarize(u) * 2 + ternarize(l);
}
uint64_t ternarize_naive2(uint64_t u, uint64_t l) {

	assert((u | l) <= 0x0000'00FF'FFFF'FFFFULL);
	assert((u & l) == 0);

	uint64_t answer = 0;

	for (uint64_t value = 1; u | l; u /= 2, l /= 2, value *= 3) {
		if (u & 1ULL) {
			answer += value * 2;
		}
		else if (l & 1ULL) {
			answer += value;
		}
	}

	return answer;
}

わかりやすさのため2通り書きましたが、常に同じ値を返します。

pshufb

この処理はSSEを使って高速化できます。具体的には、

各引数を4bit区切り10個とみなして、__m128iの8bit領域16箇所に順番に詰める。
4bit数を「3進数とみなして2進数に変換する」処理はpshufbで一発でできる。
$u*2+l$ をやる。すると結局、81進数10桁みたいなやつになる。
隣り合ってるやつを足す処理を $\lceil log_{2}10\rceil$ 回繰り返すと、 $3^{64}$ 進数1桁になる。

といった感じでやります。

uint64_t ternarize_sse(uint64_t u, uint64_t l) {

	assert((u | l) <= 0x0000'00FF'FFFF'FFFFULL);
	assert((u & l) == 0);

	const __m128i ternary_binary_table = _mm_set_epi8(40, 39, 37, 36, 31, 30, 28, 27, 13, 12, 10, 9, 4, 3, 1, 0);
	const __m128i bitmask0F = _mm_set1_epi8(0x0F);

	const __m128i lu_lo = _mm_set_epi64x(l, u);
	const __m128i lu_hi = _mm_srli_epi64(lu_lo, 4);

	const __m128i u_4bits = _mm_and_si128(_mm_unpacklo_epi8(lu_lo, lu_hi), bitmask0F);
	const __m128i l_4bits = _mm_and_si128(_mm_unpackhi_epi8(lu_lo, lu_hi), bitmask0F);

	//この時点でu_4bitsは、uを4bit区切りにして__m128iの8bit領域16箇所に順番に詰めた形になっている。l_4bitsも同様。

	const __m128i ternarized_u_4bits = _mm_shuffle_epi8(ternary_binary_table, u_4bits);
	const __m128i ternarized_l_4bits = _mm_shuffle_epi8(ternary_binary_table, l_4bits);
	const __m128i answer_base_3p4_in_i8s = _mm_add_epi8(_mm_add_epi8(ternarized_u_4bits, ternarized_u_4bits), ternarized_l_4bits);

	//この時点でanswer_base_3p4_in_i8sは、真の値を"81進数"にした状態で、各桁の値を__m128iの8bit領域16箇所に順番に詰めた形になっている。

	const __m128i tmp_mask8_lo = _mm_set_epi8(0, 0xFF, 0, 0xFF, 0, 0xFF, 0, 0xFF, 0, 0xFF, 0, 0xFF, 0, 0xFF, 0, 0xFF);
	const __m128i tmp_shuf8_hi = _mm_set_epi8(0xFF, 15, 0xFF, 13, 0xFF, 11, 0xFF, 9, 0xFF, 7, 0xFF, 5, 0xFF, 3, 0xFF, 1);
	const __m128i i16_81s = _mm_set1_epi16(81);
	const __m128i answer_tmp8_lo = _mm_and_si128(answer_base_3p4_in_i8s, tmp_mask8_lo);
	const __m128i answer_tmp8_hi = _mm_shuffle_epi8(answer_base_3p4_in_i8s, tmp_shuf8_hi);
	const __m128i answer_base3p8_in_i16s = _mm_add_epi16(answer_tmp8_lo, _mm_mullo_epi16(answer_tmp8_hi, i16_81s));

	//この時点でanswer_base3p8_in_i16sは、真の値を"6561進数"(6561=81*81=3^8)にした状態で、各桁の値を__m128iの16bit領域8箇所に順番に詰めた形になっている。

	const __m128i tmp_mask16_lo = _mm_set_epi8(0, 0, 0xFF, 0xFF, 0, 0, 0xFF, 0xFF, 0, 0, 0xFF, 0xFF, 0, 0, 0xFF, 0xFF);
	const __m128i tmp_shuf16_hi = _mm_set_epi8(0xFF, 0xFF, 15, 14, 0xFF, 0xFF, 11, 10, 0xFF, 0xFF, 7, 6, 0xFF, 0xFF, 3, 2);
	const __m128i i32_6561s = _mm_set1_epi32(6561);
	const __m128i answer_tmp16_lo = _mm_and_si128(answer_base3p8_in_i16s, tmp_mask16_lo);
	const __m128i answer_tmp16_hi = _mm_shuffle_epi8(answer_base3p8_in_i16s, tmp_shuf16_hi);
	const __m128i answer_base3p16_in_i32s = _mm_add_epi32(answer_tmp16_lo, _mm_mullo_epi32(answer_tmp16_hi, i32_6561s));

	//この時点でanswer_base3p16_in_i32sは、真の値を"43046721進数"(43046721=6561*6561=3^16)にした状態で、各桁の値を__m128iの32bit領域4箇所に順番に詰めた形になっている。

	const __m128i tmp_mask32_lo = _mm_set_epi8(0, 0, 0, 0, 0xFF, 0xFF, 0xFF, 0xFF, 0, 0, 0, 0, 0xFF, 0xFF, 0xFF, 0xFF);
	const __m128i tmp_shuf32_hi = _mm_set_epi8(0xFF, 0xFF, 0xFF, 0xFF, 15, 14, 13, 12, 0xFF, 0xFF, 0xFF, 0xFF, 7, 6, 5, 4);
	const __m128i i64_43046721s = _mm_set1_epi64x(43046721);
	const __m128i answer_tmp32_lo = _mm_and_si128(answer_base3p16_in_i32s, tmp_mask32_lo);
	const __m128i answer_tmp32_hi = _mm_shuffle_epi8(answer_base3p16_in_i32s, tmp_shuf32_hi);
	const __m128i answer_base3p32_in_i64s = _mm_add_epi64(answer_tmp32_lo, _mm_mul_epi32(answer_tmp32_hi, i64_43046721s));

	//この時点でanswer_base3p32_in_i64sは、真の値を"1853020188851841進数"にした状態で、各桁の値を__m128iの64bit領域2箇所に順番に詰めた形になっている。
	//ちなみに1853020188851841<2^63なのでepi64で扱っても問題ない。

	alignas(16) uint64_t answer[2] = {};
	_mm_storeu_si128((__m128i*)answer, answer_base3p32_in_i64s);
	return answer[1] * 1853020188851841ULL + answer[0];
}

64桁の場合

上記の処理は引数を40桁以下に限定していましたが、3進64桁の数を扱いたい場合もあります。単純に24桁と40桁に分けて処理する関数は以下のように書けます。これは、上記のSSEのやつを愚直にAVX2化して同時に計算できます。ちなみに3進24桁は2進だと39桁に収まります。AVX2版は長いのでここには貼りませんが、GitHubで公開しています。

void ternarize_naive1_full(const uint64_t u, const uint64_t l, uint64_t *answer_l, uint64_t *answer_u) {

	assert((u & l) == 0);

	*answer_l = ternarize(u & 0x0000'00FF'FFFF'FFFFULL) * 2 + ternarize(l & 0x0000'00FF'FFFF'FFFFULL);
	*answer_u = ternarize(u >> 40ULL) * 2 + ternarize(l >> 40ULL);
}

github.com

ベンチマーク結果

遺伝研のXeonGold6136 (Skylake, 3.0GHz) と EPYC7702 (Zen2, 2.0GHz) で実験しました。gcc8.2.0でg++ -std=c++1y -O3 -flto -march=nativeでコンパイルしました。xorshift64で乱数生成して計算するのを $2^{30}$ 回繰り返しました。negative controlとは本体の計算をやらずに乱数生成だけやった場合の時間です。

40bit版:

	XeonGold6136	EPYC7702
negative control (s)	2.068	2.404
naive1 (s)	49.598	70.063
naive2 (s)	33.739	45.706
SSE (s)	7.766	5.175

64bit版:

	XeonGold6136	EPYC7702
negative control (s)	2.171	2.002
naive1 (s)	76.116	97.503
naive2 (s)	52.443	62.469
AVX2 (s)	10.315	6.066

まとめ

EPYCとかいうのクロック周波数が2/3しかないのに速くてすごい。

オセロの盤面をハッシュテーブルに入れるとして、ハッシュ衝突を絶対に防ぎたいとします。普通に考えると手番側の石のビットボードとその相手の石のビットボードをコンカチした128bitをvalueに載せてチェックするでしょう。今回の方法はオセロの盤面に対する103bitの完全ハッシュとみなせます。一様にバラけるわけではないので、完全ハッシュよりは可逆圧縮と呼ぶべきかもしれません。この103bitをvalueに載せておけばハッシュ衝突を完全に防げて、空いた25bitに別のなにかを書き込む余地ができます。