String Searching Algorithms

表 1 是后文中所介绍字符串匹配算法的时间复杂度。

说明 1：表 1 中 \(n\) 为文本串长度， \(m\) 为模式串长度， \(\mid \sum \mid\) 为输入字母表大小（如文本串和模式串仅包含小写英文字母时，输入字母表大小就是 26）。
说明 2：尽管 Rabin-Karp 算法的匹配时间（最坏情况）为 \(O(nm)\) ，但它的期望匹配时间为 \(O(n)\) 。

Rabin-Karp 算法思想：首先设计一个 hash 函数，在文本串中每次取和模式串相同长度的子串，对其计算 hash，并和模式串的 hash 比较，如果 hash 相等（hash 不相同，两个串肯定不会相同），则用朴素的字符串匹配算法作进一步的验证，验证通过则说明找到匹配。

为了方便讨论，使用记号 \(T[i..j]\) 表示字符串 T 中下标从 \(i\) 到 \(j\) （包含 \(i\) 和 \(j\) ）的子串。

下面是 Rabin-Karp 算法的一个 C 语言实现。

//字符串匹配的Rabin-Karp算法实现
//代码来自http://www-igm.univ-mlv.fr/~lecroq/string/node5.html
//返回第一个匹配位置，找不到就返回NULL
//参数T为文本串指针，len1为其长度；P为模式串指针，len2为其长度。
char *Rabin_Karp(const char *T, const int len1, const char *P, const int len2) {
    assert(T != NULL && P != NULL && (len1 >= len2));
    int d, hashT, hashP, i, j;

    /* Preprocessing */
    for (d = i = 1; i < len2; ++i) { //for循环结束后，d = 2^(len2-1)，当len2很大时，d可能为0。
        d = (d << 1);
    }
    for (hashT = hashP = i = 0; i < len2; ++i) {
        hashT = ((hashT << 1) + T[i]);   //文本串中前len2个字符的hash值。
        hashP = ((hashP << 1) + P[i]);   //模式字符串的hash值。
    }

    /* Searching */
    j = 0;
    while (j <= len1 - len2) {
        if (hashP == hashT) {    //如果hash相同，就用朴素字符串匹配算法进一步验证
            for (i = 0; i < len2; i++) {
                if (P[i] != T[i + j]) {
                    break;
                }
            }
            if (i == len2) {
                return (char *) T + j;
            }
        }
        hashT = ((hashT - T[j] * d) << 1) + T[j + len2]; //计算文本串中下一个len2长的子串的hash
        ++j;
    }
    return NULL;
}

使用的 Rolling hash 函数为：
\[hash(w[0..m-1]) = \left(w[0] 2^{m-1} + w[1] 2^{m-2} + \cdots + w[m-1] 2^{0} \right) \mod INT\_MAX\]
上面 C 语言实现的代码中没有显式地使用模运算，这是因为使用 INT_MAX 作为模运算的模数时，硬件对整数操作的溢出相当于模运算。

记文本串长度为 \(n\) ，模式串长度为 \(m\) 。Rabin-Karp 算法的预处理（计算模式串 hash 值及文本串第一个 hash 值）时间为 \(\Theta(m)\) 。在最坏情况下，它的匹配时间是 \(\Theta((n-m+1)m)\) ，由于 \(m\) 一般会远小于 \(n\) ，又可记为 \(O(nm)\) ，因为 Rabin-Karp 算法和朴素字符串匹配算法一样，对每个有效偏移进行显式验证，最坏情况就是所有尝试都是 hash 的伪命中点，不得不对每次尝试都用朴素字符串匹配算法做进一步验证。

当仔细地选择一个模数时，最坏情况几乎不会出现，其伪命中点会很少，Rabin-Karp 算法的期望运算时间为 \(O(n+m)\) ，由于 \(m \le n\) ，即期望运算时间为 \(O(n)\) 。

假设文本长度为 \(n\) ，模式长度为 \(m\) ，则自动机将会有 \(0,1, \cdots, m\) 这么多种状态（初始状态为 0，接受状态为 \(m\) ）。

假设自动机已经建立好（暂时不考虑如何从模式串构造自动机），用自动机进行字符串匹配的过程如下： 从文本从左到右扫描时，对于每一个字符 \(x\) ，根据自动机当前的状态还有 \(x\) 的值可以找出自动机的下一个状态，这样一直扫描下去，当自动机状态值变为 \(m\) 时我们就可以认为在文本串中找到了一处模式串。

下面是一个简单例子，用于演示利用有限自动机进行字符串匹配的过程。

假设文本和模式只有三种字符 a,b,c（即输入字母表仅这三个字符），已知文本 T 为“abababaca”，模式 P 为“ababaca”，根据模式 P 建立的自动机如图 2 中子图(a)所示（构造过程暂不考虑，后文将介绍）。

在图 2 的子图(c)中，对照自动机状态转换表（即子图(b)），一个个地扫描文本字符，扫描前状态值初始化为 0，这样在 \(i = 9\) 的时候状态值刚好变成 \(7\) （它是一个接受状态），所以找到一个匹配。

我们再看一个例子，还是图 2 中的自动机，但把文本串变为 T=cacbababaca，这时自动机在文件串上的操作过程如图 3 所示。

匹配过程可以用下面伪代码表示：

// 自动机进行文本匹配的过程（伪代码）
Finite_Automation(T, g, m)  //T为文本字符串，g为变迁函数，m为模式P的长度
{
    n=strlen(T);
    q=0;                 //q表示状态，初始为0

    for i=1 to n {
        q=g(q, T[i]);     //通过状态转移表做状态变迁
        if (q==m) {       //判断是否为接受状态
            print "Pattern occurs with shift" i-m  //如果成功匹配，就输出初始位置
        }
    }
}

显然，它的匹配时间为 \(\Theta(n)\) 。

自动机的匹配过程介绍完了，现在问题只剩下如何根据给出的模式 P 计算出对应的自动机，请看下文。

如何模式串 P 计算出对应的自动机呢？在后面的说明中还是使用上面介绍的例子，即建立模式串 P=“ababaca”的有限自动机。

假设文本长度为 \(n\) ，模式长度为 \(m\) ，则自动机将会有 \(0,1, \cdots, m\) 这么多种状态（初始状态为 0，接受状态为 \(m\) ）。首先需要明白一点： 如果当前的状态值为 \(k\) ，其实说明当前文本的后缀与模式串的前缀的最大匹配长度为 \(k\) ，这时读入下一个文本字符时，即使该字符匹配，状态值也最多变成 \(k + 1\) 。

我们以“当前状态为 5，如何计算下一个状态”作为例子，来介绍自动机的构造过程。

假设当前状态为 5，则说明文本当前位置的最后 5 位为“ababa”，它就是模式串的前 5 位。要考虑下一位转移状态，我们需要考虑输入字母表（一般用 \(\sum\) 表示，这里仅为 a,b,c）中的所有字符。
(1) 如果下一位文本字符是“c”，则状态值就可以更新为 6。
(2) 如果下一位文本字符是“a”，这时我们 需要重新找到文本后缀与模式前缀的最大匹配长度。简单的寻找方法可以是令 \(k = 5 + 1=6\) （状态值最大的情况），判断文本后 \(k\) 位与模式串前 \(k\) 位是否相等，不相等的话就令 \(k = k - 1\) 后再继续找。 由于刚才文本后 5 位“ababa”其实就是模式的前 5 位，所以在构建自动机时不需要用到文本串。这个例子中，经过几次迭代测试后 \(k\) 最终为 1（只有“a”1 位相等）。
(3) 如果下一位文本字符是“b”，和上面的分析一样。这个例子中， \(k\) 为 4（模式串前 4 位和文本串后 4 位相等，都为“abab”）。

上面的分析，就是当前状态为 5 时，下一个状态值的转换情况。如图 4 所示。利用这个思路，我们可以计算出所有状态（即状态 \(0,1, \cdots, m\) ，共 \(m+1\) 个状态）的转移情况。

上面介绍了自动机的构造过程，这个过程可用伪代码表示如下：

// 自动机的构造过程：根据给定模式 P[1..m] 来计算转移函数 δ
Compute_Transition_Function(P, Σ)  //模式P，输入的字符表Σ
{
    m=strlen(P);
    for q=0 to m {              //从状态0一直计算到状态m
        for each character x in Σ {  //枚举每个可能的输入字符
            k=min(m, q+1) + 1;  //k的最大可能值
            repeat {
                k=k-1           //重复k=k-1，直到“P的长度为k的前缀”是“P的长度为q的前缀加上字符x”的真后缀
            } until (Pk is the suffix of Pq+x)
            δ(q, x)=k;
        }
    }
    return δ;
}

上面构造自动机算法的运行时间为 \(O(m^3 |\sum|)\) ，这是因为：

for q=0 to m                      // 消耗时间 m
for each character x in Table     // 消耗时间 |Σ|
repeat k=k-1                      // 消耗时间 m
test 'Pk is the suffix of Pq+x'   // 消耗时间 m

上面各个部分都是嵌套的，所以运行时间为 \(O(m^3 |\sum|)\) 。有改进的方法可以在 \(O(m |\sum|)\) 时间内根据模式串构造好自动机，这里不介绍（参考《算法导论(第 3 版)》练习 32.4-8）。

什么是模式串的前缀函数呢？

考察一下朴素算法的匹配过程。在图 5 中，T为文本串，P为模式串，匹配的字符被打上阴影且用垂直直线连接。如果发现不匹配后，则在文本串的 \(s+1\) 处进行下一轮的匹配，当然马上发会现不能匹配。

其实，我们根据前面 5 个匹配字符的已有信息，可以提前知道文本串 \(s+1\) 处的匹配肯定不会成功，这样就根本用不着在 \(s+1\) 处去进行匹配，可以直接在 \(s+2\) 处匹配。更重要的是，我们还可以根据已有信息提前知道在 \(s+2\) 处匹配时，其前 3 个字符能肯定匹配上（因为在模式串的前 5 个字符构造的子串中，其前 3 个字符和后 3 个字符是相同的），也可以跳过，如图 6 所示。

综上所述， 我们利用已有信息可以跳过 5 个（2+3，2 表示下一个可能的有效位移，3表示在这个有效位移上还可以进一步忽略比较的字符数）字符的比较，可以直接拿下一个文本字符和第 4 个（即 3 的后一个）模式字符进行比较。

模式串子串的“真后缀的最长前缀长度”就是当前模式串子串的前缀函数值，所有前缀函数值组成了前缀函数，一般记为 \(\pi\) ，可保存在数组中。上面例子中，模式串子串“ababa”的前缀函数值为 3，记为 \(\pi[5]=3\) （这里，约定首位下标从 1 开始，而不是从 0 开始）。

我们知道，对于朴素算法，文本串和模式串在某位失配后，文本串退回到当次比较最开始时位置的后一位进行下一次尝试，模式串从首位重新开始。

假定已经求得了前缀函数，KMP 算法的匹配过程和朴素算法非常类似，区别在于失配后如何处理：对于 KMP 算法，假设在模式串 P 的 \(k\) 位失配，如果 \(k\) 位是模式串首位，则文本串向右滑动一位后重新开始和模式串比较；如果 \(k\) 位不是模式串首位，则文本串不动（文本串不用回溯！），模式串的下一次尝试的位置是 \(P[1..k-1]\) 对应的前缀值的后面一位。 对于图 6 所示例子，失配位置是模式串第 6 位（不是模式串首位），文本串下一次尝试还是当前失配时的位置，而模式串的下一次尝试是 \(P[1..5]\) 对应的前缀值（即 3）的后一位，即第 4 个模式字符。

经过上面分析知，通过修改朴素算法失配时文本串和模式串的位移情况即可得到 KMP 算法实现。参见下文。

以模式串“ababaca”的子串“ababa”为例。下面演示如何得到 \(\pi[5]=3\) ，进而得到所有的前缀函数值。

计算过程如下： “ababa”的真后缀（有点类似集合里真子集的概念）有Φ、“a”、“ba”、“aba”、“baba”，这些真后缀当中同时又是“ababa”的前缀的有Φ、“a”、“aba”，其中最长的前缀的长度为 3，所以 \(\pi[5]=3\) （约定首位下标从 1 开始，而不是从 0 开始）。

模式串“ababaca”的前缀有：“a”、“ab”、“aba”、“abab”、“ababa”、“ababac”、“ababaca”（不算空前缀），用上面介绍的计算方法可以求得每个前缀的“真后缀的最长前缀长度”分别为：0、0、1、2、3、0、1。所以模式串“ababaca”的完整的前缀函数 \(\pi\) 用表格表示如表 2 所示。

利用节 4.2 介绍的方法，可以求得模式串“ababaca”的 DFA 如图 12 所示。

模式串“ababaca”的 DFA 还可以表示为图 13 所示。

图 13 中，红线表示失配时状态转移方向，每个状态（从状态 1 开始）失配时转换到的状态依次为“0012301”。由前面的计算知，它实际上就是“ababaca”的前缀函数（参见表 2 ）。

下面讨论对 KMP 算法的两个小改进。