ハッカーのたのしみ第7章7.5の解説とちょっとした改善

algorithm-study/hackers_delight_7_5 at master · eukaryo/algorithm-study · GitHub

ハッカーのたのしみ第7章7.5では、「ワードの一般的順列演算」と呼ばれている演算の実装方法について論じられています。「ワードの一般的順列演算」は、愚直に書くと以下のように書けます。引数xのi桁目のビットをs[i]桁目に移すというものです。

uint64_t permutation_naive(const uint64_t x, const int s[64]) {

	//引数sは0～63の順列だとする。
	for (int i = 0; i < 64; ++i)
		assert(0 <= s[i] && s[i] < 64);
	for (int i = 0; i < 64; ++i)for (int j = i + 1; j < 64; ++j)
		assert(s[i] != s[j]);

	uint64_t answer = 0;

	for (int i = 0; i < 64; ++i) {
		const uint64_t bit = 1ULL << i;
		if (x & bit)answer |= 1ULL << s[i];
	}

	return answer;
}

さて、ハッカーのたのしみ第7章の7.4では、compressと呼ばれる別の演算について説明されています。これはx86ではpextという名前で実装されています。（というか、ハッカーのたのしみがpextについて論じていることを最近まで忘れていました）

「ワードの一般的順列演算」を実装するにあたり、そのcompress演算を使ってsag関数という部品をまず定義します。sag関数は、pextのmaskにて立っている場所のビットを上位桁側にまとめて、それ以外を下位桁側にまとめるというものです。言い換えると、xとmaskを1ビットの64要素の配列として見たとき、maskの値に基づきxを安定ソートしていると解釈できます。

uint64_t compress(uint64_t x, uint64_t mask) {
	return _pext_u64(x, mask);
}
uint64_t sag(uint64_t x, uint64_t mask) {
	return (compress(x, mask) << (_mm_popcnt_u64(~mask)))
		| compress(x, ~mask);
}

次に、順列を定義している64要素の配列を「ビットごとに行列転置」して、64ビット整数6個に変換します。

uint64_t p[6] = {};//さしあたりグローバル変数

void init_p_array(const int x[64]) {

	//引数xは0～63の順列だとする。
	for (int i = 0; i < 64; ++i)
		assert(0 <= x[i] && x[i] < 64);
	for (int i = 0; i < 64; ++i)for (int j = i + 1; j < 64; ++j)
		assert(x[i] != x[j]);

	for (int i = 0; i < 6; ++i)p[i] = 0;

	//2進6桁の値64個からなる配列xを「ビットごとに行列転置」して、
	//2進64桁の値6個の配列pに格納する。
	for (int i = 0; i < 64; ++i) {
		for (uint64_t b = x[i], j = 0; b; ++j, b >>= 1) {
			if (b & 1ULL) {
				p[j] |= 1ULL << i;
			}
		}
	}

	//ハッカーのたのしみ133ページの事前計算
	p[1] = sag(p[1], p[0]);
	p[2] = sag(sag(p[2], p[0]), p[1]);
	p[3] = sag(sag(sag(p[3], p[0]), p[1]), p[2]);
	p[4] = sag(sag(sag(sag(p[4], p[0]), p[1]), p[2]), p[3]);
	p[5] = sag(sag(sag(sag(sag(p[5], p[0]), p[1]), p[2]), p[3]), p[4]);
}

所望の順列を意味する配列をinit_p_arrayの引数に与えると、グローバル変数の配列pに、いい感じのマジックナンバーが格納されます。

そのうえで、sag関数を逐次的に6回呼ぶことで「ワードの一般的順列演算」ができます。（ハッカーのたのしみでは32bitなので5回でしたが、64bitだと6回になります）

uint64_t permutation(uint64_t x) {
	x = sag(x, p[0]);
	x = sag(x, p[1]);
	x = sag(x, p[2]);
	x = sag(x, p[3]);
	x = sag(x, p[4]);
	x = sag(x, p[5]);
	return x;
}

ハッカーのたのしみにも書かれていることですが、このpermutation関数は実質的にradix sortをやっています。これは安定ソートです。

ここからはハッカーのたのしみに書かれていないのですが、所望の順列が良い性質を持っていれば、より良いマジックナンバーを生成できて、permutation関数内でsag関数を実行する回数を5回以下にできます。まず以下の関数によって、所望の順列を意味する配列xを別の配列x_stableに変換します。

void func(const int x[64], int x_stable[64]) {
	int rep = 0;
	for (int count = 0; count < 64; ++rep) {
		for (int i = 0; i < 64; ++i) {
			if (x[i] == count) {
				x_stable[i] = rep;
				count++;
			}
		}
	}
}

以下の表はランダムなxにこの関数を作用させた例です。

添字	0	1	2	3	4	5	6	7
x[]	3	0	4	7	1	5	6	2
x_stable[]	1	0	1	2	0	1	1	0

配列xとx_stableは「ソート対象物を順序付けるもの」なので、ソートが安定であればどちらの順序付けに基づいても同一のソート結果をもたらすのがポイントです。

前述のinit_p_array関数は以下のように書き換える必要があります。

void init_p_array_better(const int x[64]) {

	for (int i = 0; i < 64; ++i)
		assert(0 <= x[i] && x[i] < 64);
	for (int i = 0; i < 64; ++i)for (int j = i + 1; j < 64; ++j)
		assert(x[i] != x[j]);

	int x_stable[64] = {};
	func(x, x_stable);

	for (int i = 0; i < 6; ++i)p[i] = 0;
	for (int i = 0; i < 64; ++i) {
		for (uint64_t b = x_stable[i], j = 0; b; ++j, b >>= 1) {
			if (b & 1ULL) {
				p[j] |= 1ULL << i;
			}
		}
	}

	p[1] = sag(p[1], p[0]);
	p[2] = sag(sag(p[2], p[0]), p[1]);
	p[3] = sag(sag(sag(p[3], p[0]), p[1]), p[2]);
	p[4] = sag(sag(sag(sag(p[4], p[0]), p[1]), p[2]), p[3]);
	p[5] = sag(sag(sag(sag(sag(p[5], p[0]), p[1]), p[2]), p[3]), p[4]);
}

permutation関数内でsag関数を実行すべき回数は、x_stable配列の最大値をmとすると、ceiling(log2(m+1))回となります。

例えば、所望の順列が「16bitごとに逆転させる」というものだったならば、sagは2回で済みます。配列pがp[2]以降全てゼロになるのです。

int main(void) {

	int x[64];
	for (int i = 0; i < 4; ++i)for (int j = 0; j < 16; ++j) {
		x[i * 16 + j] = (3 - i) * 16 + j;
	}

	init_p_array(x);

	std::cout << "p = (";
	for (int i = 0; i < 6; ++i)std::cout << std::hex << p[i] << (i != 5 ? ", " : ")");
	std::cout << std::endl;

	const uint64_t a = 0x1234'5678'9ABC'DEF0ULL;
	const uint64_t b0 = permutation_naive(a, x);
	const uint64_t b1 = permutation(a);

	init_p_array_better(x);

	std::cout << "p = (";
	for (int i = 0; i < 6; ++i)std::cout << std::hex << p[i] << (i != 5 ? ", " : ")");
	std::cout << std::endl;

	const uint64_t b2 = permutation(a);

	std::cout << "a  = " << a << std::endl;
	std::cout << "b0 = " << b0 << std::endl;
	std::cout << "b1 = " << b1 << std::endl;
	std::cout << "b2 = " << b2 << std::endl;

	return 0;
}

以上のプログラムを走らせると、出力は以下の通りになります。

p = (aaaaaaaaaaaaaaaa, aaaaaaaaaaaaaaaa, aaaaaaaaaaaaaaaa, aaaaaaaaaaaaaaaa, 5555555555555555, 5555555555555555)
p = (ffff0000ffff, ffff0000ffff, 0, 0, 0, 0)
a  = 123456789abcdef0
b0 = def09abc56781234
b1 = def09abc56781234
b2 = def09abc56781234

「16bitごとに逆転させる」というのはデモのための例に過ぎなくて、それ自体はビットマスクとシフトとかでやったほうが当然速いです。

今回のsagベースのアプローチに加えて、ビットシフトとかもアリだとしたときに、所望の順列への変換を実現する最速の機械語列を本当は求めたいのですが、これは難しい問題です。