bitDPの話

[ビット演算テクニック Advent Calendar 2016](http://www.adventar.org/calendars/1403) 12日目の記事です。

みんな大好きbitDPの話です。編集距離 (Edit Distance: ED) をbitDPを使って計算する話を少々。

## 編集距離とは

まずは定義から。

編集距離 (ED; レーベンシュタイン距離: Levenshtein Distanceとも) は、2つの文字列の間で定義され、片方の文字列をもう一方の文字列に変換する場合の最小の {置換, 挿入, 削除} 回数を表す。

例えば、`abcdefg`と`abxdeg`の編集距離は2で、うちわけは置換が1回 (`c` → `x`)、挿入が0回、削除が1回 (`f`)となります。

## ふつうに解く

2次元のDPで解けることが知られています。ここでは漸化式だけ示します。

![ed.png](https://qiita-image-store.s3.amazonaws.com/0/153117/8434cb0c-11ef-2a1d-6371-20b078fbcf90.png)


## bitDPで解く

[Myersのアルゴリズム](http://dl.acm.org/citation.cfm?id=316550)という、ビット並列のアルゴリズムが知られています。日本語の解説記事も多く書かれており、知っている人も多いかと思います。

* [ビット並列アルゴリズムを使った編集距離](http://d.hatena.ne.jp/tsubosaka/20091102/1257174409)
* [C++: 編集距離を求めるアルゴリズム](http://handasse.blogspot.com/2009/04/c_29.html)
* [高速な編集距離の計算方法](http://qiita.com/aflc/items/f4299700471cc11f1d1c)

DP変数1セルに対し1bitのカラムを割り当て、64bitのレジスタを使って64セル並列に漸化式を計算するアルゴリズムとなっています (すごい)。もちろん上に示した漸化式のままでは1セル1bitで表現することはできないので、このアルゴリズムでは、セル値の代わりに隣接するセル値の差分を表現することで必要なビット幅を圧縮する工夫をしています。編集距離の漸化式ではこの差分が {-1, 0, 1} の3通りとなるため (証明は論文にあります)、pとmの2つの変数を用意し {(p, m)} = {(0, 1), (0, 0), (1, 0)} に対応させることで、1セル1bit x レジスタ4本 (p, mを縦横それぞれ) を実現しています。

1セル1bitなので、そこらのSIMD並列より並列度が高いことと、依存関係のチェインを加算命令のキャリーを使って解決しているところがうまい点です。

## LLCSも解けるはず

最長共通部分列長 (Length of Longest Common Substring: LLCS) も類似の2次元DPで解くことができます。これも隣り合うDP変数の差分が {0, 1} のいずれかとなるため、編集距離の場合と同様にbitDPに帰着させることができます。現在もっとも綺麗なアルゴリズムは[Hyyröによるもの](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.59.6975&rep=rep1&type=pdf)のようです。発想はEDのアルゴリズムと同様ですが、片方の変数 (水平差分) を反転させることで演算の回数を減らしています。

## とりあえず実装しよう

このアルゴリズムが強力なのは、必要なデータをすべてレジスタにのせたまま64セルの更新が可能な点です。この特徴を損ねないよう注意深く記述します。

しかし任意長の入力を受け取れるようにするのはなかなか面倒なので、入力は最長64文字までという仕様にします。 `uint64_t editdist(char const *a, uint64_t alen, char const *b, uint64_t blen);` とすればそこそこに実用的に使えるかなと。例えば文章中一部を切り出してEDを計算するような場合に、zero-copyで実現できますね。

### コード

[完全版はこちらに](https://gist.github.com/ocxtal/8c87a60d12da58dc275f62a7cbe54807)

```C
uint64_t editdist(
	char const *a,
	uint64_t alen,
	char const *b,
	uint64_t blen)
{
	v32i8_t a1 = _loadu_v32i8(a), a2 = _loadu_v32i8(a + 32);
	uint64_t ph = 0xffffffffffffffff, mh = 0;
	for(char const *p = b; p < b + blen; p++) {
		v32i8_t b = _set_v32i8(*p);
		uint64_t eq = (_m(_eq_v32i8(b, a2))<<32) | _m(_eq_v32i8(b, a1));
		uint64_t xh = eq | mh, xv = (((eq & ph) + ph) ^ ph) | eq;
		uint64_t pv = mh | ~(xv | ph), mv = ph & xv;
		pv<<=1; mv<<=1; pv |= 0x01ULL;
		ph = mv | ~(xh | pv); mh = pv & xh;
	}
	uint64_t mask = (0x01ULL<<alen) - 1;
	return(blen + popcnt(ph & mask) - popcnt(mh & mask));
}
```

```C
uint64_t llcs(
	char const *a,
	uint64_t alen,
	char const *b,
	uint64_t blen)
{
	v32i8_t a1 = _loadu_v32i8(a), a2 = _loadu_v32i8(a + 32);
	uint64_t h = 0xffffffffffffffff;
	for(char const *p = b; p < b + blen; p++) {
		v32i8_t b = _set_v32i8(*p);
		uint64_t eq = (_m(_eq_v32i8(b, a2))<<32) | _m(_eq_v32i8(b, a1));
		uint64_t v = h & eq; h = (h + v) | (h - v);
	}
	uint64_t mask = (0x01ULL<<alen) - 1;
	return(alen - popcnt(h & mask));
}
```

### コンパイル結果

clang-3.9 -O3 -msse4.1 on macです。EDの実装の再内のループのみ示します。意図通りメモリへの退避/リロードはなくせました。

```asm
LBB0_2:                                 ## =>This Inner Loop Header: Depth=1
	movzbl	(%rdx), %ebx
	movd	%ebx, %xmm5
	pshufb	%xmm4, %xmm5
	movdqa	%xmm5, %xmm6
	pcmpeqb	%xmm2, %xmm6
	pmovmskb	%xmm6, %ebx
	movdqa	%xmm5, %xmm6
	pcmpeqb	%xmm3, %xmm6
	pmovmskb	%xmm6, %ecx
	shll	$16, %ecx
	orl	%ebx, %ecx
	shlq	$32, %rcx
	movdqa	%xmm5, %xmm6
	pcmpeqb	%xmm0, %xmm6
	pmovmskb	%xmm6, %edi
	pcmpeqb	%xmm1, %xmm5
	pmovmskb	%xmm5, %ebx
	shll	$16, %ebx
	orl	%edi, %ebx
	orq	%rcx, %rbx
	movq	%rbx, %rcx
	andq	%rax, %rcx
	addq	%rax, %rcx
	xorq	%rax, %rcx
	orq	%rbx, %rcx
	orq	%r13, %rbx
	movq	%rcx, %rdi
	orq	%rax, %rdi
	notq	%rdi
	orq	%r13, %rdi
	andq	%rax, %rcx
	addq	%rcx, %rcx
	leaq	1(%rdi,%rdi), %rdi
	movq	%rdi, %rax
	orq	%rbx, %rax
	notq	%rax
	orq	%rcx, %rax
	andq	%rdi, %rbx
	incq	%rdx
	cmpq	%r12, %rdx
	movq	%rbx, %r13
	jb	LBB0_2
```

## おわりに

* ベンチマークとらないとね
* 配列のロード、ページ境界をまたぐとSEGVする可能性があるね......
* UTF-8を食えるようにするにはどうすればいいだろうか...

というわけで、これらの問題が解決したら随時アップデートします。

明日 (今日) はそすうさん ([@_primenumber](http://www.adventar.org/users/7698)) の並列Bit Scan Reverseだそうです。

## 参考文献

* [Gene Myers, A fast bit-vector algorithm for approximate string matching based on dynamic programming., Journal of the ACM (JACM). 1999](http://dl.acm.org/citation.cfm?id=316550)
* [Heikki Hyyrö, Bit-Parallel LCS-length Computation Revisited., Proc. 15th Australasian Workshop on Combinatorial Algorithms (AWOCA 2004). 2004.](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.59.6975&rep=rep1&type=pdf)