后缀数组

字符串的后缀是指从字符串的某个位置开始到其末尾的字符串子串。我们认为原串和空串也是后缀。比如字符串 "c", "abc"，""，"cabc" 都是字符串 "cabc" 的后缀。

后缀数组是指将某个字符串的所有后缀按照字典序排序后得到的数组。举个例子，字符串 "abcabc" 的后缀数组为 "", "abc", "abcabc", "bc", "bcabc", "c", "cabc"，注意排列的顺序。然而在实际应用中，并不需要保存所有的后缀字符串，只要记录对应的起始位置即可。也就是说对 "abcabc" 来说，对应的后缀数组为 sa = [6, 3, 0, 4, 1, 5, 2]

后缀数组的计算使用 Manber 和 Myers 发明的算法该算法的时间复杂度是 \( O(n log^{2} n) \)，其中 n 是字符串的长度。

算法的基本思想是倍增。首先计算从每个位置开始的长度为 2 的子串顺序，再利用这个结果计算长度为 4 的子串顺序，接下来计算长度为 8 的子串顺序，不断倍增，直到长度大于等于串长 n，这样就得到了后缀数组。

字符串"abcabc" 后缀数组计算

假设已经求得了长度为 k 的子串顺序，要求长度为 2k 的子串顺序。设 rankk[i] 在排序好后的 rankk 数组中的位置代表了 S[i:i+k] 在长度为 k 的所有子串中的排序好后的位置。举个例子对 S = "abcabc"，k = 1。则 rank1 可以为 [97,98,99,97,98,99,-1] ，及空串取 -1, 其它字符取 ascii 码值。则 rankk[1] = 98 在 [-1,97,97,98,98,99,99] 中排在第四位，则代表 "b" 在 ["", "a","a","b","b","c","c"] 中排在第四位。假设我们已经得到了 rankk ，则可以通过对数对 ( rankk[i], rankk[i+k]) 和 数对 (rankk[j], rankk[j+k]) 的比较替代子串 S[i:i+2k] 和子串 S[j:j+2k] 的比较。

from functools import cmp_to_key
import numpy as np

def construct_sa(S):
    n = len(S)

    sa = np.zeros(n+1,dtype=int)
    sa[-1] = n
    rank = np.zeros(n+1,dtype=int)
    rank[-1] = -1
    #初始长度为1，rank 直接取字符串的编码
    for i in range(0,n):
        rank[i] = ord(S[i])
        sa[i] = i

    def cmpBaseRank(n,rank,k,i,j):
        if rank[i] != rank[j]:
            return 1 if rank[i] > rank[j] else -1
        else:
            ri = rank[i+k] if (i+k<=n) else -1
            rj = rank[j+k] if (j+k<=n) else -1
            return 1 if ri > rj else -1

    k = 1
    tmp = np.zeros(n+1,dtype=int)
    while k <= n:
        sa = sorted(sa, key=cmp_to_key(lambda i,j:cmpBaseRank(n, rank, k,i,j)))

        #由长度 k 的 rank 计算长度 2*k 的 rank
        tmp[sa[0]] = 0
        for i in range(1,n+1):
            tmp[sa[i]] = tmp[sa[i-1]] + (1 if (cmpBaseRank(n,rank,k,sa[i-1],sa[i]) == -1) else 0)
        rank = np.tile(tmp,1)
        k = k * 2

    return sa

高度数组

高度数组指的是由后缀数组中相邻两个后缀的最长公共前缀(LCP)的长度组成的数组。记后缀数组为 sa，高度数组为 lcp，则有后缀 S[sa[i]:] 和 S[sa[i+1]:] 的最长公共前缀的长度为 lcp[i]。求高度数组的时间复杂度为 O(n)。下表是字符串 "abcabc" 对应的后缀数组和高度数组。

i sa[i] lcp[i] S[sa[i]:] S[sa[i+1]:]
0 6 0 (空串) abc
1 3 3 abc abcabc
2 0 0 abcabc bc
3 4 2 bc bcabc
4 1 0 bcabc c
5 5 1 c cabc
5 2 - cabc - 
def construct_lcp(S, sa):
    n = len(S)
    rank = np.zeros(n+1, dtype=int)
    for i in range(0,n+1):
        rank[sa[i]] = i

    h = 0
    lcp = np.zeros(n,dtype=int)
    for i in range(0,n):
        j = sa[rank[i]-1]

        if h>0:
            h -= 1
        while j + h < n and i + h < n:
            if S[j+h] != S[i+h]:
                break
            h += 1
        lcp[rank[i] - 1] = h
    return lcp

最长公共子串

将要求最长公共子串的串 S 和 串 T 用不出现在串 S 和 串 T 的字符 '\0' 连接起来得到串 ST，则通过 ST 的高度数组，可以很方便的得到最长公共子串。时间复杂度为 \( O((|T| + |S|)log^{2}(|T| + |S|)) \) ，主要时间花费在求 ST 的后缀数组上。

def getCommanString(S,T):
    ST = S + '\0' + T

    sa = construct_sa(ST)
    lcp = construct_lcp(ST, sa)

    ans = 0
    pos = len(ST)
    n1 = len(S)
    for i in range(0,pos):
        if (sa[i] < n1) != (sa[i+1] < n1):
            if lcp[i] > ans:
                ans = lcp[i]
                pos = sa[i]
    #print(pos,ans)
    return ST[pos:pos+ans]

测试结果

def test():
    S = "abcdAB"
    T = "abcdDD"
    print(getCommanString(S,T), getCommanString(S,T) == "abcd")

    S = "KLabcBDAD"
    T = "abcdBDAD"
    print(getCommanString(S,T), getCommanString(S,T) == "BDAD")

    S = "BDeios"
    T = "abcdBDAD"
    print(getCommanString(S,T), getCommanString(S,T) == "BD")

    S = ""
    T = ""
    print(getCommanString(S,T), getCommanString(S,T) == "")

    S = "ABDD"
    T = ""
    print(getCommanString(S,T), getCommanString(S,T) == "")

    S = ""
    T = "AB"
    print(getCommanString(S,T), getCommanString(S,T) == "")

    S = "AB"
    T = "AB"
    print(getCommanString(S,T), getCommanString(S,T) == "AB")

test()

最长公共子串测试结果

原创文章，转载请注明出处：https://daiyufish.com/article/zuichanggonggongzichuan/