KMP算法与一个经典概率问题

    考虑一个事件,它有两种概率均等的结果。比如掷硬币,出现正面和反面的机会是相等的。现在我们希望知道,如果我不断抛掷硬币,需要多长时间才能得到一个特定的序列。

序列一:反面、正面、反面
序列二:反面、正面、正面

    首先,我反复抛掷硬币,直到最近的三次抛掷结果形成序列一,然后我记下这次我抛掷了多少次才得到了我要的序列。重复执行这个过程,我可以算出得到序列一平均需要的抛掷次数。同样地,反复抛掷硬币直到序列二产生,它所需要的次数也有一个平均值。你认为这两个平均值哪一个大哪一个小?换句话说,出现序列一平均所需的抛掷次数少还是出现序列二平均需要的次数少?

    大多数人会认为,两个序列会以同样快的速度出现,因为在所有“正”和“反”的8种三元组合里,“反正反”和“反正正”各占1/8,其概率是均等的。而事实上,我们将会看到掷出序列二所需的次数更少一些。不妨考虑这样一个问题:在由“正”和“反”构成的n位01序列中,有多少个序列以序列一结尾但之前不曾出现过序列一?有多少个序列以序列二结尾但之前不曾出现过序列二?当n比较小时,两者答案是一样的(例如n=3时符合要求的情况都是唯一的),但到后来n越大时,两者的差距越明显:后者的个数总比前者的个数要多一些。不妨看一看n=6的情况。对于序列一,只有以下5个序列是符合要求的:

  • 反反反反正反
  • 反正正反正反
  • 正正正反正反
  • 正反反反正反
  • 正正反反正反

    但对于序列二来说,符合条件的序列就有7个:

  • 反反反反正正
  • 反正反反正正
  • 反反正反正正
  • 正反反反正正
  • 正正反反正正
  • 正正正反正正
  • 正反正反正正

    你可以通过计算机编程枚举,计算一下n为其它值的情况。计算结果和刚才也一样:在n位01序列中,以序列二结尾但之前不含序列二的情况不会少于以序列一结尾但之前不含序列一的情况。这说明,抛掷第n次硬币后恰好出现了序列二,其概率不会小于恰好出现序列一的概率。显然,当n渐渐增大时,这个概率应该呈下降趋势;同时,随着n的增长,两个序列各自出现的概率由相等开始慢慢拉开差距,第n次抛掷产生序列二的概率下降得要缓慢一些,或者说更多的情况集中发生在n更小的时候。因此总的来说,出现序列二所需要的抛掷硬币次数的期望值更小。
    虽然我们通过一系列的观察验证了这个结论,并且我们也相信这个结论是正确的(虽然没有严格的证明),但我们仍然不是很接受这个结论。这种情况是有悖于我们的直觉的,它与我们的生活经验不相符合。此刻,我们迫切需要一个解释,来说明这种出人意料的反常现象产生的原因。

    如果不亲自做几次试验的话,你很难体会到这种微妙的差距。考虑整个游戏的实际过程,“反正正”序列显然会出现得更早一些。假如某一次我们得到了序列“反正”。如果我们需要的是“反正反”序列,那么下一次抛掷结果为反面将结束本轮的抛掷,而下一次是正面则前功尽弃,你必须再次从零开始。如果我们需要的是“反正正”序列,那么下一次抛掷结果为正面将结束本轮的抛掷,而下一次是反面的话我至少不会惨到一切归零,这相当于我已经有了一个反面作为新的开头,只需再来两个正面即可。这样看的话,提前掷出“反正正”的可能性更大一些。
    反复体会上面的想法,了解KMP算法的网友会恍然大悟:这就是KMP算法的基本思路!考虑这样一个问题:我们在当前字串中寻找子串“反正正”第一次出现的位置。假如当前已经能匹配模式串的前两个字“反正”,主串中的下一个字是“正”则匹配成功,主串的下一个字是“反”则将使模式串的当前匹配位置退到第一个字。考虑一个更复杂的例子:我们希望在主串中寻找子串abbaba,现在已经在主串中找到了abbab。如果主串下一个字符是a,则成功匹配;如果主串下一个字符是b,则模式串最多能匹配到的位置退到了第三个字符,我只需要从abb开始继续匹配,而不必一切从头再来。
    我们可以用KMP算法完美地解决上面的问题。首先预处理出一个数组c,c[i,0]表示模式串匹配到了第i个字符,主串下一个字符为0(反)时,模式串的匹配位置将退到哪里;同样地,c[i,1]表示模式串匹配到了第i个字符,主串下一个字符为1(正)时,新的模式串匹配位置在什么地方。设f[i,j]表示第i次抛掷硬币后恰好匹配到模式串第j位有多少种情况,则f[i,j]=Σf(i-1,k) + Σf(i-1,l),其中k满足c[k,0]=j,l满足c[l,1]=j。将f[i,j]除以2的i次方,我们就得到了相应的概率值。或者更直接地,设P[i,j]表示第i次抛掷硬币后,最远能匹配到的模式串位置是第j位的概率,则P[i,j]=Σ( P(i-1,k)/2 ) + Σ( P(i-1,l)/2 )。注意,我们还应该添加一种特殊的概率值P[i,*],它表示在主串第i个字符以前已经成功匹配过的概率,这样的话下表中每一列的和才能为1。

    来看一看程序的输出结果:
Pattern 1: s[]="aba"
主串位置       1        2       3       4       5       6       7       8       9      10
匹配到s[0]  0.5000  0.2500  0.2500  0.2500  0.2188  0.1875  0.1641  0.1445  0.1270  0.1113
匹配到s[1]  0.5000  0.5000  0.3750  0.3125  0.2813  0.2500  0.2188  0.1914  0.1680  0.1475
匹配到s[2]  0.0000  0.2500  0.2500  0.1875  0.1563  0.1406  0.1250  0.1094  0.0957  0.0840
匹配到s[3]  0.0000  0.0000  0.1250  0.1250  0.0938  0.0781  0.0703  0.0625  0.0547  0.0479
已找到匹配  0.0000  0.0000  0.0000  0.1250  0.2500  0.3438  0.4219  0.4922  0.5547  0.6094

Pattern 2: s[]="abb"
主串位置       1        2       3       4       5       6       7       8       9      10
匹配到s[0]  0.5000  0.2500  0.1250  0.0625  0.0313  0.0156  0.0078  0.
0039  0.0020  0.0010
匹配到s[1]  0.5000  0.5000  0.5000  0.4375  0.3750  0.3125  0.2578  0.2109  0.1719  0.1396
匹配到s[2]  0.0000  0.2500  0.2500  0.2500  0.2188  0.1875  0.1563  0.1289  0.1055  0.0859
匹配到s[3]  0.0000  0.0000  0.1250  0.1250  0.1250  0.1094  0.0938  0.0781  0.0645  0.0527
已找到匹配  0.0000  0.0000  0.0000  0.1250  0.2500  0.3750  0.4844  0.5781  0.6563  0.7207

    这下我们可以清楚地看到,序列二提前出现的概率要大得多。注意到,根据我们的概率定义,表格中每一列的数字之和都是1。同时,倒数第二行的数字之和(有无穷多项)也应该为1,因为最后一行的概率就是倒数第二行的概率值累加的结果,而根据最后一行概率的定义,主串无穷长时已找到匹配的概率应该为1。因此,我们也可以把倒数第二行看作是模式串在主串第i个位置首次匹配成功的概率。我们可以根据这一结果近似地计算出抛掷次数的期望值。

Matrix67原创
转贴请注明出处

KMP算法详解

    如果机房马上要关门了,或者你急着要和MM约会,请直接跳到第六个自然段。

    我们这里说的KMP不是拿来放电影的(虽然我很喜欢这个软件),而是一种算法。KMP算法是拿来处理字符串匹配的。换句话说,给你两个字符串,你需要回答,B串是否是A串的子串(A串是否包含B串)。比如,字符串A="I'm matrix67",字符串B="matrix",我们就说B是A的子串。你可以委婉地问你的MM:“假如你要向你喜欢的人表白的话,我的名字是你的告白语中的子串吗?”
    解决这类问题,通常我们的方法是枚举从A串的什么位置起开始与B匹配,然后验证是否匹配。假如A串长度为n,B串长度为m,那么这种方法的复杂度是O (mn)的。虽然很多时候复杂度达不到mn(验证时只看头一两个字母就发现不匹配了),但我们有许多“最坏情况”,比如,A= "aaaaaaaaaaaaaaaaaaaaaaaaaab",B="aaaaaaaab"。我们将介绍的是一种最坏情况下O(n)的算法(这里假设 m<=n),即传说中的KMP算法。
    之所以叫做KMP,是因为这个算法是由Knuth、Morris、Pratt三个提出来的,取了这三个人的名字的头一个字母。这时,或许你突然明白了AVL 树为什么叫AVL,或者Bellman-Ford为什么中间是一杠不是一个点。有时一个东西有七八个人研究过,那怎么命名呢?通常这个东西干脆就不用人名字命名了,免得发生争议,比如“3x+1问题”。扯远了。
    个人认为KMP是最没有必要讲的东西,因为这个东西网上能找到很多资料。但网上的讲法基本上都涉及到“移动(shift)”、“Next函数”等概念,这非常容易产生误解(至少一年半前我看这些资料学习KMP时就没搞清楚)。在这里,我换一种方法来解释KMP算法。

    假如,A="abababaababacb",B="ababacb",我们来看看KMP是怎么工作的。我们用两个指针i和j分别表示,A[i-j+ 1..i]与B[1..j]完全相等。也就是说,i是不断增加的,随着i的增加j相应地变化,且j满足以A[i]结尾的长度为j的字符串正好匹配B串的前 j个字符(j当然越大越好),现在需要检验A[i+1]和B[j+1]的关系。当A[i+1]=B[j+1]时,i和j各加一;什么时候j=m了,我们就说B是A的子串(B串已经整完了),并且可以根据这时的i值算出匹配的位置。当A[i+1]<>B[j+1],KMP的策略是调整j的位置(减小j值)使得A[i-j+1..i]与B[1..j]保持匹配且新的B[j+1]恰好与A[i+1]匹配(从而使得i和j能继续增加)。我们看一看当 i=j=5时的情况。

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B = a b a b a c b
    j = 1 2 3 4 5 6 7

    此时,A[6]<>B[6]。这表明,此时j不能等于5了,我们要把j改成比它小的值j'。j'可能是多少呢?仔细想一下,我们发现,j'必须要使得B[1..j]中的头j'个字母和末j'个字母完全相等(这样j变成了j'后才能继续保持i和j的性质)。这个j'当然要越大越好。在这里,B [1..5]="ababa",头3个字母和末3个字母都是"aba"。而当新的j为3时,A[6]恰好和B[4]相等。于是,i变成了6,而j则变成了 4:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =     a b a b a c b
    j =     1 2 3 4 5 6 7

    从上面的这个例子,我们可以看到,新的j可以取多少与i无关,只与B串有关。我们完全可以预处理出这样一个数组P[j],表示当匹配到B数组的第j个字母而第j+1个字母不能匹配了时,新的j最大是多少。P[j]应该是所有满足B[1..P[j]]=B[j-P[j]+1..j]的最大值。
    再后来,A[7]=B[5],i和j又各增加1。这时,又出现了A[i+1]<>B[j+1]的情况:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =     a b a b a c b
    j =     1 2 3 4 5 6 7

    由于P[5]=3,因此新的j=3:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =         a b a b a c b
    j =         1 2 3 4 5 6 7

    这时,新的j=3仍然不能满足A[i+1]=B[j+1],此时我们再次减小j值,将j再次更新为P[3]:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =             a b a b a c b
    j =             1 2 3 4 5 6 7

    现在,i还是7,j已经变成1了。而此时A[8]居然仍然不等于B[j+1]。这样,j必须减小到P[1],即0:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =               a b a b a c b
    j =             0 1 2 3 4 5 6 7

    终于,A[8]=B[1],i变为8,j为1。事实上,有可能j到了0仍然不能满足A[i+1]=B[j+1](比如A[8]="d"时)。因此,准确的说法是,当j=0了时,我们增加i值但忽略j直到出现A[i]=B[1]为止。
    这个过程的代码很短(真的很短),我们在这里给出:

j:=0;
for i:=1 to n do
begin
   while (j>0) and (B[j+1]<>A[i]) do j:=P[j];
   if B[j+1]=A[i] then j:=j+1;
   if j=m then
   begin
      writeln('Pattern occurs with shift ',i-m);
      j:=P[j];
   end;
end;

    最后的j:=P[j]是为了让程序继续做下去,因为我们有可能找到多处匹配。
    这个程序或许比想像中的要简单,因为对于i值的不断增加,代码用的是for循环
。因此,这个代码可以这样形象地理解:扫描字符串A,并更新可以匹配到B的什么位置。

    现在,我们还遗留了两个重要的问题:一,为什么这个程序是线性的;二,如何快速预处理P数组。
    为什么这个程序是O(n)的?其实,主要的争议在于,while循环使得执行次数出现了不确定因素。我们将用到时间复杂度的摊还分析中的主要策略,简单地说就是通过观察某一个变量或函数值的变化来对零散的、杂乱的、不规则的执行次数进行累计。KMP的时间复杂度分析可谓摊还分析的典型。我们从上述程序的j 值入手。每一次执行while循环都会使j减小(但不能减成负的),而另外的改变j值的地方只有第五行。每次执行了这一行,j都只能加1;因此,整个过程中j最多加了n个1。于是,j最多只有n次减小的机会(j值减小的次数当然不能超过n,因为j永远是非负整数)。这告诉我们,while循环总共最多执行了n次。按照摊还分析的说法,平摊到每次for循环中后,一次for循环的复杂度为O(1)。整个过程显然是O(n)的。这样的分析对于后面P数组预处理的过程同样有效,同样可以得到预处理过程的复杂度为O(m)。
    预处理不需要按照P的定义写成O(m^2)甚至O(m^3)的。我们可以通过P[1],P[2],…,P[j-1]的值来获得P[j]的值。对于刚才的B="ababacb",假如我们已经求出了P[1],P[2],P[3]和P[4],看看我们应该怎么求出P[5]和P[6]。P[4]=2,那么P [5]显然等于P[4]+1,因为由P[4]可以知道,B[1,2]已经和B[3,4]相等了,现在又有B[3]=B[5],所以P[5]可以由P[4] 后面加一个字符得到。P[6]也等于P[5]+1吗?显然不是,因为B[ P[5]+1 ]<>B[6]。那么,我们要考虑“退一步”了。我们考虑P[6]是否有可能由P[5]的情况所包含的子串得到,即是否P[6]=P[ P[5] ]+1。这里想不通的话可以仔细看一下:

        1 2 3 4 5 6 7
    B = a b a b a c b
    P = 0 0 1 2 3 ?

    P[5]=3是因为B[1..3]和B[3..5]都是"aba";而P[3]=1则告诉我们,B[1]、B[3]和B[5]都是"a"。既然P[6]不能由P[5]得到,或许可以由P[3]得到(如果B[2]恰好和B[6]相等的话,P[6]就等于P[3]+1了)。显然,P[6]也不能通过P[3]得到,因为B[2]<>B[6]。事实上,这样一直推到P[1]也不行,最后,我们得到,P[6]=0。
    怎么这个预处理过程跟前面的KMP主程序这么像呢?其实,KMP的预处理本身就是一个B串“自我匹配”的过程。它的代码和上面的代码神似:

P[1]:=0;
j:=0;
for i:=2 to m do
begin
   while (j>0) and (B[j+1]<>B[i]) do j:=P[j];
   if B[j+1]=B[i] then j:=j+1;
   P[i]:=j;
end;

    最后补充一点:由于KMP算法只预处理B串,因此这种算法很适合这样的问题:给定一个B串和一群不同的A串,问B是哪些A串的子串。

    串匹配是一个很有研究价值的问题。事实上,我们还有后缀树,自动机等很多方法,这些算法都巧妙地运用了预处理,从而可以在线性的时间里解决字符串的匹配。我们以后来说。

    昨天发现一个特别晕的事,知道怎么去掉BitComet的广告吗?把界面语言设成英文就行了。
    还有,金山词霸和Dr.eye都可以去自杀了,Babylon素王道。

Matrix67原创
转贴请注明出处