Rabin-Karp算法 (拉宾-卡普)

Rabin-karp算法是朴素字符串匹配算法的一个特例。当字母表∑为d进制数时，即∑={0,1,2,…d-1}。如当d=10时字母表中的每个字符都是一个十进制数。我们在比较两个长度为m的子串时，可以把这两个子串当作整数进行比较，而不用逐个字符比较，从而在某种程度上减少算法时间。

把一个由d进制数字组成的字符串转换成相应的十进制整数，这是大家曾经都写过的东西，一个可能的简单实现如下：

int ConvertToInt(const char* str, int d){int ans =0;int m = strlen(str);for(int i=0; i<m; i++){ans = ans*d + str[i]-'0';}return ans;}

这就是所谓的霍纳法则(Horner’s rule)

对一个长度为m的子串，求其对应的整数值，其复杂度为O(m)，如果不做一些特殊的处理，每次把子串和模式P比较时，先求子串对应的整数t，再与模式P对应的整数值p比较，算法的复杂度并没有得到改进。如果用t(s)表示当前位移下子串T[s+1..s+m]的值，我们注意到把模式向右滑动一个窗口之后，下一个m位的子串T[s+2..s+m+1]对应的整数值t(s+1)和t(s)相比，只是去掉了最高位数字T[s+1]，增加了一个最低位的数字T[s+m+1]。因而t(s+1)和t(s)之间有如下关系式：

t(s+1)=d*(t(s)-T[s+1]*d^(m-1))+T[s+m+1]

例如对于文本十进制数字组成的文本”2359023141526739921”，d=10，m=5，当s=6时，t(s)=”31415”，T[s+1]=T[7]=3，T[s+m+1]=T[12]=2，所以

t(s+1) =10*(31415-3*10^4)+2=14152

因而我们初始时，只用求p=P[1..m],t0=T[1..m],t(s+1)的求值可以在每次循环迭代中去求，避免了多次函数调用的开销。因而一个可能的实现如下：

void Bad_Rabin_Karp_Matcher(const char* T, const char* P, int d){int n = strlen(T);int m = strlen(P);int h = pow(static_cast<double>(d),m-1);//h=d^(m-1)int p=0;int t=0;for(int i=0; i<m; i++){p = p*d + P[i]-'0';//p=P[0,m-1]t = t*d + T[i]-'0';//t=T[0,m-1]}for(int s=0; s<=n-m; s++){if(p==t){cout<<"Pattern occurs with shift"<<s<<endl;}if(s<n-m){t = d*(t-h*(T[s]-'0')) + T[s+m]-'0';}}}

注：

(1) 在前面讲解过程中数组下标从1开始，而在实际编程中数组下标是从0开始的，所以代码和描述的有点差异，但是原理是一样的。

(2)注意到我们求p,t,h的过程，不管用的是int还是long long型，当m增大时，肯定会产生溢出，比如说一个32位的int型表示的最大整数也就是20多亿。

(3)一般做过一些ACM题目的同学，知道在这种情况我们可以通过对一个数求模来防止溢出，记其为q。

(4)这里又出现了一个问题，当p≠t (mod q)时，的确可以推出p≠q；但是当p＝t (mod q)时，推不出p=t。即可能出现所谓的伪命中点，即p＝t (mod q)但是p≠q。所以当p＝t (mod q)时，我们要进行额外检查，依次比较这m个字符，看其是否真的命中。

(5)通过把q设置为一个较大的素数，可以有效的减少伪命中点数，因而可以减小额外检查的开销。

通过以上分析，改进后一个可能的算法实现如下：

void Rabin_Karp_Matcher(const char* T, const char* P, int d, int q){int n = strlen(T);int m = strlen(P);int p = 0;int t = 0;int h = d;//当m=1时有问题，此处应该为h=1,下面循环中k初始为1for(int k=2; k<m; k++){h = h*d % q;}for(int i=0; i<m; i++){p = (p*d + P[i]-'0') % q;t = (t*d + T[i]-'0') % q;}for(int s=0; s<=n-m; s++){//cout<<"p="<<p<<" t="<<t<<endl;if(p==t){if(strncmp(T+s,P,m)==0){cout<<"Pattern occurs with shift"<<s<<endl;}}if(s<n-m){t = (d*(t-h*(T[s]-'0'))+T[s+m]-'0') % q;//此处也有问题，t可能小于0}}}

但是运行结果和和预期不符，比如说你输入P=”111111”，T=”1”，结果只输出了一个位移值0，输入T=”1234”，P=”34”居然没有有效的位移输出。

分析代码，唯一可能的原因就是求余过程出现到了问题，因此我们在循环中比较p和t的值之前先输出它们的值，便于分析。果不其然，拿P=”111111”，T=”1”的输出作为示例，除了第一个位移p和t相等，后面的位移t居然是负值。也就是说对q求余的结果可能出现[-q+1,-1]范围的负值，为了计算的正确性，我们要保证迭代的过程中t都为正值。当t<0时，只需加上q即可。当然初始计算p和t的值时，p和t是不可能为负值的。修改完代码后，上面举的第二个输入用例运行正确，但是第一个仍然有问题，分析发现，原来当模式P只含一个字符时m=1，h=d^(m-1)=d^0=1，而前面我在求h时，只想到h是m-1个的相乘，把h初始化为了d，疏忽了m可能取1，所以出现了此种情况。到此分析结束，完整正确的函数及测试代码如下：

#include <cstdlib>#include <iostream>#include <cstring>#include <cmath>using namespace std;void Rabin_Karp_Matcher(const char* T, const char* P, int d, int q){int n = strlen(T);int m = strlen(P);int p = 0;int t = 0;int h = 1;for(int k=1; k<m; k++) //预处理{h = h*d % q;}for(int i=0; i<m; i++){p = (p*d + P[i]-'0') % q;t = (t*d + T[i]-'0') % q;}for(int s=0; s<=n-m; s++)//匹配{cout<<"p="<<p<<" t="<<t<<endl;if(p==t){if(strncmp(T+s,P,m)==0){cout<<"Pattern occurs with shift"<<s<<endl;}}if(s<n-m){t = (d*(t-h*(T[s]-'0'))+T[s+m]-'0') % q;if(t<0){t = t+q;}}}}int main(int argc, char *argv[]){const int Max_Length = 1000;const int d = 10;const int q = 13;char T[Max_Length];char P[Max_Length];while(gets(T)){gets(P);Rabin_Karp_Matcher(T,P,d,q);cout<<"next case:"<<endl;}    system("PAUSE");    return EXIT_SUCCESS;}

算法分析：预处理时间Θ(m)，即求h,p,t的时间为，匹配时间在最坏情况下为Θ((n-m-1)m)，因为可能出现每次都是可能命中点的情况。如T=a^n,P=a^m，此种情况下验证时间为Θ((n-m-1)m)。当然实际中，可能的命中点一般很少。假设有c个，则算法的期望匹配时间为O(n-m+1 +cm)=O(m+n)，当m<<n时，期望匹配时间为O(n).