String Matching Michael Tsai 2013/05/28

問題: 字串比對 陣列T[1..n]中有一個長度為n的字串 陣列P[1..m]中有一個長度為m的字串 𝑚≤𝑛 要在T中找P是否出現 P和T的字串從一個字元的集合Σ中拿出 如: Σ={0,1} 或 Σ={a,b,…,z} Pattern P occurs with shift s in text T (Pattern P occurs beginning at position s+1 in text T) if 𝑇 𝑠+𝑗 =𝑃[𝑗], for 1≤𝑗≤𝑚. If P occurs with shift s in T, we call s a valid shift. Otherwise, we call s an invalid shift. 字串比對問題是要在T中間找到所有P出現的位置(valid shift)

一些定義 Σ ∗ : 所有使用Σ中字元組成的有限長度字串(包括長度為0的空字串) |𝑥|: 字串x的長度 xy: 把字串x和y接起來 (concatenation) 𝑤⊏𝑥: 字串w是字串x的prefix (也就是x=wy, y ∈Σ ∗ ) (𝑤⊏𝑥 表示 𝑤 ≤|𝑥|) 𝑤⊐𝑥: 字串w是字串x的suffix (也就是x=yw, y ∈Σ ∗ ) (𝑤⊐𝑥 表示 𝑤 ≤|𝑥|) 例如: ab ⊏abcca, cca ⊐ abcca 空字串𝜖為任何字串的prefix & suffix 對任何字串x, y和字元a, 𝑥⊐𝑦 iff 𝑥𝑎⊐𝑦𝑎 ⊏和⊐為transitive(具遞移律)的operator

方法一: 笨蛋暴力法 Native-String-Matcher(T, P) N=T.length M=P.length for s=0 to n-m if P[1..m]==T[s+1..s+m] print “Pattern occurs with shift” s n-m+1次 𝑂(𝑚) 𝑂 𝑛−𝑚+1 𝑚

為什麼不好? 因為每次for執行比對, 如果錯了, 這回合的資訊完全丟掉. 例: P=aaab 如果我們發現s=0是valid shift (表示T開頭為aaab), 那麼從之前的結果應該可以知道shift 1, 2, 3 都可以直接跳過, 不需要一一比對. T a b … P a b a b a b a b

方法二: The Rabin-Karp Algorithm 假設Σ= 0,1,…,9 那麼每個長度為k的字串可以想成是一個k位數的十進位數 例如字串”31415”可以想成是十進位數31415 把 𝑡 𝑠 設為代表T[s+1..s+m]的十進位數 p設為代表P的十進位數 那麼 𝑡 𝑠 =𝑝 iff T[s+1..s+m]=P[1..m], 也就是s是valid shift 怎麼從P計算p呢? 𝑝=𝑃 𝑚 +10 𝑃 𝑚−1 +10 𝑃 𝑚−2 +…+10 𝑃 2 +10𝑃 1 … 𝑃 𝑚−1 +10 𝑃 𝑚−2 +…+10 𝑃 2 +10𝑃 1 … P[1] P[2] P[3] … P[m-1] P[m]

方法二: The Rabin-Karp Algorithm 怎麼從P計算p呢? 𝑡 0 =𝑇 𝑚 +10 𝑇 𝑚−1 +10 𝑇 𝑚−2 +…+10 𝑇 2 +10𝑇 1 … 𝑇 𝑚−1 +10 𝑇 𝑚−2 +…+10 𝑇 2 +10𝑇 1 … 𝑡 𝑠+1 =10 𝑡 𝑠 − 10 𝑚−1 𝑇 𝑠+1 +𝑇[𝑠+𝑚+1] 𝑠=0 T[1] T[2] T[3] … T[m-1] T[m] T[m+1] T[m+2] 𝑠=1 整個往右移一格 拿掉最左邊那一格 加上最右邊那一格

方法二: The Rabin-Karp Algorithm 那麼用這個方法要多花少時間呢? (簡易分析版) 𝑝=𝑃 𝑚 +10 𝑃 𝑚−1 +10 𝑃 𝑚−2 +…+10 𝑃 2 +10𝑃 1 … 𝑃 𝑚−1 +10 𝑃 𝑚−2 +…+10 𝑃 2 +10𝑃 1 … 𝑡 0 =𝑇 𝑚 +10 𝑇 𝑚−1 +10 𝑇 𝑚−2 +…+10 𝑇 2 +10𝑇 1 … 𝑇 𝑚−1 +10 𝑇 𝑚−2 +…+10 𝑇 2 +10𝑇 1 … 然後用 𝑡 𝑠+1 =10 𝑡 𝑠 − 10 𝑚−1 𝑇 𝑠+1 +𝑇[𝑠+𝑚+1], 算 𝑡 1 , 𝑡 2 ,…, 𝑡 𝑛−𝑚 的值 (每次都是constant time, 共n-m次) 所以總共: 𝑂(𝑚) preprocessing時間, 𝑂(𝑛−𝑚)比對時間 𝑂(𝑚) 𝑂(𝑚)

方法二: The Rabin-Karp Algorithm 之前的兩個問題: Σ如果是general的character set, 怎麼辦? (不再是{0,1,…,9}) 如何解決呢? 假設 Σ =d, Σ={ 𝑐 1 , 𝑐 2 , 𝑐 3 ,…, 𝑐 𝑑 }. 可以把之前的式子改成 𝑝=𝑖𝑑𝑥(𝑃[𝑚])+𝑑 𝑖𝑑𝑥 𝑃 𝑚−1 +𝑑 𝑖𝑑𝑥(𝑃 𝑚−2 )+…+𝑑 𝑖𝑑𝑥 𝑃 2 +𝑑⋅𝑖𝑑𝑥(𝑃 1 ) … 𝑖𝑑𝑥 𝑃 𝑚−1 +𝑑 𝑖𝑑𝑥(𝑃 𝑚−2 )+…+𝑑 𝑖𝑑𝑥 𝑃 2 +𝑑⋅𝑖𝑑𝑥(𝑃 1 ) … (a) 把整個string看成一個d進位的數. (b)字元在Σ 中index當作該字元所代表的的值

方法二: The Rabin-Karp Algorithm 當m比較大的時候, p和 𝑡 𝑠 將很難用電腦直接處理 (用long long也存不下) 加一加乘一乘最後總是會overflow 如何解決? 利用同餘理論. Michael Rabin Richard Karp

同餘理論 (Modular Arithmetic) 假設a, b都為整數 𝑎≡𝑏 𝑚𝑜𝑑 𝑛 : 表示a和b除以n的餘數相等 例如: 38≡14 𝑚𝑜𝑑 12 0≡12 𝑚𝑜𝑑 12 −13≡−3 𝑚𝑜𝑑 5 更棒的性質: 𝑎 1 ≡ 𝑏 1 𝑚𝑜𝑑 𝑛 𝑎 2 ≡ 𝑏 2 𝑚𝑜𝑑 𝑛 則 𝑎 1 + 𝑎 2 ≡ 𝑏 1 + 𝑏 2 𝑚𝑜𝑑 𝑛 𝑎 1 − 𝑎 2 ≡ 𝑏 1 − 𝑏 2 𝑚𝑜𝑑 𝑛 𝑎 1 𝑎 2 ≡ 𝑏 1 𝑏 2 𝑚𝑜𝑑 𝑛

同餘理論 (Modular Arithmetic) 證明: 𝑎 1 ≡ 𝑏 1 𝑚𝑜𝑑 𝑛 表示 𝑎 1 = 𝑐 𝑎 1 𝑛+ 𝑟 1 𝑏 1 = 𝑐 𝑏 1 𝑛+ 𝑟 1 𝑎 2 ≡ 𝑏 2 𝑚𝑜𝑑 𝑛 表示 𝑎 2 = 𝑐 𝑎 2 𝑛+ 𝑟 2 𝑏 2 = 𝑐 𝑏 2 𝑛+ 𝑟 2 所以 𝑎 1 + 𝑎 2 = 𝑐 𝑎 1 + 𝑐 𝑎 2 𝑛+ 𝑟 1 + 𝑟 2 𝑏 1 + 𝑏 2 = 𝑐 𝑏 1 + 𝑐 𝑏 2 𝑛+ 𝑟 1 + 𝑟 2 兩者餘數相同! 𝑎 1 + 𝑎 2 ≡ 𝑏 1 + 𝑏 2 𝑚𝑜𝑑 𝑛 得證. 𝑎 1 ≡ 𝑏 1 𝑚𝑜𝑑 𝑛 𝑎 2 ≡ 𝑏 2 𝑚𝑜𝑑 𝑛 則 𝑎 1 + 𝑎 2 ≡ 𝑏 1 + 𝑏 2 𝑚𝑜𝑑 𝑛

同餘理論 (Modular Arithmetic) 證明: 𝑎 1 ≡ 𝑏 1 𝑚𝑜𝑑 𝑛 表示 𝑎 1 = 𝑐 𝑎 1 𝑛+ 𝑟 1 𝑏 1 = 𝑐 𝑏 1 𝑛+ 𝑟 1 𝑎 2 ≡ 𝑏 2 𝑚𝑜𝑑 𝑛 表示 𝑎 2 = 𝑐 𝑎 2 𝑛+ 𝑟 2 𝑏 2 = 𝑐 𝑏 2 𝑛+ 𝑟 2 所以 𝑎 1 𝑎 2 =( 𝑐 𝑎 1 𝑛+ 𝑟 1 ) 𝑐 𝑎 2 𝑛+ 𝑟 2 = 𝑐 𝑎 1 𝑐 𝑎 2 𝑛+ 𝑐 𝑎 1 𝑟 2 + 𝑐 𝑎 2 𝑟 1 𝑛+ 𝑟 1 𝑟 2 𝑏 1 𝑏 2 = 𝑐 𝑏 1 𝑐 𝑏 2 𝑛+ 𝑐 𝑏 1 𝑟 2 + 𝑐 𝑏 2 𝑟 1 𝑛+ 𝑟 1 𝑟 2 兩者餘數相同! 得證! 𝑎 1 ≡ 𝑏 1 𝑚𝑜𝑑 𝑛 𝑎 2 ≡ 𝑏 2 𝑚𝑜𝑑 𝑛 則 𝑎 1 𝑎 2 ≡ 𝑏 1 𝑏 2 𝑚𝑜𝑑 𝑛

The Rabin-Karp Algorithm 修正版 取q使得dq可以用一個電腦word (32-bit or 64-bit)來表示 既然mod後再加, 減, 乘也會保持原本的關係, 我們可以把這些operation都變成mod版本的 𝑝=𝑃 𝑚 +𝑑 𝑃 𝑚−1 +𝑑 𝑃 𝑚−2 +…+𝑑 𝑃 2 +𝑑𝑃 1 … 𝑃 𝑚−1 +𝑑 𝑃 𝑚−2 +…+𝑑 𝑃 2 +𝑑𝑃 1 … 𝑚𝑜𝑑 𝑞 𝑡 0 =𝑇 𝑚 +𝑑 𝑇 𝑚−1 +𝑑 𝑇 𝑚−2 +…+𝑑 𝑇 2 +𝑑𝑇 1 … 𝑇 𝑚−1 +𝑑 𝑇 𝑚−2 +…+𝑑 𝑇 2 +𝑑𝑇 1 … (𝑚𝑜𝑑 𝑞) 𝑡 𝑠+1 =𝑑 𝑡 𝑠 − 𝑑 𝑚−1 𝑇 𝑠+1 +𝑇 𝑠+𝑚+1 (𝑚𝑜𝑑 𝑞)

The Rabin-Karp Algorithm 修正版 新的mod版algorithm會造成一個問題: 雖然 𝑡 𝑠 =𝑝⇒ 𝑡 𝑠 ≡𝑝 𝑚𝑜𝑑 𝑞 但 𝑡 𝑠 =𝑝⇍ 𝑡 𝑠 ≡𝑝 𝑚𝑜𝑑 𝑞 例如38≡14 𝑚𝑜𝑑 12 , 38≠14 但是如果 𝑡 𝑠 ≢𝑝 𝑚𝑜𝑑 𝑞 ⇒ 𝑡 𝑠 ≠𝑝 所以演算法變成這樣: 如果 𝑡 𝑠 ≢𝑝 𝑚𝑜𝑑 𝑞 , 那麼現在這個s為invalid shift 如果 𝑡 𝑠 ≡𝑝 𝑚𝑜𝑑 𝑞 , 那麼必須額外檢查 (直接比對範圍內的字串花很多時間) 當 𝑡 𝑠 ≡𝑝 𝑚𝑜𝑑 𝑞 , 但是 𝑡 𝑠 ≠𝑝時, 稱為spurious hit 當q夠大的時候, 希望spurious hit會相當少

例子: Rabin-Karp Algorithm

Pseudo Code: Rabin-Karp Rabin-Karp-Matcher(T,P,d,q) n=T.length m=P.length h= 𝑑 𝑚−1 mod q p=0 t=0 for i=1 to m p=(dp+P[i]) mod q t=(dt+T[i]) mod q for s=0 to n-m if p==t if P[1..m]==T[s+1..s+m] print “Pattern occurs with shift” s if s<n-m t=(d(t-T[s+1]h)+T[s+m+1]) mod q T: string to be searched P: pattern to be matched d: size of the character set q: max number Pre-processing: 𝑂 𝑚 迴圈跑n-m+1次 Hit的時候比對: 𝑂 𝑚

Worst-case Running Time Rabin-Karp-Matcher(T,P,d,q) n=T.length m=P.length h= 𝑑 𝑚−1 mod q p=0 t=0 for i=1 to m p=(dp+P[i]) mod q t=(dt+T[i]) mod q for s=0 to n-m if p==t if P[1..m]==T[s+1..s+m] print “Pattern occurs with shift” s if s<n-m t=(d(t-T[s+1]h)+T[s+m+1]) mod q Worst case的時候: T= 𝑎 𝑛 (n個a) P= 𝑎 𝑚 (m個a) 比對的時間為 O(m(n-m+1)) Pre-processing: 𝑂 𝑚 迴圈跑n-m+1次 Hit的時候比對: 𝑂 𝑚

Average Running Time Rabin-Karp-Matcher(T,P,d,q) n=T.length m=P.length h= 𝑑 𝑚−1 mod q p=0 t=0 for i=1 to m p=(dp+P[i]) mod q t=(dt+T[i]) mod q for s=0 to n-m if p==t if P[1..m]==T[s+1..s+m] print “Pattern occurs with shift” s if s<n-m t=(d(t-T[s+1]h)+T[s+m+1]) mod q 平常的時候, valid shift很少 (假設有c個) 不會每次都有modulo的hit. 假設字串各種排列組合出現的機率相等，則spurious hit的機率可當成1/𝑞. Pre-processing: 𝑂 𝑚 則比對花的時間: Spurious hit共花O( (n-m+1)/q)=O(n/q)次 總共比對花的時間為 O((n-m+1)+(m(c+n/q))) 迴圈跑n-m+1次 If c=O(1) and q≥m,  O(n+m)=O(n) Hit的時候比對: 𝑂 𝑚

方法三: The Knuth-Morris-Pratt Algo.

Knuth-Morris-Pratt Don Knuth James Morris Vaughan Pratt

正式一點的說法 If 𝑃 𝑞 是 𝑇 𝑠+𝑞 的suffix 找 𝑃 𝑞 的最長prefix使得它是 𝑇 𝑠+𝑞 的suffix 假設P[1..q]和T[s+1..s+q]已經match了 要找出最小的shift s’使得某個k<q可以滿足 𝑃 1..𝑞 =𝑇 𝑠 ′ +1.. 𝑠 ′ +𝑘 , 𝑠 ′ +𝑘=𝑠+𝑞 If 𝑃 𝑞 是 𝑇 𝑠+𝑞 的suffix s’=s+(q-k) 找 𝑃 𝑞 的最長prefix使得它是 𝑇 𝑠+𝑞 的suffix

最好的狀況下: 沒有重複的pattern … a b c f a b c d s q … a b c f a b c d s’=s+q k=0

先處理P來取得”重複pattern”的資訊 找 𝑃 𝑞 的最長prefix使得它是 𝑇 𝑠+𝑞 的suffix 找 𝑃 𝑞 的最長prefix 𝑃 𝑘 使得它是 𝑃 𝑞 的suffix 定義Prefix function (failure function) 𝜋: Input: {1,2,…,m} Output: {0,1,…,m-1} 𝜋 𝑞 = max 𝑘:𝑘<𝑞 and 𝑃 𝑘 is a suffix of 𝑃 𝑞 (也就是前面例子中的k值, 可以想成最長的重複pattern的長度)

Prefix function example 𝜋 𝑞 = max 𝑘:𝑘<𝑞 and 𝑃 𝑘 is a suffix of 𝑃 𝑞 i 1 2 3 4 5 6 7 P[i] A B C 𝜋[𝑖] 𝜋[𝑖] 1 2 3 i 1 2 3 4 5 6 7 8 9 10 11 12 P[i] A B C 𝜋[𝑖] 𝜋[𝑖] 1 2 3 4 5

Pseudo-code: Prefix function Compute-Prefix-Function(P) m=P.length let 𝜋[1..𝑚] be a new array 𝜋 1 =0 k=0 for q=2 to m while k>0 and P[k+1]!=P[q] k=𝜋 𝑘 if P[k+1]==P[q] k=k+1 𝜋 𝑞 =𝑘 return 𝜋

例子: Matching T=BACBABABAABCBAB 實際的Matching Pseudo Code和計算prefix function非常像 請見Cormen p. 1005 KMP-Matcher i 1 2 3 4 5 6 7 P[i] A B C 𝜋[𝑖]

算Prefix function花多少時間? Compute-Prefix-Function(P) m=P.length let 𝜋[1..𝑚] be a new array 𝜋 1 =0 k=0 for q=2 to m while k>0 and P[k+1]!=P[q] k=𝜋 𝑘 if P[k+1]==P[q] k=k+1 𝜋 𝑞 =𝑘 return 𝜋 共𝑂 𝑚 次 麻煩的是這邊: 總共會跳多少次呢?

算Prefix function花多少時間? Compute-Prefix-Function(P) m=P.length let 𝜋[1..𝑚] be a new array 𝜋 1 =0 k=0 for q=2 to m while k>0 and P[k+1]!=P[q] k=𝜋 𝑘 if P[k+1]==P[q] k=k+1 𝜋 𝑞 =𝑘 return 𝜋 進入迴圈的時候k<q, 且q每次增加, k有時候不增加 所以k<q永遠成立 所以𝜋 𝑞 =𝑘<q. 所以每執行一次迴圈就減少k一次 且k永遠不是負的 Total:O(m) k只會在這邊增加, 因此最多總共增加m-1次(迴圈執行次數) 最後: 既然有增加才有得減少, while loop總共執行的次數不會超過O(m)

KMP執行時間 類似的方法可以證明比對的部分執行時間為O(n) 所以總和來看: Preprocessing時間 O(m) 比對時間 O(n)

Today’s Reading Assignment Cormen ch. 32, 32.1, 32.2, 32.4 (正確性的證明略為複雜)