Boyer-Moore算法

Algorithm

August 01, 2020

author: minghu6

(2nd) revised at 2024-10-29

本章节内容需要以《前缀函数与KMP算法》作为前置章节。

之前的KMP算法将前缀匹配的信息用到了极致，

而 BM 算法背后的基本思想是通过 后缀匹配 获得比前缀匹配更多的信息来实现更快的跳转。

基础介绍

想象一下，如果我们的的模式字符串 $\texttt{pat}$ ，被放在文本字符串 $\texttt{string}$ 的头部，使它们的第一个字符对齐。

\[\begin{aligned} \qquad \texttt{pat}:\qquad &\texttt{EXAMPLE} \\ \texttt{string}:\qquad &\texttt{HERE IS A SIMPLE EXAMPLE} \dots \\ &\qquad\ \ \ \, \Uparrow \end{aligned}\]

$\texttt{pat}$ 的长度为 $\texttt{patlen}$ ，$\texttt{string}$ 的长度 $\texttt{stringlen}$。

按照算法惯例，设定字符串的基底为 $1$

追踪下 $\texttt{string}$ 上的第 $\texttt{patlen}$ 个字符 $\texttt{char}$ ，考虑我们能得到什么信息：

观察 1

如果我们知道 $\texttt{char}$ 不在 $\texttt{pat}$ 中, 可以直接将 $\texttt{pat}$ 和 $\texttt{string}$ 指针向下滑动 $\texttt{patlen}$ 个字符（然后重新开始后缀匹配）。

观察 2

如果我们知道 $\texttt{char}$ 在 $\texttt{pat}$ 中，并且离尾端至少有 $\Delta_1$ 个字符，可以直接将 $\texttt{pat}$ 和 $\texttt{string}$ 指针向下滑动 $\Delta_1$ 个字符。

以失配字符 $\texttt{char}$ 为参数，在 $\texttt{pat}$ 上构建计算 $\Delta_1$ 的函数，：

\[\begin{array}{ll} \textbf{int}\ \texttt{delta1} (\textbf{char}\ \texttt{char}) \\ \qquad \textbf{if}\ \texttt{char} \notin \texttt{pat}\\ \qquad\qquad\textbf{return}\ \texttt{patlen} \\ \qquad \textbf{else} \\ \qquad\qquad\textbf{return}\ \texttt{patlen} - i\ \ \text{/* } i = \max\lbrace i\vert\ \texttt{pat}[i] = \texttt{char} \rbrace \text{ */} \end{array}\]

观察 3-a

假设不是在 $\texttt{pat}$ 的最后一个字符，而是倒数第 $m+1$ 个字符上失配，那么 $\texttt{pat}$ 向下滑动 $k = \Delta_1 - m$ ，而 $\texttt{string}$ 指针仍然向下滑动 $k + m = \Delta_1$ （因为要重新从尾部开始匹配）。

观察 3-b

假设在倒数第 $m+1$ 个字符上失配，已匹配完的 $m$ 个字符的后缀子串为 $\texttt{subpat}$，作者把 $\texttt{pat}$ 上所有与后缀子串 $\texttt{subpat}$ 相等的子串称为 plausible reoccurrence，其中最右边的 $\texttt{subpat}$ 称为 rightmost plausible reoccurrence ，简称 $\texttt{rpr}$ 。

定义 $\texttt{rpr}(j)$ 表示后缀子串 $\texttt{pat}[j+1\dots]$ 的 $\texttt{rpr}$ 的起始位置，得到 $\Delta_2$ ：

\[\begin{array}{ll} \textbf{int}\ \texttt{delta2}(\textbf{int}\ j) \ \ \text{/* j为失配字符索引 */} \\ \qquad\textbf{return}\ \texttt{patlen} + 1 -\texttt{rpr}(j) \\ \end{array}\]

这样在失配时， $\texttt{pat}$ 向下滑动 $\max(\Delta_1,\Delta_2) - m$ ，$\texttt{string}$ 上的指针向下滑动 $\max(\Delta_1,\Delta_2)$ （全都需要重新从尾部开始匹配）。

实例说明

箭头指向失配字符 $\texttt{char}$：

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad\ \ \ \Uparrow \end{aligned}\]

$\texttt{F}$ 不在 $\texttt{pat}$ 中，根据 观察 1 ，把 $\texttt{pat}$ 直接向下滑动 $\texttt{patlen}$ 个字符，也就是 $7$ 个字符；

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\qquad\quad\ \ \, \texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad\ \ \, \ \qquad\quad\ \ \ \Uparrow \end{aligned}\]

根据 观察 2，我们需要将 $\texttt{pat}$ 向下滑动 $4$ 个字符使得 $\texttt{_}$ 字符对齐：

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\qquad\qquad\quad\ \ \, \texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad\ \ \, \qquad\qquad\qquad \Uparrow \end{aligned}\]

现在 $\texttt{char: T}$ 匹配了，把指针向前一步，继续匹配：

在 $\texttt{L}$ 上失配，因为 $\texttt{L}$ 不在 $\texttt{pat}$ 中，根据 观察 3-a ，$\texttt{pat}$ 向下滑动 $k= \Delta_1-m=7-1=6$ 个字符，而 $\texttt{string}$ 上指针向下移动 $\Delta_1=7$ 个字符：

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\qquad\qquad\qquad\qquad\ \ \,\, \texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad\ \ \ \; \qquad\qquad\qquad\qquad\ \ \ \Uparrow \end{aligned}\]

这时 $\texttt{char}$ 又一次匹配到了 $\texttt{pat}$ 的最后一个字符 $\texttt{T} $，继续向前匹配，匹配到了 $\texttt{A}$ ，再向前匹配，发现在字符 $\texttt{-}$ 失配：

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\qquad\qquad\qquad\qquad\ \ \,\, \texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad\ \, \ \qquad\qquad\qquad\quad\ \ \ \,\, \Uparrow \end{aligned}\]

此时根据 观察 3-b ，$\Delta_1=4$，$\Delta_2=7-0 = 7$，即 $\texttt{string}$ 上指针向下滑动 $\max(\Delta_1,\Delta_2)= 7$， $\texttt{pat}$ 向下滑动 $\max(\Delta_1,\Delta_2)-m= 7-2=5$，使得后缀 $\texttt{AT}$ 对齐：

\[\begin{aligned} \texttt{pat}:\qquad\qquad &\qquad\qquad\qquad\qquad\qquad\quad \;\, \texttt{AT-THAT} \\ \texttt{string}:\ \ \ \dots\ &\texttt{WHICH-FINALLY-HALTS.--AT-THAT-POINT} \dots \\ &\qquad \ \, \ \quad\qquad\qquad\qquad\qquad\quad \ \ \; \Uparrow \end{aligned}\]

现在我们在 $\texttt{string}$ 上找到了一个 $\texttt{pat}$ 的匹配。而我们只花费了$14$ 次对 $\texttt{string}$ 的引用，其中 $7$ 次是完成一个成功的匹配所必需的比较次数（ $\texttt{patlen}=7$ ），另外 $7$ 次让我们跳过了 $22$ 个字符！

算法设计

匹配算法

\[\begin{array}{ll} i \gets \texttt{patlen}. \\ \\ \textbf{loop}\\ \qquad \textbf{if}\ i > \texttt{stringlen} \\ \qquad \qquad \textbf{return}\ \texttt{false} \\ \\ \qquad j \gets \texttt{patlen} \\ \\ \qquad \textbf{loop} \\ \qquad\qquad \textbf{if}\ j < 1 \\ \qquad\qquad\qquad \textbf{return}\ i+1 \\ \\ \qquad\qquad\textbf{if}\ \texttt{string}[i]=\texttt{pat}[j] \\ \qquad\qquad\qquad j \gets j-1 \\ \qquad\qquad\qquad i \gets i-1 \\ \\ \qquad\qquad\qquad \textbf{continue} \\ \\ \qquad\qquad\textbf{break}\\ \\ \qquad i \gets i+\max(\texttt{delta1}(\texttt{string}[i]), \texttt{delta2}(j)) \\ \\ \end{array}\]

$i$ 表示 $\texttt{string}$ 上的全局指针，$j$ 表示 $\texttt{pat}$ 局部指针。

$\textbf{return}\ \texttt{false}$ 表示 $\texttt{pat}$ 不在 $\texttt{string}$ 中，返回数字表示 $\texttt{pat}$ 在 $\texttt{string}$ 左起第一次出现的位置。

$\texttt{rpr}(j)$

让我们完整定义下 $\texttt{rpr}(j)$ 函数。

根据前文定义，$\texttt{rpr}(j)$ 表示 $\texttt{subpat}=\texttt{pat}[j+1\dots]$ 在 $\texttt{pat}$ 上最右边相等子串的起始位置。

假设找到的 $\texttt{rpr}$ 字串为 $\texttt{pat}[k\dots k+\texttt{patlen}-j-1]=\texttt{subpat}$，首先允许拓展 $k$ 到 $k \lt 1$ ，相当于在 $\texttt{pat}$ 前面补充一段虚拟前缀，这就向下兼容了所有可能的 plausible reoccurrence 情况。

另外当 $k>1$ 时，考虑一种优化，如果 $\texttt{pat}[k]=\texttt{pat}[j]$ ，那么显然 $\texttt{pat}[k\dots k+\texttt{patlen}-j-1]$ 就是一个无效的 plausible reoccurrence ，因为 $pat[j]$ 本身是失配字符，那么向下滑动后重新从尾部开始匹配，在 $\texttt{pat}[k]$ 处还是会失配。

最后考虑到 $\texttt{delta2}(\texttt{patlen}) = 0$ ，所以设置 $\texttt{rpr}(\texttt{patlen}) = \texttt{patlen}$ 。

实例说明

\[\begin{aligned} \textit{j}:& &\texttt{1 2 3 4 5 6 7 8 9} \\ \texttt{pat}:& &\texttt{A B C X X X A B C} \\ \texttt{rpr}(j):& &\texttt{4 3 2 1 0 2 1 0 9} \\ \text{sign}:& &\texttt{- - - - - - - - +} \end{aligned}\]

$j = 1$ ，$\texttt{subpat = BCXXXABC}$， $\texttt{[(BCXXX)ABC]XXXABC}$，$\texttt{rpr}(1)=-4$；

$j = 6$ ，$\texttt{subpat = ABC}$，$\texttt{[ABC]XXXABC}$，$\texttt{rpr}(6)=1$；

$j = 7$ ，$\texttt{subpat = BC}$，$\texttt{A[BC]XXXABC}$，又因为 $\texttt{pat}[1]=\texttt{pat}[7]$，所以改为 $\texttt{[(BC)]ABCXXXABC}$ ，所以 $\texttt{rpr}(j)=-1$；

$j = 9 = \texttt{patlen}$，根据定义，$\texttt{rpr}(9)= \texttt{patlen}=9$。

以及另外一个例子：

\[\begin{aligned} \textit{j}:& &\texttt{1 2 3 4 5 6 7 8 9} \\ \texttt{pat}:& &\texttt{A B Y X C D E Y X} \\ \texttt{rpr}(j):& &\texttt{7 6 5 4 3 2 3 0 9} \\ \text{sign}:& &\texttt{- - - - - - + - +} \end{aligned}\]

匹配算法的改进

最后，实践过程中考虑到搜索过程中估计有 $80\%$ 的时间用在了基于 观察 1 和 观察 2 的跳转上，也就是 $\texttt{string}[i]$ 和 $\texttt{pat}[\texttt{patlen}]$ 不匹配，然后跳过整个 $\texttt{patlen}$ 的情况。

为此进行优化，定义一个 $\texttt{delta0}$：

\[\begin{array}{ll} \textbf{int}\ \texttt{delta0}(\textbf{char}\ \texttt{char}) \\ \qquad \textbf{if}\ \texttt{char}=\texttt{pat}[\texttt{patlen}] \\ \qquad\qquad \textbf{return}\ \texttt{large} \\ \qquad \textbf{return}\ \texttt{delta1}(\texttt{char}) \end{array}\]

其中 $\texttt{large}$ 是一个足够大的数，满足 $\texttt{large} \gt \texttt{stringlen + patlen}$ ，用作信号功能。

用 $\texttt{delta0}$ 代替 $\texttt{delta1}$ ，得到改进后的匹配算法：

\[\begin{array}{ll} i \gets \texttt{patlen} \\ \\ \textbf{if} \ i > \texttt{stringlen} \\ \qquad\textbf{return}\ \texttt{false}\\ \\ \textbf{loop} \\ \qquad\ \text{/* fast */}\\ \\ \qquad\ \textbf{while}\ i \leqslant \texttt{stringlen} \\ \qquad\qquad i \gets i + \texttt{delta0}(\texttt{string}[i])\\ \\ \qquad\ \text{/* undo */}\\ \\ \qquad\ \textbf{if}\ i \leqslant\ \text{large}\\ \qquad\qquad\textbf{return}\ \texttt{false}\\ \\ \qquad i \gets i-\texttt{large} - 1 \\ \qquad j \gets \texttt{patlen} - 1. \\ \\ \qquad\ \text{/* slow */}\\ \\ \qquad\textbf{loop}\\ \qquad\qquad \textbf{if}\ j < 1 \\ \qquad\qquad \qquad \textbf{return}\ i+1 \\ \\ \qquad\qquad\textbf{if}\ \texttt{string}[i]=\texttt{patlen}[j]\\ \qquad\qquad\qquad j \gets j-1 \\ \qquad\qquad\qquad i \gets i-1 \\ \\ \qquad\qquad\qquad\textbf{continue} \\ \\ \qquad\qquad\textbf{break}\\ \\ \qquad i \gets i + \max(\texttt{delta1}(\texttt{string}[i]), \texttt{delta2}(j)) \\ \\ \end{array}\]

$\text{fast：}$ 只要一直在 $\texttt{pat}$ 尾字符上失配，就不需要考虑 $\Delta_2$ 。

$\text{undo：}$ 从 $\text{fast}$ 滑落有两种情况，一种是正常的发现尾字符匹配，由于 $\texttt{+large}$ 导致滑落；

还有一种是 $\texttt{string}$ 指针走完了，于是滑落下来，对于此种情况 $0\leqslant i \leqslant \texttt{large} $ ，直接返回 $\texttt{false}$ ，表明子串不存在。

在继续下落前应该还原 $i$ 指针，并把它额外向前拨 $1$ ，因为尾字符已经匹配了。

$\text{slow:}$ 原算法的主要匹配代码。

经过改进，可以跳过 $\Delta_2$ 的多余计算，使得在通常字符集下搜索字符串的性能有明显的提升。

$\texttt{delta2}$

探寻 $\texttt{delta2}$ 实现，也是不易，要进行一番技术史的考据和对 KMP 算法融会贯通基础上进行重新发明。

历史细节

发表在 1977 年 10 月的 Communications of the ACM 上的在 Boyer、Moor 的论文¹只描述了一个静态表，并没有说明如何生成它。

对 $\texttt{delta2}$ 具体实现的讨论反而出现在 1977 年 6 月 Knuth、Morris、Pratt 在 SIAM Journal on Computing 上正式联合发表KMP算法的 Fast Pattern Matching in Strings² （除 KMP 外，还有若干字符串匹配的算法构想和介绍，其中就包括了 BM 算法），这听力来有点儿魔幻：

1969 年夏天 Morris 为某个大型机编写文本编辑器时利用有限自动机的理论发明了等价于 KMP 算法的字符串匹配算法，而他的算法由于过于复杂，被不理解他算法的同事当做 bug 修改得一团糟。
1970 年 KMP 中的 K，Kruth 在研究 Cook 的 two-way deterministic pushdown automaton 的理论时受到启发，也独立发明了 KMP 算法的雏形，并把它展示给他的同事 P，Pratt ，Pratt 改进了算法的数据结构。
1974 年 Boyer、Moor 发现通过更快地跳过不可能匹配的文本能实现比 KMP 更快的字符串匹配（ Gosper 也独立地发现了这一点），而一个只有原始 $\texttt{delta1}$ 定义的匹配算法是 BM 算法的最初版本。
1975 年 Boyer、Moor 提出了原始的 $\texttt{delta2}$ 表，而这个版本的 $\texttt{delta2}$ 表不仅不会对性能有所改善，还会在处理小字符表时拖累性能表现，而同年 MIT 人工智能实验室的 Kuipers 和 Knuth 向他们提出了类似的关于 $\texttt{delta2}$ 的改进建议，于是 Boyer、Moor 在论文的下一次修改中提到了这个建议，并设想用一个二维表代替 $\texttt{delta1}$ 和 $\texttt{delta2}$ 。
1976 年 1 月 Knuth 证明了有关 $\texttt{delta2}$ 的改进会得到更好的性能，于是 Boyer、Moor两人又一次修改了论文，得到了现在版本的 $\texttt{delta2}$ 定义。同年 4 月，斯坦福的 Floyd 又发现了 Boyer、Moor 两人第一版本的公式中的严重的统计错误，并给出了现在版本的公式。
Standish 又为 Boyer、Moor 提供了现在的匹配算法的改进。
1977 年 6 月 Knuth、Morris、Pratt 正式联合发表了 KMP 算法的论文，其中在提及比 KMP 表现更好的算法中提出了 $\texttt{delta2}$ 的构建方式。（其中也感谢了 Boyer、 Moor 对于证明线性定理（linearity theorem）提供的帮助）

这是一个团结、友谊和协作的故事，也是一个大佬带我飞的故事。

朴素算法

在介绍 Knuth 的 $\texttt{delta2}$ 构建算法之前，根据定义，我们会有一个 $O(n^3)$ 朴素算法：

虽然时间复杂度看起来糟糕透了，但考虑到 $\texttt{delta2}$ 的只构建一次而且 $\texttt{pat}$ 通常长度有限，

很多时候这已经够用了，除非是搜索长度成百上千的不同 $\texttt{pat}$ 。

对于每个 $j \in 1\dots \texttt{patlen}$ 在 $-\texttt{patlen}-j-1\dots j$ 上从右到左暴力寻找第一个等于 $\texttt{subpat}$ 的子串；
逐字符进行比较，当指针 $i\lt 1$ 时直接返回字符匹配成功（因为是虚拟前缀）；
当子串匹配成功时检查是否 $i\lt 1$ 或 $\texttt{pat}[i] \neq \texttt{pat}[j] $；
最后令 $\texttt{pat}[\texttt{patlen}] = 0 $ 。

注意实现代码里的字串的基底为 $0$

use std::cmp::Eq;

pub fn build_delta2_table_naive(p: &'a [impl Eq]) -> Vec<usize> {
    let patlen = p.len();
    let lastpos = patlen - 1;
    let mut delta2 = vec![];

    for j in 0..patlen {
        let subpatlen = (lastpos - j) as isize;

        if subpatlen == 0 {
            delta2.push(0);
            break;
        }

        for i in (-subpatlen..=j as isize).rev() {
            // subpat 匹配
            if (i..i + subpatlen)
            .zip(j + 1..patlen)
            .all(|(rpr_index, subpat_index)| {
                if rpr_index < 0 {
                    return true;
                }

                if p[rpr_index as usize] == p[subpat_index] {
                    return true;
                }

                false
            })
            && (i <= 0 || p[(i - 1) as usize] != p[j])
            {
                delta2.push((lastpos as isize - i) as usize);
                break;
            }
        }
    }

    delta2
}

特别地，对Rust语言特性进行必要地解释，下不赘述：

usize 和 isize 是和内存指针同字节数的无符号整数和有符号整数。
索引数组、向量、分片时使用 usize 类型的数字（因为在做内存上的随机访问并且下标不能为负值），所以如果需要处理负值要用isize，而进行索引时又要用usize，这就看到使用as关键字进行二者之间的数值转换。
impl PartialEq 用作泛型，这样可以同时支持Unicode编码的char和二进制的u8 和其他所有实现了 PartialEq 的数据类型。

这是个显见的时间复杂度为 $O(n^3)$ 的暴力算法。

高效算法

下面我们要介绍的是时间复杂度为 $O(n)$，但是需要额外 $O(n)$ 空间复杂度的高效算法。

需要指出的是，虽然 1977 年 Knuth 提出了这个构建方法，然而他的原始版本的构建算法存在一个缺陷，会产生不出符合定义的 $\texttt{delta2}$。

Rytter 在 1980 年 SIAM Journal on Computing 上发表的文章 ³ 对此提出了修正，但是 Rytter 这篇文章在细节上还有令人疑惑的地方，包括不限于

示例中奇怪的 $\Delta_2$ 数值（不清楚他依据的 $\texttt{delta2}$ 是否和最终版 $\texttt{delta2}$ 的定义有差别，但我实在不想因为这事儿继续考古了😱）
明显地复述 Knuth 算法时的笔误、算法上错误的缩进（可能是文章录入时的问题？）
可读性很差的变量命名（考虑到那是个 goto 语句、汇编语言和大型机的时代，随性的变量命名倒也很合理）

总之你绝不想看他那个修正算法的具体实现，不过好在他的文字描述比伪代码清晰多了，现在我们用更清晰的思路和代码结构去重新发明 $\texttt{delta2}$ 的构建算法：

$\texttt{delta2}$ 的定义比较复杂，高效实现的关键，是按照 $\texttt{subpat}$ 的 plausible reoccurrence（以下简称重现）位置分类处理。

按照重现位置从前到后，分为三类：

重现完全在 $\texttt{pat}$ 左边，比如 $\texttt{[(EYX)]ABYXCDEYX}$ ，此时 $\texttt{rpr}(j) = 1-\texttt{subpatlen}$ ，而根据定义 $\texttt{subpatlen}=\texttt{patlen} - j$ ，于是得到 $\texttt{delta2}(j) = \texttt{patlen} - \texttt{rpr}(j) + 1 = 2 \times \texttt{patlen} - j$；
重现有一部分在 $\texttt{pat}$ 左边，有一部分是 $\texttt{pat}$ 前缀，比如$\texttt{[(XX)ABC]XXXABC}$，此时 $\texttt{patlen} < \texttt{delta2}(j) < 2 \times \texttt{patlen} - j$；
重现完全在 $\texttt{pat}$ 中，比如 $\texttt{AB[YX]CDEYX} $ ，此时 $\texttt{delta2}(j) \leqslant \texttt{patlen} $。

接下来讨论如何高效地计算这三种情况：

case-1

这是最简单的情况， $\texttt{delta2}(j)$ 是确定值。

这也是重现的兜底情况，要通过一次遍历，在这里完成 $\texttt{delta2}$ 的初始化，之后看是否在 case-2 , case-3 上有更新的机会。

case-2

此时应该是 $\texttt{subpat}$ 的某个后缀和 $\texttt{pat}$ 的某个前缀相等。

比如之前的例子：

\[\begin{aligned} \textit{j}:& &\texttt{1 2 3 4 5 6 7 8 9} \\ \texttt{pat}:& &\texttt{A B C X X X A B C} \\ \end{aligned}\]

$\texttt{delta2}(4)$ 的重现 $\texttt{[(XX)ABC]XXXABC}$，在 $\texttt{subpat} = \texttt{XXABC}$ 的后缀与 $\texttt{pat}$ 的前缀中，有一个（最长的）相等子串 $\texttt{ABC}$。

这就来到了《前缀函数与KMP算法》所涉及的领域，实际上对第二和第三种情况的计算都离不开前缀函数，或者更一般地讲，DP 思想。

计算此种情况下的 $\texttt{delta2}(j)$：

设相等的前后缀长度为 $\texttt{prefixlen}$。

于是有 $\texttt{pat}$ 左边部分的长度为 $\texttt{subpatlen}-\texttt{prefixlen} = \texttt{patlen} - j - \texttt{prefixlen}$；

$\texttt{rpr}= \texttt{pat}[1-( \texttt{patlen} - j - \texttt{prefixlen})\dots \texttt{prefixlen}]$；

$\texttt{delta2}(j) = \texttt{patlen} - \texttt{rpr}(j) = 2 \times \texttt{patlen} - j - \texttt{prefixlen}$。

实际上要计算所有长度不同的 $\texttt{subpat}$ 后缀与 $\texttt{pat}$ 前缀相等的情况。

也就是计算 $\texttt{pat}$ 所有真（真就真在不包括它自己）后缀和真前缀相等的情况，然后按照长度从大到小，将 $\texttt{delta2}(j)$ 分区间进行计算。

假如 $\texttt{pat}$ 上有三组长度不同的相等前后缀，记为 $\texttt{pat}[i_1\dots]$、 $\texttt{pat}[i_2\dots]$ 和 $ \texttt{pat}[i_3\dots]$ 。

于是有：

$1\leqslant j \lt i_1 $ ， $\texttt{prefixlen} = \texttt{patlen} - i_1 + 1$
$i_1\leqslant j \lt i_2 $ ， $\texttt{prefixlen} = \texttt{patlen} - i_2 + 1$
$i_2\leqslant j \lt i_3 $ ， $\texttt{prefixlen} = \texttt{patlen} - i_3 + 1$
$i_3\leqslant j $ ， $\texttt{prefixlen} =0$ ，实际上是 case-1 的情况

如何计算 $\texttt{pat}$ 所有相等的真前缀和真后缀（的长度）呢？

可以利用 KMP 里的前缀函数。

让我们回顾一下前缀函数 $\pi[i]$ ，它表示 $\texttt{pat}[\dots i]$ 最长的那对儿相等真前缀和真后缀的长度。

$\pi[\texttt{patlen}]$ 就是 $\texttt{pat}$ 上最长的一对儿相等真前缀和真后缀。

那么如何计算其他次长、次次长、次次次长…… 的相等真前缀和真后缀长度呢？

回想一下前缀函数的计算过程，$\pi[i]$ 首先检查下 $\texttt{pat}[\pi[i-1] + 1]$ 是否等于 $\texttt{pat}[\texttt{patlen}]$ ；

否则再检查 $\texttt{pat}[\pi[\pi[i-1]] + 1]$ ，$\texttt{pat}[\pi[\pi[\pi[i-1]]] + 1]$ ，等等，直到等于 $\texttt{pat}[\texttt{patlen}]$ ，或者滑落到 $\texttt{pat}[1]$ 为止。

利用同样的规律，反推回去：

$\texttt{prefixlen}_1 = \pi[\texttt{patlen}]$；

$\texttt{prefixlen}_2 = \pi[\pi[\texttt{patlen}]]$；

$\texttt{prefixlen}_3 = \pi[\pi[\pi[\texttt{patlen}]]]$；

直到 $\texttt{prefixlen}_n = 0$ 。

这样就完成了对 $\texttt{delta2}(j)$ 的计算

case-3

重现完全在 $\texttt{pat}$ 中。

需要按照从右到左的顺序，在 $\texttt{pat}[0\dots \texttt{patlen}-1]$ 中搜索符合条件的 $\texttt{subpat}$ 。

开启脑洞：既然是个字符串搜索的问题，那么当然可以用 BM 算法本身递归地解决，设置结束条件是 $\texttt{patlen} = 1 $。

Knuth 发现如果从尾部开始建立“前缀”函数，失配时相等的子串恰好就是一个保证下一个字符和作为 $\texttt{pat}$ 后缀的 $\texttt{subpat}$ 的下一个字符不相同的（因为失配了）合法重现。

这个过程中只要检查每一个遇到的重现，就可以完成 case-3 的 $\texttt{delta2}$ 更新。

具体实现

pub fn compute_pi(pat: &[impl Eq]) -> Vec<usize> {
    let mut pi = vec![0usize; pat.len()];

    for i in 1..pat.len() {
        let mut j = pi[i - 1] as isize;

        while j > 0 && pat[i] != pat[j as usize] {
            j = pi[j as usize - 1] as isize;
        }

        if pat[i] == pat[j as usize] {
            j += 1;
        }

        pi[i] = j as usize;
    }

    pi
}

pub fn build_delta2_table_improved_minghu6(pat: &'a [impl Eq]) -> Vec<usize> {
    if pat.is_empty() {
        return Vec::new();
    }

    let patlen = pat.len();
    let lastpos = patlen - 1;
    let mut delta2 = Vec::with_capacity(patlen);

    /* case-1: delta2[j] = 2 * lastpos - j */

    for i in 0..patlen {
        delta2.push(lastpos * 2 - i);
    }

    /* case-2: lastpos <= delata2[j] < 2 * lastpos - j */

    let pi = compute_pi(pat);
    let mut prefixlen = pi[lastpos];
    let mut i = 0;

    while prefixlen > 0 {
        for j in i..(patlen - prefixlen) {
            delta2[j] = 2 * lastpos - j - prefixlen;
        }

        i = patlen - prefixlen;
        prefixlen = pi[prefixlen - 1];
    }

    /* case-3: delata2[j] < lastpos */

    let mut suffix = pi;
    suffix[lastpos] = patlen;

    let mut j = lastpos;

    for i in (0..patlen - 1).rev() {
        suffix[i] = j;

        while j < patlen && pat[i] != pat[j] {
            if delta2[j] > lastpos - i {
                delta2[j] = lastpos - i
            }

            j = suffix[j];
        }

        j -= 1;
    }

    delta2[lastpos] = 0;

    delta2
}

复杂度分析

case-1 初始化的时候是一个 $O(N)$ ；

case-2 计算前缀函数是一个 $O(N)$ ，和遍历 $\texttt{prefixlen}$ 不超过 $O(N)$；

case-3 从尾部计算前缀函数是一个 $O(N)$ ，回退过程顺带更新也不超过 $O(N)$；

加起来不超过 $O(5\cdot N)$ ，另外需要一个和 $\texttt{pat}$ 等长的辅助空间，用来存储前缀函数，

因此这个构建算法的时间复杂度是 $O(N)$ ，空间复杂度是 $O(N)$ 。

Galil规则

最坏情况

在不覆盖（onoverlapping）的多次匹配情况下，

如果 $\texttt{string}$ 的字母表很小（$\Delta_1$ 失效），并且 $\texttt{pat}$ 由非常短的子串循环构成（$\Delta_2$ 失效），最坏时间复杂度会回到 $O(m\cdot n) $ 。

一个极端的例子是 $\texttt{pat} = \texttt{AAA}$， $\texttt{string} = \texttt{AAAAA}\dots$。

最短周期 $k$

假定一个 $\texttt{pat}$ ，它是某个子串 $U$ 重复 $n$ 次构成的字符串 $UUUU\dots$ 的前缀，那么我们称 $U$ 为 $\texttt{pat}$ 的一个周期。

比如，对于 $\texttt{pat} = \texttt{ABCABCAB}$ ，$\texttt{ABC}$ 就是这个它的一个周期，当然 $\texttt{ABCABC}\dots$ 也是它的一个周期，但我们只关注最短的那个。

显然，$\texttt{pat}$ 至少拥有一个长度为它自身的周期，我们规定它的最短周期为 $k$，$k\leqslant \texttt{patlen}$。

在搜索过程中，假如我们的 $\texttt{pat}$ 成功地完成了一次匹配，那么依照周期的特点，实际上只需将 $\texttt{string}$ 上的指针向下滑 $k$ 个字符，比较这 $k$ 个字符是否等于 $\texttt{pat}[\texttt{patlen} - k + 1\dots ]$ 就可以直接判断是否存在 $\texttt{pat}$ 的又一个匹配。

而如何计算这个最短周期的长度呢，假如我们知道 $\texttt{pat}$ 的相等的一对儿前缀后缀，设它们的长度为 $\texttt{prefixlen}$，那么有 $\texttt{pat}[i] = \texttt{pat}[i+(patlen-\texttt{prefixlen})]$。

而从数学的角度看这个公式，显然我们已经有了长度为 $\texttt{patlen}-\texttt{prefixlen}$ 的周期，而当我们知道 $\texttt{pat}$ 最长的那一对相等的前缀后缀，我们就得到了 $\texttt{pat}$ 最短的周期。

而这个最长相等的前后缀长度，就是 $\pi[\texttt{patlen}]$ ，在我们在构建 $\texttt{delta2}$ 的时候已经计算过了。

实践评价

从实践的角度上说，理论上的最坏情况并不容易影响性能表现，哪怕是很小的只有 $4$ 的字符集的随机文本测试下这种最坏情况的影响也小到难以观察。

也因此如果没有很好地设计，使用 Galil 法则反而会拖累一点平均的性能表现，但对于一些极端特殊的比如上述例子中的 Gulil 规则的应用确实会使得性能表现提高数倍。

具体实现

字符类型

上世纪70年代那个时候，人们考虑字符，默认的前提是它是ASCII码，但现在的多字节编码方案不管是 GB18030 还是 Unicode ，采用的都是变长的字节编码方案，但是：

字符串匹配算法高效的关键在于字符索引的快速跳转
字符索引一定要建立在等宽字符的基础上，

所以匹配的对象要么是二进制字节，要么是固定码点的列表。

pub trait GetDelta1 {
    type Item;

    fn get_delta1(&self, char: &Self::Item) -> usize;
}


pub trait BMMatch: Sized {
    fn build_delta1_table<'a>(
        pat: &'a [Self],
    ) -> Box<dyn GetDelta1<Item = Self> + 'a>
    where
        Self: Eq + Hash,
    {
        let patlen = pat.len();
        let lastpos = patlen - 1;

        let mut map = HashMap::new();

        for (i, v) in pat.iter().enumerate() {
            map.insert(v, lastpos - i);
        }

        Box::new(HashMapWrapper { patlen, map })
    }
}


impl GetDelta1 for [usize; 256] {
    type Item = u8;

    fn get_delta1(&self, char: &Self::Item) -> usize {
        self[*char as usize] as _
    }
}


impl BMMatch for u8 {
    fn build_delta1_table<'a>(
        pat: &'a [Self],
    ) -> Box<dyn GetDelta1<Item = Self> + 'a> {
        let mut delta1 = [pat.len(); 256];
        let lastpos = pat.len() - 1;

        for i in 0..lastpos {
            delta1[pat[i] as usize] = lastpos - i;
        }

        Box::new(delta1)
    }
}


impl<'a> From<&'a str> for BMPattern<'a, u8> {
    fn from(pat: &'a str) -> Self {
        pat.as_bytes().into()
    }
}


impl<'a, M> From<&'a [M]> for BMPattern<'a, M>
where
    M: BMMatch + Eq + Hash,
{
    fn from(pat: &'a [M]) -> Self {
        let _delta1 = M::build_delta1_table(pat);
        let delta2 = build_delta2_table_improved_minghu6(pat);
        let k = compute_k(pat);

        BMPattern {
            pat,
            _delta1,
            delta2,
            k,
        }
    }
}


impl<'a, M: BMMatch + Eq> BMPattern<'a, M> {
    fn delta0(&self, char: &M) -> usize {
        if char == self.pat.last().unwrap() {
            LARGE
        }
        else {
            self.delta1(char)
        }
    }

    fn delta1(&self, char: &M) -> usize {
        self._delta1.get_delta1(char)
    }

    fn find(&self, string: &[M]) -> Option<usize> {
        if self.pat.is_empty() {
            return Some(0);
        }

        if string.is_empty() {
            return None;
        }

        let pat = self.pat;
        let patlastpos = pat.len() - 1;
        let stringlastpos = string.len() - 1;

        let mut i = patlastpos;

        if patlastpos == 0 {
            return string.iter().enumerate().find_map(|(i, v)| {
                if *v == self.pat[0] {
                    Some(i)
                }
                else {
                    None
                }
            });
        }

        if i > stringlastpos {
            return None;
        }

        loop {
            while i <= stringlastpos {
                i += self.delta0(&string[i]);
            }

            if i < LARGE {
                return None;
            }

            i -= LARGE + 1;

            let mut j = patlastpos - 1;

            loop {
                if string[i] == pat[j] {
                    if j == 0 {
                        return Some(i);
                    }

                    j -= 1;
                    i -= 1;

                    continue;
                }

                break;
            }

            i += max(self.delta1(&string[i]), self.delta2[j]);
        }
    }

    /// Find overlapping
    fn find_all(&'a self, string: &'a [M]) -> impl Iterator<Item = usize> + 'a
    where
        M: Debug,
    {
        std::iter::from_coroutine(
            #[coroutine]
            move || {
                let patlen = self.pat.len();
                let k = self.k;

                let mut suffix = string;
                let mut base = 0;

                while let Some(mut i) = self.find(suffix) {
                    yield i + base;

                    loop {
                        if suffix.len() < i + patlen + k {
                            return;
                        }

                        if suffix[i + patlen..i + patlen + k] == self.pat[patlen - k..] {
                            i += k;

                            yield i + base;
                        }
                        else {
                            break;
                        }
                    }

                    let shift = i + 1;

                    base += shift;

                    suffix = &suffix[shift..];
                }
            },
        )
    }
}


/// Galil rule shortest period k
///
/// k = patlen - prefixlen
pub fn compute_k(p: &[impl Eq]) -> usize {
    let patlen = p.len();
    let lastpos = patlen - 1;

    let pi = compute_pi(p);

    patlen - pi[lastpos]
}

算法变种

也可以叫做只使用 $\Delta_1$ 的算法变种。

Simplified Boyer-Moore 算法

BM 算法最复杂的地方就在于构建 $\texttt{delta2}$ 表（有一个通俗的名字，好后缀 表）的构建，而实践中一般字符集上的匹配性能只在于 $\Delta_1$ （通俗的名字是 坏字符 表）（后面的理论计算也证明了这一点），于是出现了仅仅使用 $\texttt{delta1}$ 表的简化版BM 算法。

Boyer-Moore-Horspol 算法

Horspol 算法同样是基于坏字符的规则，不过是在与 $\texttt{pat}$ 尾部对齐的字符上应用 $\Delta_1$ ，这个效果类似于前文对匹配算法的改进，所以它的通常表现优于原始 BM、和匹配算法改进后的 BM 差不多。

impl<'a, M> From<&'a [M]> for HorspoolPattern<'a, M>
where
    M: BMMatch + Eq + Hash,
{
    fn from(pat: &'a [M]) -> Self {
        let _delta1 = M::build_delta1_table(pat);

        Self { pat, _delta1 }
    }
}


impl<'a, M: Eq> HorspoolPattern<'a, M> {
    fn delta1(&self, char: &M) -> usize {
        self._delta1.get_delta1(char)
    }

    pub fn find(&self, string: &'a [M]) -> Option<usize> {
        if self.pat.is_empty() {
            return Some(0);
        }

        let stringlen = string.len();
        let patlastpos = self.pat.len() - 1;

        let mut i = patlastpos;

        while i < stringlen {
            if &string[i - patlastpos..=i] == self.pat {
                return Some(i - patlastpos);
            }

            i += self.delta1(&string[i]);
        }

        None
    }
}

Boyer-Moore-Sunday 算法

Sunday 算法同样是利用坏字符规则，只不过相比 Horspool 更进一步，属于是不演了：

既然绝大多数情况下匹配算法做的工作就是在跳适配字符，那么匹配失败后直接关注 $\texttt{pat}$ 尾部对齐的那个字符的下一个字符（因为下一个字符一定也在比较的序列里），从而可以更快地失败。

需要稍微修改一下 $\texttt{delta1}$ 表，偏移量 $+1$ ，并且将计算尾字符也计算在内。

Sunday 算法通常用作一般情况下实现最简单而且平均表现最好之一的实用算法，通常表现比 Horspool、BM 都要快一点。

pub struct SundayPattern<'a, M> {
    pat: &'a [M],
    _delta1: Box<dyn GetDelta1<Item = M> + 'a>,
}


impl<'a, M> From<&'a [M]> for SundayPattern<'a, M>
where
    M: BMMatch + Eq + Hash,
{
    fn from(pat: &'a [M]) -> Self {
        let _delta1 = M::build_delta1_table(pat);

        Self { pat, _delta1 }
    }
}

impl<'a, M: Eq> SundayPattern<'a, M> {
    fn delta1(&self, char: &M) -> usize {
        if char == self.pat.last().unwrap() {
            1
        }
        else {
            self._delta1.get_delta1(char) + 1
        }
    }

    pub fn find(&self, string: &'a [M]) -> Option<usize> {
        if self.pat.is_empty() {
            return Some(0);
        }

        let stringlen = string.len();
        let patlastpos = self.pat.len() - 1;

        let mut i = patlastpos;

        while i < stringlen {
            if &string[i - patlastpos..=i] == self.pat {
                return Some(i - patlastpos);
            }

            if i + 1 == stringlen {
                break;
            }

            i += self.delta1(&string[i + 1]);
        }

        None
    }
}

BMHBNFS 算法

是 CPython 实现 stringlib 模块时用到的 find 的算法⁴，似乎国内更有名气，不清楚为何叫这个名字，怎么就“AKA”了，以下简称 B5S。

既要快

该算法属于是还有高手，结合了 Horspool 和 Sunday ，不仅追求快速失败，还要尽量一口气跳过整个 $\texttt{pat}$ ：

先看尾端的下一个字符，如果不能靠这个字符跳过整个 $\texttt{pat}$ ，就再试下尾端的字符。

又要省

对于 Horspool 只保留尾字符失配时的 $\Delta_1$ 叫做 $\texttt{skip}$ ；

对于 Sunday 只用一个 $\texttt{pat}$ 上构建的字母表，字母表还要用一个 概率数据结构 来尽可能节省空间，当末尾的下一个字符不在字母表里时直接跳过 $\texttt{patlen} + 1$ 。

Bloom 过滤器

Bloom 过滤器是一种 概率数据结构 ，是通过牺牲准确率（和运行时间）来节省存储空间的 Set 类型的数据结构，存在假阳性（False Positive，FP），也就是集合中不存在的项会被误判为存在，但不存在假阴性（False Negatives，FN），也就是不会把存在的项误判为不存在。

这就像一个筛子一样，把一些不符合条件的值过滤出去。

一个标准的 Bloom 过滤器，要达到最佳 FP 概率，需要使用多个哈希函数做多次检查，但由此带来的时间常数是不可接受的。

这里简单地通过一个掩码把一部分二进制位映射到数据里，而目前已知最快的非加密哈希算法 xxHash⁵，所需时间都要高一个数量级。

考虑做一个 $8$ bit 数据的字母表，一般情况需要 $2^8 = 256$ 个 bool 空间，用 bit 代替 byte 做存储，可以省掉 $3$ 个信息位，如果再舍弃其中 $2$ 个信息位，那么就只需要 $2^{8-5} = 2^3 = 8$ 个字节也就是一个 u64 。

#[repr(transparent)]
pub struct SimpleBloomFilter {
    data: u64,
}


impl SimpleBloomFilter {
    pub fn new() -> Self {
        SimpleBloomFilter { data: 0 }
    }

    #[inline]
    pub fn insert(&mut self, elem: &u8) {
        (self.data) |= 1u64 << (elem & 63);
    }

    #[inline]
    pub fn contains(&self, elem: &u8) -> bool {
        (self.data & (1u64 << (elem & 63))) != 0
    }
}

这里做得实际就是将一个字节的信息映射为它的前 $6$ 位，原作者建议至少要映射 $5$ 位。

具体实现

pub struct B5SSpacePattern<'a, M> {
    pat: &'a [M],
    patalphabet: SimpleBloomFilter,
    /// patpatlastpos delta1
    skip: usize,
}


impl<'a> B5SSpacePattern<'a, u8> {
    fn build(pat: &'a [u8]) -> (SimpleBloomFilter, usize) {
        let mut patalphabet = SimpleBloomFilter::new();
        //let mut alphabet = FastBloomFilter::with_rate(p.len(), 0.15);
        let lastpos = pat.len() - 1;
        let mut skip = pat.len();

        for i in 0..pat.len() - 1 {
            patalphabet.insert(&pat[i]);

            if pat[i] == pat[lastpos] {
                skip = lastpos - i;
            }
        }

        patalphabet.insert(&pat[lastpos]);

        (patalphabet, skip)
    }

    pub fn find(&self, string: &'a [u8]) -> Option<usize> {
        if self.pat.is_empty() {
            return Some(0);
        }

        let stringlen = string.len();
        let patlen = self.pat.len();
        let patlastpos = patlen - 1;

        let mut i = patlastpos;

        while i < stringlen {
            if string[i] == self.pat[patlastpos] {
                if string[i-patlastpos..i] == self.pat[..patlastpos] {
                    return Some(i - patlastpos);
                }

                if i + 1 == stringlen {
                    break;
                }

                if !self.patalphabet.contains(&string[i + 1]) {
                    // sunday
                    i += patlen + 1;
                }
                else {
                    // horspool
                    i += self.skip;
                }
            }
            else {
                if i + 1 == stringlen {
                    break;
                }

                if !self.patalphabet.contains(&string[i + 1]) {
                    // sunday
                    i += patlen + 1;
                }
                else {
                    i += 1;
                }
            }
        }

        None
    }
}

这个算法在通常情况下，在 $O(1)$ 的空间复杂度下却有着最快的时间性能，简直神中神！

只是在特殊的小字符集上的性能表现会回退到 $O(n\cdot m)$，而这是原作者所能接受的。

理论分析

现在我们通过一个简单的概率模型来做一些绝不枯燥的理论分析，借此可以发现一些有趣而更深入的事实。

建立模型

想象一下，我们滑动字符串 $\texttt{pat}$ 到某个新的位置，这个位置还没有完成匹配，我们可以用发现失配所需要的代价与失配后 $\texttt{pat}$ 能够向下滑动的距离的比值来衡量算法的平均性能表现。

假如这个代价是用对 $\texttt{string}$ 的引用来衡量，那么我们就可以知道平均每个字符需要多少次 $\texttt{string}$ 的引用，这是在理论上衡量算法表现的关键指标；

而如果这个代价是用机器指令衡量，那我们可以知道平均每个字符需要多少条机器指令；

当然也可以有其他的衡量方式，这并不影响什么，这里我们采用对 $\texttt{string}$ 的引用进行理论分析。

同时为我们的概率模型提出一个假设：

$\texttt{pat}$，$\texttt{string}$ 中的每个字符是独立随机变量，它们出现的概率相等，为 $p$ ，$p$ 取决于全字母表的大小。

显然，假如全字母表的大小为 $q$ ，则 $p=\dfrac{1}{q}$ ，例如我们基于字节的串匹配算法，可以近似为 $q=\dfrac{1}{256}$。

现在可以更准确地刻画这个比率，$\texttt{rate}(\texttt{patlen}, p)$：

\[\frac{\sum_{m=0}^{\texttt{patlen}-1} \texttt{cost}(m) \times \texttt{prob}(m)}{\sum_{m=0}^{\texttt{patlen}-1} \texttt{prob}(m) \times \sum_{k=1}^{\texttt{patlen}} \texttt{skip}(m,k) \times k }\]

其中，$\texttt{cost}(m)$ 为前面讨论到的在匹配成功了 $m$ 个字符后失配时的代价：

\[\texttt{cost}(m) = m+1\]

$\texttt{prob}(m)$ 为匹配成功 $m$ 个字符后失配的概率（其中 $1-p^{\texttt{patlen}}$ 排除掉 $\texttt{pat}$ 全部匹配的情况）：

\[\texttt{prob}(m) = \frac{p^m(1-p)}{1-p^{\texttt{patlen}}}\]

$\texttt{skip}(m,k)$ 为发生失配时 $\texttt{pat}$ 向下滑动 $k$ 个字符的概率。

实际上所有字符串匹配算法的核心就是 $\texttt{skip}(m,k)$，下面我们会通过分析 $\Delta_1$ 和 $\Delta_2$ 来计算 BoyerMoore 算法的$\texttt{skip}(m,k)$。

计算 $\texttt{skip}(m,k)$

$\Delta_1$

首先考虑 $\Delta_1$ 不起作用的情况，也就是发现失配字符在 $\texttt{pat}$ 上重现的位置在已经匹配完的 $m$ 个字符中，这种情况的概率$\texttt{probdelta1-worthless}$为：

\[\texttt{probdelta1-worthless}(m) = 1 - (1-p)^m\]

而对于 $\Delta_1$ 起作用的情况，可以根据 $k$ 的范围分为四种情况进行讨论：

$k = 1$

1.1 失配字符对应位置的下一个字符恰好等于失配字符；

1.2 失配字符已经是 $\texttt{pat}$ 右手起最后一个字符。
$1 \lt k \lt \texttt{patlen}-m$，$\texttt{pat}$ 在失配字符对应位置的左边还有与失配字符相等的字符，并且不满足 1.；
$k = \texttt{patlen} - m$，在失配字符对应位置左边找不到另一个与失配字符相等的字符，并且不满足 1. ，这时 $\texttt{pat}$ 有最大的向下滑动距离；
$k \gt \texttt{patlen} - m$，对于 $\Delta_1$，这是不可能存在的情况。

于是有计算 $\Delta_1$ 的概率函数：

\[probdelta1(m,k) = \left\{ \begin{array}{lcl} (1-p)^m\times \left\{ \begin{array}{lcl} 1 & \text{for} &m+1=\texttt{patlen} \\ p &\text{for}&m+1\neq \texttt{patlen} \\\end{array}\right. & \text{for} & k = 1 \\ (1-p)^{m+k-1}\times p & \text{for} & 1 < k < \texttt{patlen} - m\\ (1-p)^{\texttt{patlen}-1} & \text{for} & k = \texttt{patlen} - m\\ 0 & \text{for} & \texttt{patlen} - m < k \leqslant \texttt{patlen} \end{array}\right.\]

$\Delta_2$

对于 $\Delta_2$ 概率的计算，根据定义，首先计算某个 $\texttt{subpat}$ 的重现概率，只要考虑该重现左边还有没有字符来提供额外的判断与失配字符是否相等的检查：

\[\texttt{probpr}(m,k) = \left\{ \begin{array}{lcl} (1-p)\times p^m & \text{for} & 1 \leqslant k < \texttt{patlen} - m\\ p^{\texttt{patlen}-k} & \text{for} & \texttt{patlen} - m \leqslant k \leqslant \texttt{patlen} \end{array}\right.\]

于是 $\texttt{delta2}(m,k)$ 就可以通过保证 $\texttt{pr}(m,k)$ 存在并且 $k$ 更小的 $\Delta_2$ 不存在，来递归计算：

\[\texttt{probdelta2}(m,k) = \texttt{probpr}(m,k)(1-\sum_{n=1}^{k-1} \texttt{probdelta2}(m, n))\]

$\Delta_1 \cup \Delta_2 - \Delta_1 \cap \Delta_2$

前面已经独立讨论了 $\Delta_1$ ，$\Delta_2$ 的概率函数，不过还需要额外考虑一下这两个概率函数之间相互影响的情况，虽然只是一个很少数的情况：

当 $\texttt{delta2}$ 计算的 $k$ 为 $1$ 的时候，根据 $\Delta_2$ 定义我们就知道

$\texttt{pat}[\texttt{patlen}-m] = \texttt{pat}[\texttt{patlen}-m + 1]\dots =\texttt{pat}[\texttt{patlen}]$

而这种情况已经排除了 $\Delta_1$ 不起作用的情况，因为当如前文讨论的，$\Delta_1$ 不起作用要求与失配字符 $\texttt{pat}[\texttt{patlen}-m]$ 相等的字符出现在 $\texttt{pat}[\texttt{patlen}-m+1]\dots \texttt{pat}[\texttt{patlen}]$ 中，这就产生了不可能在倒数 $m+1$ 个字符上失配的矛盾。

因此针对 $\Delta_1$ 不起作用的情况需要一个稍微修改过的 $\Delta_2$ 概率函数：

\[\texttt{probdelta2'}(m,k) = \left\{ \begin{array}{lcl} 0 & \text{for} & k = 1\\ \texttt{probpr}(m,k)(1-\sum_{n=2}^{k-1} \texttt{probdelta2}'(m, n)) & \text{for} & 1 \leqslant k \leqslant \texttt{patlen} \end{array}\right.\]

于是通过组合 $\Delta_1$ 和 $\Delta_2$ 起作用的情况，我们就得到了 BoyerMoore 算法的 $\texttt{skip}$ 概率函数：

\[\texttt{skip}(m,k) = \left\{\begin{array}{lcl} \texttt{probdelta1}(m, 1) \times \texttt{probdelta2}(m, 1) & & k = 1\\ \\ \\ \\ \texttt{prodelta1 - worthless}(m)\times \texttt{probdelta2'}(m, 1)\\ \quad +\ \texttt{probdelta1}(m, k)\times \sum_{n=1}^{k-1} \texttt{probdelta2}(m, n)\\ \quad +\ \texttt{probdelta2}(m, k)\times \sum_{n=1}^{k-1} \texttt{probdelta1}(m,n)\\ \quad +\ \texttt{probdelta1}(m, k)\times \texttt{probdelta2}(m, k) & & 1 < k \leqslant \texttt{patlen} \end{array}\right.\]

分析比较

为了结构清晰、书写简单、演示方便，我们使用 Python 平台的 Lisp 方言 Hy (v0.13) 来进行实际计算：

myprob.hy

(require [hy.contrib.sequences [defseq seq]])

(import [hy.contrib.sequences [Sequence end-sequence]])
(import [hy.models [HySymbol]])


(defmacro simplify-probfn [patlen p probfn-list]
    "(prob-xxx patlen p m k) -> (prob-xxx-s m k)"
    (lfor probfn probfn-list
        [(setv simplified-probfn-symbol (HySymbol (.format "{}-s" (name probfn))))
        `(defn ~simplified-probfn-symbol [&rest args] (~probfn patlen p #*args))]))

(defn map-sum [range-args func]
    (setv [start end] range-args)
    (-> func (map (range start (inc end))) sum))

(defn cost [m] (+ m 1))

(defn prob-m [patlen p m]
    (/
        (* (** p m) (- 1 p))
        (- 1 (** p patlen))))

(defn prob-delta_1 [patlen p m k]
    (cond [(= 1 k) (* (** (- 1 p) m)
                      (if (= (inc m) patlen) 1 p))]
          [(< k (- patlen m)) (* p (** (- 1 p) (dec (+ k m))))]
          [(= k (- patlen m)) (** (- 1 p) (dec patlen))]
          [(> k (- patlen m)) 0]))


(defn prob-delta_1-worthless [p m] (- 1 (** (- 1 p) m)))


(defn prob-pr [patlen p m k] (if (< patlen (+ m k))
                                (* (- 1 p) (** p m))
                                (** p (- patlen k))))


(defn prob-delta_2 [patlen p m]
    "prob-delta_2(_, k) = prob-delta_2-seq[k]"
    (defseq prob-delta_2-seq [n]
        (cond [(< n 1) 0]
              [(= n 1) (prob-pr patlen p m 1)]
              [(> n 1) (* (prob-pr patlen p m n) (- 1 (sum (take n prob-delta_2-seq))))]))
    prob-delta_2-seq)


(defn prob-delta_2-1 [patlen p m]
    (defseq prob-delta_2-1-seq [n]
        (cond [(< n 2) 0]
              [(>= n 2) (* (prob-pr patlen p m n) (- 1 (sum (take n prob-delta_2-1-seq))))]))
    prob-delta_2-1-seq)


(defn skip [patlen p m k prob-delta_2-seq prob-delta_2-1-seq]
    (simplify-probfn patlen p [prob-delta_1])
    (if (= k 1) (* (prob-delta_1-s m 1) (get prob-delta_2-seq 1))
        (sum [(* (prob-delta_1-worthless p m) (get prob-delta_2-1-seq k))
              (* (prob-delta_1-s m k) (sum (take k prob-delta_2-seq)))
              (* (get prob-delta_2-seq k) (map-sum [1 (- k 1)]
                                                  (fn [n] (prob-delta_1-s m n))))
              (* (prob-delta_1-s m k) (get prob-delta_2-seq k))])))


(defn bm-rate [patlen p]
    (simplify-probfn patlen p [prob-m prob-delta_2 prob-delta_2-1 skip])
    (/
        (map-sum [0 (dec patlen)]
            (fn [m] (* (cost m) (prob-m-s m))))

        (map-sum [0 (dec patlen)]
            (fn [m]
                (setv prob-delta_2-seq (prob-delta_2-s m)
                      prob-delta_2-1-seq (prob-delta_2-1-s m))
                (* (prob-m-s m) (map-sum [1 patlen]
                                    (fn [k] (* k (skip-s m k prob-delta_2-seq prob-delta_2-1-seq)))))))))

并且为了进行比较，还额外计算了简化BM算法:

myprob.hy

(defn simplified-bm-skip [patlen p m k]
    (simplify-probfn patlen p [prob-delta_1])
    (if (= k 1) (+ (prob-delta_1-s m 1) (prob-delta_1-worthless p m))
        (prob-delta_1-s m k)))

(defn sbm-rate [patlen p]
    (simplify-probfn patlen p [prob-m simplified-bm-skip])
    (/
        (map-sum [0 (dec patlen)]
            (fn [m] (* (cost m) (prob-m-s m))))

        (map-sum [0 (dec patlen)]
            (fn [m] (* (prob-m-s m) (map-sum [1 patlen]
                                        (fn [k] (* k (simplified-bm-skip-s m k)))))))))

和 KMP 算法：

myprob,hy

(defn prob-pi [patlen p]
    "prob-pi-s(m, l) = prob-pi-seq[m][l]"
    (defseq prob-pi-seq [n]
        (cond [(= n 0) []]
              [(= n 1) [1]]
              [(= n 2) [(- 1 p) p]]
              [(> n 2) (lfor i (range n)
                            (+
                                (* (getone prob-pi-seq (dec n) i :default 0) (- 1 p))
                                (* (get prob-pi-seq (dec n) (dec i)))))]))
    prob-pi-seq)


(defn skip [m k prob-pi-seq]
    (if (= m 0) (if (= k 1) 1 0)
        (get prob-pi-seq m (- m k))))

(defn at-least-1 [n]
    (if (< n 1) 1 n))


(defn kmp-rate [patlen p]
    (simplify-probfn patlen p [prob-m prob-pi])
    (setv prob-pi-seq (prob-pi-s))
    (/
        (map-sum [0 (dec patlen)]
            (fn [m] (* (cost m) (prob-m-s m))))

        (map-sum [0 (dec patlen)]
            (fn [m] (* (prob-m-s m) (map-sum [1 (at-least-1 (dec m))]
                                        (fn [k] (* k (skip m k prob-pi-seq)))))))))

然后我们就可以通过Python上的plotnine图形包看一下计算的数据（并用高斯过程回归拟合曲线）：

import pandas as pd
from plotnine import *
import hy

from my_prob import bm_rate, sbm_rate, kmp_rate


theme_update(text=element_text(family="SimHei"))

def plot(p, title, N=30):
    model_range = range(1, N+1)
    data = {'rate':[], 'alg':[], 'patlen':[]}
    categories_list = [(bm_rate, 'BoyerMoore'),
                       (sbm_rate, 'S BoyerMoore'),
                       (kmp_rate, 'KMP'),
                       (lambda patlen, p: 1/patlen, '$\\frac{1}{patlen}$')]

    for alg_fun, label in categories_list:
        data['rate'].extend([alg_fun(\texttt{patlen}, p) for patlen in model_range])
        data['alg'].extend([label for _ in model_range])
        data['patlen'].extend(model_range)

    df = pd.DataFrame(data)

    return (ggplot(df, aes(x='patlen', y='rate', color='alg'))
            + geom_point()
            + geom_smooth(method='gpr')
            + labs(color='Algs', title=title, x='$patlen$', y='$\\frac{cost}{skip}$'))

plot(1/256, '$p= \\frac{1}{256}$')：

观察这个图像，令人印象深刻的首先就是抬头的一条天蓝色线，几乎笔直地画出了算法性能的下限，不愧是 KMP 算法，$O(n)$ 的时间复杂度，一看就很真实。( ﾟ▽ﾟ)/

接着会发现 BoyerMoore 算法与简化版 BoyerMoore 算法高度重叠的这条红绿紫曲线，同时也是 $\dfrac{1}{\texttt{patlen}}$ ，

这就是在一般字符集下随机文本搜索能达到的$O(\dfrac{n}{m})$的强力算法吗？ (ﾟ△ﾟ;ﾉ)ﾉ

另外此时可以绝大多数的字符跳转依靠 $\Delta_1$（比 $\Delta_2$ 高几个数量），这也是基于 $\texttt{delta1}$ 表的 BM 变种算法最佳的应用场景！

接着我们可以看一下在经典的小字符集，比如在 DNA { A, C, T, G } 碱基对序列中算法的性能表现（plot(1/4, '$p= \\frac{1}{4}$')）：

曲线出现了明显的分化，当然 KMP 还是一如既往地稳定。 ( ﾟ▽ﾟ)/

如果此时在测试中监控一下 $\texttt{delta1}$ 表和 $\texttt{delta2}$ 表作用情况，会发现 $\texttt{delta2}$ 起作用的次数超过了 $\texttt{delta1}$ ，而且贡献的跳过字符数更是远超 $\texttt{delta1}$ ，思考下，这件事其实也很好理解。

总结一下，通过概率模型的计算，一方面看到了在较大的字符集，比如日常搜索的过程中 BoyerMoore 系列算法的优越表现，其中主要依赖 $\texttt{delta1}$ 表实现字符跳转；另一方面，在较小的字符集里，$\texttt{delta1}$ 的作用下降，而 $\texttt{delta2}$ 的作用得到体现。

如果有一定富裕空间的情况下，使用完整的空间复杂度为 $O(m)$ 的 BoyerMoore 算法应该是一种适用各种情况、综合表现都很优异的算法选择。

引用

https://dl.acm.org/doi/10.1145/359842.359859 ↩
https://epubs.siam.org/doi/abs/10.1137/0206024 ↩
https://epubs.siam.org/doi/10.1137/0209037 ↩
https://web.archive.org/web/20040907193754/http://effbot.org/zone/stringlib.htm ↩
https://cyan4973.github.io/xxHash/ ↩