第16章
Boyer–Moore アルゴリズム
(やさしい版)
Pearls of Functional Algorithm Design(関数プログラミングによるアルゴリズム設計の真珠)
どんな問題?
ある短い文字列(パターン)が、別の長い文字列(テキスト)の中でどこに現れているかをすべて見つけたい、というのがこの章のテーマです。ワープロの「検索」機能を思い浮かべてください。
たとえば、パターンが "abcab"、テキストが "ababcabcab" のとき、パターンが終わる位置(1 から数える)は 7 と 10 です。
matches "abcab" "ababcabcab" = [7, 10]
この章で扱うこと
素朴に照合すると、テキストの長さ n、パターンの長さ m に対して Θ(mn) の時間がかかります。この章では、有名な Boyer–Moore(BM)アルゴリズムを導き、時間を Θ(m + n) まで縮めます。次章では同じ計算量を持つ Knuth–Morris–Pratt(KMP)アルゴリズムを扱います。
まずは仕様
Haskell 風に書くと、matches は次のように仕様化できます。
matches :: Eq a ⇒ [a] → [a] → [Int]
matches ws = map length · filter (endswith ws) · inits
Dart
// 仕様:xs のすべての接頭辞のうち ws を接尾辞に持つものの長さ
List<int> matches<T>(List<T> ws, List<T> xs) {
final List<int> result = [];
for (int n = 0; n <= xs.length; n++) {
final List<T> prefix = xs.sublist(0, n);
if (endswith(ws, prefix)) result.add(n);
}
return result;
}
bool endswith<T>(List<T> ws, List<T> xs) {
if (ws.length > xs.length) return false;
final int start = xs.length - ws.length;
for (int i = 0; i < ws.length; i++) {
if (ws[i] != xs[start + i]) return false;
}
return true;
}
- inits xs:xs の接頭辞を長さの短い順に全部並べたリスト
- endswith ws xs:ws が xs の接尾辞(末尾)になっているか
- matches ws xs:ws がちょうど位置 p で終わるような p のリスト
照合は要素どうしの等号 (==) しか使わないので、要素の型は何でも構いません(つまり多相です)。等号判定を定数時間と仮定すると、素朴な実装の最悪計算量は Θ(mn) になります。
scan の補題 — 高速化のカギ
この補題が主役
inits が絡む問題では、次の scan の補題がほとんどいつも主役になります。
map (foldl op e) · inits = scanl op e
ポイント
左辺は各接頭辞ごとに畳み込みをやり直すので Θ(n2) 回の op 評価が必要です。右辺の scanl は途中結果を使い回すので Θ(n) 回で済みます。この差が、線形時間アルゴリズムの源です。
filter を map の隣に運ぶ
matches の定義には map だけでなく filter もあるので、まず次の法則で並び替えます。
map f · filter p = map fst · filter snd · map (fork (f, p)) (16.1)
ここで fork (f, p) x = (f x, p x) です。これを matches に適用すると、
matches ws
= map fst · filter snd · map (fork (length, endswith ws)) · inits
となり、map が inits の隣に来ました。あとは中身の fork (length, endswith ws) を foldl の形にできれば、scan の補題が使えます。
foldl の形にできるか?
length は簡単で、length = foldl count 0(count n x = n + 1)で書けます。同じように endswith ws も
endswith ws = foldl op e (16.2)
という形にできれば、foldl のタプリング則
fork (foldl op1 e1, foldl op2 e2) = foldl op (e1, e2)
を使って一気に整理できます。この式は「2 つの畳み込みを 1 本の畳み込みにまとめられる」という便利な法則です。うまくいけば
matches ws = map fst · filter snd · scanl step (0, e)
という線形時間の形になります。
壁にぶつかる
ところが、endswith ws は真偽値(Bool)を返すだけなので、そのままでは foldl にはできません。畳み込みの途中状態として保持すべき情報が足りないのです。
合成関数として書き直す
次善の策として、endswith ws を「畳み込みで状態を作り、その状態から真偽値を取り出す」形に分解します。
endswith ws = p · foldl op e (16.3)
これに合わせて (16.1) を少し一般化した法則
map f · filter (p · g) = map fst · filter (p · snd) · map (fork (f, g)) (16.4)
を使うと、
matches ws = map fst · filter (p · snd) · scanl step (0, e)
となります。p と op が償却定数時間(ならして定数時間)で計算できれば、matches は線形時間です。
endswith の 2 つの定義 — 分岐点
あとは p、op、e を決めればよいのですが、そもそも endswith にはもっともらしい定義が 2 通りあります。
endswith ws xs = reverse ws ⊑ reverse xs
endswith ws xs = ws ∈ tails xs
ここで us ⊑ vs は「us が vs の接頭辞である」の意味です。「ws が末尾に来る」ことは「両方を逆順にしたときに ws の逆が先頭に来る」ことと同じなので、上の 1 番目の定義は成り立ちます。接頭辞判定は素直に書けます。
[ ] ⊑ vs = True
(u : us) ⊑ [ ] = False
(u : us) ⊑ (v : vs) = (u == v ∧ us ⊑ vs)
Dart
// us が vs の接頭辞か(us ⊑ vs)
bool isPrefix<T>(List<T> us, List<T> vs) {
if (us.length > vs.length) return false;
for (int i = 0; i < us.length; i++) {
if (us[i] != vs[i]) return false;
}
return true;
}
大事な分岐
2 つの定義は結果としては同じ関数を表します。しかし、
導出の道筋を決めるのは「関数の中身」ではなく「式の形」です。
- 1 番目の定義(reverse と接頭辞)を選ぶと → BM アルゴリズム
- 2 番目の定義(tails と要素判定)を選ぶと → KMP アルゴリズム
この章では 1 番目の道を進みます。
BM アルゴリズムの基本形
右から左に照合するのはなぜ?
endswith の 1 番目の定義は、次のように合成として書けます。
endswith ws = (reverse ws ⊑ ) · reverse
これを (16.4) に当てはめると
matches ws
= map fst · filter ((sw ⊑) · snd) · map (fork (length, reverse)) · inits
where sw = reverse ws
となります。ここで reverse = foldl (flip (:)) [ ] なので、再びタプリング則と scan の補題を使って、
matches ws = map fst · filter ((sw ⊑) · snd) · scanl step (0, [ ])
where sw = reverse ws
step (n, sx) x = (n + 1, x : sx)
Dart
// BM の基本形:窓 (n, sx) を左から作り、sw ⊑ sx で判定
// sx は「テキストを 1 文字ずつ先頭に積んだ逆順リスト」
List<int> matchesBasic<T>(List<T> ws, List<T> xs) {
final List<T> sw = ws.reversed.toList();
final List<int> result = [];
int n = 0;
List<T> sx = [];
// 初期窓 (0, []) から順に scan
if (isPrefix(sw, sx)) result.add(n);
for (final T x in xs) {
n = n + 1;
sx = [x, ...sx]; // step: x を先頭に付ける
if (isPrefix(sw, sx)) result.add(n);
}
return result;
}
とまとまります。これが BM アルゴリズムの基本形です。
「窓」と「シフト」
scanl はテキストを左から 1 文字ずつ読みながら、その時点までの反転(sx)と現在位置(n)を組にした窓 (window) を次々に作ります。窓は 1 文字ずつずれるので、1 段階ごとのシフト (shift) は 1 です。各窓では、パターン ws を右から左にあてて照合します。この「右から左」という向きは、endswith の定義選びから自然に出てきたものです。
シフトを大きく — 最悪計算量を減らす
なぜまだ遅いのか
基本形はまだ最悪 Ω(mn) 時間です。判定 (sw ⊑) が最悪 Ω(m) かかるからです(以下 m = length ws、m ≠ 0 と仮定します)。たとえばパターンが "aaaaa"、テキストが "aaaaa...a" のような場合が最悪です。
照合失敗を活かして飛ばす
そこで、現在の窓で「どこまで一致したか」を利用して、続く何個かの窓をまとめて飛ばすことを考えます。飛ばした窓に照合結果があってはならないので、飛ばし方には注意が必要です。
llcp sw sx を「sw と sx の最長共通接頭辞の長さ」とします(前章で登場した関数です)。sw ⊑ sx は llcp sw sx = m と同値です。
いま現在の窓 (n, sx) で i = llcp sw sx が分かっているとしましょう。次にマッチしうる窓の位置 n + k(0 < k ≤ m)を考えます。もし次の窓 (n+k, ys ++ sx)(k = length ys)でマッチするなら、take k sw = ys かつ drop k sw ⊑ sx でなければなりません。この条件から、次の重要な等式が導けます。
llcp sw (drop k sw) = min i (m−k) (16.5)
(16.5) の意味
重要な観察
k ステップ分ずらしたときにマッチしうるための必要条件が (16.5) です。k は「sw を k ずらしても、まだ現在得ている情報 i と整合する」最小の値でなければならないのです。
証明の骨子は、次の 2 つの場合分けです。
- i < m−k のとき:take i (drop k sw) と take i sw が等しく、take (i+1) では等しくないことを追いかけると llcp sw (drop k sw) = i が出ます。
- i ≥ m−k のとき:drop k sw が sw の接頭辞となり、その長さは m−k。よって llcp sw (drop k sw) = m−k。
この 2 つを合わせれば min i (m−k) という形になります。
シフト量を式で書く
いま任意の i(0 ≤ i ≤ m)に対して、(16.5) を満たす最小の正の k(1 ≤ k ≤ m)を shift sw i と定義します。
shift sw i = head [k | k ← [1 .. m], llcp sw (drop k sw) == min i (m−k)]
Dart
// llcp: 最長共通接頭辞の長さ
int llcp<T>(List<T> a, List<T> b) {
int i = 0;
while (i < a.length && i < b.length && a[i] == b[i]) {
i++;
}
return i;
}
// shift sw i:(16.5) を満たす最小の正の k (1 <= k <= m)
int shift<T>(List<T> sw, int i) {
final int m = sw.length;
for (int k = 1; k <= m; k++) {
final int need = i < m - k ? i : m - k; // min i (m-k)
if (llcp(sw, sw.sublist(k)) == need) return k;
}
return m; // 理論上 k=m で必ず条件を満たす
}
m ≠ 0 のかぎり、より小さな k がなくても k = m は (16.5) を満たすので、必ず答えが見つかります。
まとめ
現在の窓で長さ i だけ一致したら、続く shift sw i − 1 個の窓は無視して構わない——マッチを見逃す心配はありません。
この式のままでは shift sw i の計算に最悪 Ω(m2) かかりますが、後の節で map (shift sw) [0 .. m] を O(m) で計算する方法を示します。
この時点の matches
matches ws = test · scanl step (0, [ ])
where
test [ ] = [ ]
test ((n, sx) : nxs) = if i == m
then n : test (drop (k−1) nxs)
else test (drop (k−1) nxs)
where i = llcp sw sx
k = shift sw i
(sw, m) = (reverse ws, length ws)
Dart
// シフトを使う中間版:一致長 i に応じて k-1 個の窓を飛ばす
List<int> matchesShift<T>(List<T> ws, List<T> xs) {
final List<T> sw = ws.reversed.toList();
final int m = sw.length;
// 全窓 [(n, sx), ...] を素直に scanl step (0, []) で作る
final List<(int, List<T>)> windows = [(0, <T>[])];
int n = 0;
List<T> sx = <T>[];
for (final T x in xs) {
n = n + 1;
sx = [x, ...sx];
windows.add((n, sx));
}
final List<int> result = [];
int idx = 0;
while (idx < windows.length) {
final (int wn, List<T> wsx) = windows[idx];
final int i = llcp(sw, wsx);
final int k = shift(sw, i);
if (i == m) result.add(wn);
idx += k; // 現在の窓 + (k-1) 個スキップ
}
return result;
}
ただし、この形と元の形が等価と言えるのは m ≠ 0 の場合に限られます。
最後の一手 — 前回の情報を持ち越す
次の窓で最初の k 文字だけを見ればよい場合
もう 1 つ改良の余地があります。i = llcp sw sx、k = shift sw i とし、さらに m−k ≤ i という条件が成り立つとします。このとき llcp (drop k sw) sx = m−k、つまり次の窓のうち、後半 m−k 文字はすでに一致することが確定しています。だから比較する必要はありません。
次の窓 (n+k, ys ++ sx)(length ys = k)について次のように計算できます。
llcp sw (ys ++ sx)
= {i' = llcp sw ys, i' ≤ k}
if i' == k then k + llcp (drop k sw) sx else i'
= {m−k ≤ i のとき llcp (drop k sw) sx = m−k}
if i' == k then m else i'
節約のしくみ
m−k ≤ i のときは、次の窓で先頭の k 文字だけを比較すれば十分です。m−k > i のときは節約できず、最大 m 回の比較が要ります。
最終プログラム
この節約を test に「次の窓のどこまで検査するか」を表すパラメータ j を加えて組み込むと、図 16.1 のプログラムが得られます(llcp と shift の中身以外は完成しています)。これは Galil (1979) 版の BM アルゴリズムです。shift の計算時間を除けば、長さ n のテキストに対する実行時間は O(m + n) です(証明は Gusfield (1997) の定理 3.2.3 参照)。
matches ws = test m · scanl step (0, [ ])
where
test j [ ] = [ ]
test j ((n, sx) : nxs) | i == m = n : test k (drop (k−1) nxs)
| m−k ≤ i = test k (drop (k−1) nxs)
| otherwise = test m (drop (k−1) nxs)
where i' = llcp sw (take j sx)
i = if i' == j then m else i'
k = shift sw i
(sw, m) = (reverse ws, length ws)
Dart
// 図 16.1:Galil 版 BM。j は「次窓のうち先頭何文字を検査するか」
// shift sw i の代わりに前計算した表 a[i] を参照する
List<int> matchesBM<T>(List<T> ws, List<T> xs) {
final List<T> sw = ws.reversed.toList();
final int m = sw.length;
if (m == 0) {
return List<int>.generate(xs.length + 1, (i) => i);
}
final List<int> a = buildShiftTable(sw); // a[i] == shift sw i
// すべての窓を作る
final List<(int, List<T>)> windows = [(0, <T>[])];
int n = 0;
List<T> sx = <T>[];
for (final T x in xs) {
n = n + 1;
sx = [x, ...sx];
windows.add((n, sx));
}
final List<int> result = [];
int j = m; // 次窓で検査する先頭長
int idx = 0;
while (idx < windows.length) {
final (int wn, List<T> wsx) = windows[idx];
final int iPrime = llcp(sw, wsx.sublist(0, j < wsx.length ? j : wsx.length));
final int i = (iPrime == j) ? m : iPrime;
final int k = a[i];
if (i == m) {
result.add(wn);
j = k;
} else if (m - k <= i) {
j = k;
} else {
j = m;
}
idx += k;
}
return result;
}
図 16.1 最終プログラム
シフト表を線形時間で作る
目標
shift sw i をそのまま毎回計算すると立方時間になってしまいます。そこで、shifts sw = map (shift sw) [0 .. m] をあらかじめ線形時間で計算して配列 a に格納し、shift sw i の呼び出しを a ! i に置き換えます。BM アルゴリズムでいちばん微妙な部分です。
2 種類の項に分解する
簡潔さのため f(k) = llcp sw (drop k sw) と書きます。f(m) = 0、f(k) ≤ m−k です。0 ≤ i ≤ m について次のように整理できます。
shift sw i
= {定義}
minimum [k | k ← [1 .. m], f(k) == min i (m−k)]
= {min の場合分け}
minimum ([k | k ← [1 .. m−i], f(k) == i] ++
[k | k ← [m−i+1 .. m], f(k)+k == m])
= {f(k) = i ⇒ k ≤ m−i}
minimum ([k | k ← [1 .. m], f(k) == i] ++
[k | k ← [m−i+1 .. m], f(k)+k == m])
2 つの内包表記の和になっています。これらを 1 つの配列にまとめて扱いたいのです。
accumArray を使う
ここで Haskell ライブラリ Data.Array の accumArray を持ち込みます(第 1 章で登場した関数です)。定義から次が直ちに従います。
(accumArray op e (0, m) vks) ! i = foldl op e [v | (v, k) ← vks, v == i]
ただし添字が範囲外だと未定義になるので、map fst vks ⊆ [0 .. m] が必要です。特に、
a = accumArray min m (0, m) vks
vks = [(f(k), k) | k ← [1 .. m]]
とすれば、
a ! i = minimum ([k | k ← [1 .. m], f(k) == i] ++ [m])
となり、これで最初の項が処理できました。
2 番目の項を混ぜ込む
次に、a の材料に vks' を追加します。
a = accumArray min m (0, m) (vks ++ vks')
ここで vks' は
[(i, minimum [k | k ← [m−i+1 .. m], f(k)+k == m]) | i ← [1 .. m]]
の並べ替えなら何でもよいです。すると shift sw i = a ! i となります。
vks' を右から作る
次の定義(リストを逆順に生成)が要件を満たします。
vks' = zip [m, m−1 .. 1] (foldr op [ ] vks)
vks = [(f(k), k) | k ← [1 .. m]]
op (v, k) ks = if v + k == m then k : ks else head ks : ks
Dart
// vks' を右から作る(foldr op [] vks を Dart で表現)
// xs[i] = f(k)+k==m を満たす k のうち i より大きい最小のもの
List<int> buildXs(List<(int, int)> vks, int m) {
List<int> ks = <int>[];
for (int idx = vks.length - 1; idx >= 0; idx--) {
final (int v, int k) = vks[idx];
if (v + k == m) {
ks = [k, ...ks];
} else {
ks = [ks.first, ...ks];
}
}
return ks;
}
op(f(m), m) [ ] = [m] となることに注意(f(m) = 0 なので)。たとえば xs = foldr op [ ] vks のとき次のような対応になります。
f 2 4 0 5 2 3 0 2 0
k 1 2 3 4 5 6 7 8 9
xs 4 4 4 4 6 6 9 9 9
読み方
xs の i 番目(0 始まり)は、「f(k)+k = m を満たす k のうち i より大きい最小のもの」です。vks' ではこの値が添字 m−i と組になるので、結局 i と xs !! (m−i) がペアになり、目的の値が入ります。
allcp と組み合わせて完成
前章の allcp を思い出しましょう。
allcp xs = [llcp xs (drop k xs) | k ← [0 .. length xs − 1]]
Dart
// allcp xs:xs と各 drop k xs の最長共通接頭辞長
// (簡潔さのため素朴実装。前章の O(n) アルゴリズムに差し替え可能)
List<int> allcp<T>(List<T> xs) {
return List<int>.generate(
xs.length,
(int k) => llcp(xs, xs.sublist(k)),
);
}
これは線形時間で計算できました。今回は先頭を落として末尾に llcp xs [ ](= 0)を加えた変種が必要なので、
allcp' xs = tail (allcp xs) ++ [0]
Dart
// allcp' xs:先頭を落として末尾に 0 を加えた変種
List<int> allcpPrime<T>(List<T> xs) {
final List<int> a = allcp(xs);
return [...a.sublist(1), 0];
}
と定義します。すると
[(f(k), k) | k ← [1 .. m]]
= {f の定義}
[(llcp sw (drop k sw), k) | k ← [1 .. m]]
= {zip の定義}
zip [llcp sw (drop k sw) | k ← [1 .. m]] [1 .. m]
= {allcp' の定義}
zip (allcp' sw) [1 .. m]
とまとまります。以上を組み合わせた最終形が次のとおりです。
a = accumArray min m (0, m) (vks ++ vks')
where
m = length sw
vks = zip (allcp' sw) [1 .. m]
vks' = zip [m, m−1 .. 1] (foldr op [ ] vks)
op (v, k) ks = if v + k == m then k : ks else head ks : ks
Dart
// shift 表 a を線形時間で構築する(a[i] == shift sw i)
// 初期値 m、演算 min、範囲 [0, m]、材料は vks ++ vks'
List<int> buildShiftTable<T>(List<T> sw) {
final int m = sw.length;
final List<int> a = List<int>.filled(m + 1, m);
// vks = zip (allcp' sw) [1..m]
final List<int> ap = allcpPrime(sw);
final List<(int, int)> vks = List<(int, int)>.generate(
m,
(int j) => (ap[j], j + 1),
);
// vks 由来:a[v] = min(a[v], k)
for (final (int v, int k) in vks) {
if (0 <= v && v <= m && k < a[v]) a[v] = k;
}
// vks' = zip [m, m-1, ..., 1] (foldr op [] vks)
final List<int> xs = buildXs(vks, m);
// i 番目(0 始まり)の添字は m - i、値は xs[i]
for (int i = 0; i < xs.length; i++) {
final int idx = m - i;
final int k = xs[i];
if (0 <= idx && idx <= m && k < a[idx]) a[idx] = k;
}
return a;
}
図 16.1 の shift sw i を a ! i に置き換えれば、matches 全体が線形時間で走ります。
Dart
void main() {
final List<String> ws = 'abcab'.split('');
final List<String> xs = 'ababcabcab'.split('');
print(matches(ws, xs)); // [7, 10] 仕様
print(matchesBasic(ws, xs)); // [7, 10] BM 基本形
print(matchesShift(ws, xs)); // [7, 10] シフト付き
print(matchesBM(ws, xs)); // [7, 10] 図 16.1(Galil 版)
// shift 表の中身を確認
final List<String> sw = ws.reversed.toList();
print(buildShiftTable(sw)); // 各 i に対する shift sw i
}
結び
BM アルゴリズムは Boyer と Moore (1977) が最初に述べたもので、Cormen 他 (2001)、Crochemore と Rytter (2003)、Gusfield (1997) にも詳しい議論があります。教科書ではしばしば不一致文字規則(bad character rule)と良い接尾辞規則(good suffix rule)という 2 つの規則で説明されますが、上の導出ではどちらも明示的には出てきません。
この章のねらい
- 基本形の導出は、式の形に応じて効率化の法則(scan の補題と foldl のタプリング則)を機械的にあてはめただけ。
- 「パターンを右から左に照合する」という BM の鍵となる着想は、endswith の自然な定義から自動的に浮かび上がる。
- その後の高速化(シフト量、次窓の再利用、シフト表の線形時間計算)は式の内容に踏み込む必要があるが、微妙な着想を含むアルゴリズムでは当然のこと。
参考文献
Boyer, R. S. and Moore, J. S. (1977). A fast string searching algorithm. Communications of the ACM 20, 762–72.
Cormen, T. H., Leiserson, C. E., Rivest, R. L. and Stein, C. (2001). Introduction to Algorithms, second edition. Cambridge, MA: The MIT Press.
Crochemore, M. and Rytter, W. (2003). Jewels of Stringology. Hong Kong: World Scientific.
Galil, Z. (1979). On improving the worst cast of the Boyer–Moore string matching algorithm. Communications of the ACM 22 (9), 505–8.
Gusfield, D. (1997). Algorithms on Strings, Trees and Sequences. Cambridge, UK: Cambridge University Press.
Lecroq, T. (2003). Experimental results on string matching algorithms. Software – Practice and Experience 25 (7), 727–65.