bitDPの話 #ビット演算

ビット演算テクニック Advent Calendar 2016 12日目の記事です。

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

編集距離とは

まずは定義から。

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

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

ふつうに解く

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

bitDPで解く

Myersのアルゴリズムという、ビット並列のアルゴリズムが知られています。日本語の解説記事も多く書かれており、知っている人も多いかと思います。

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öによるもののようです。発想はEDのアルゴリズムと同様ですが、片方の変数 (水平差分) を反転させることで演算の回数を減らしています。

とりあえず実装しよう

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

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

コード

完全版はこちらに

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));
}

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の実装の再内のループのみ示します。意図通りメモリへの退避/リロードはなくせました。

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) の並列Bit Scan Reverseだそうです。

bitDPの話