后缀数组【原理+python代码】

 

后缀数组

参考：https://blog.csdn.net/a1035719430/article/details/80217267
https://blog.csdn.net/YxuanwKeith/article/details/50636898

主要由三个数组组成，一部分是sa和rank，另一部分是height

sa与rank的产生：倍增法

还有其它方法，但是倍增法虽然时间复杂度低，但是已经算好理解的了……

当i在[0,n)之间变化时，s[i:]就代表着数组的n个后缀。
现在我们要对这n个后缀按字典序进行排序，比如对于"banana$"，排好序之后长这样：

现在我们用两个数组来记录这种排序，那就是sa和rank

• sa[i]表示排名为i的后缀的起始位置，它的下标是排名，值是位置，意思就是【s[sa[i]:]排在第i个】。
• rank[i]表示以i起始的后缀的排名，它的下标是位置，值是排名，意思就是【s[i:]的排名为rank[i]】。

那么显然会有这样的等量关系：

• rank[sa[i]] = i，这里的i表示的是排名
• sa[rank[i]] = i，这里的i表示的是位置

接下来看看这两数组怎么求。
想象一下这么一个过程：

1. 首先，字符串里的每个字符都代表着它的对应后缀的第一个字母，假如我们只按这个首字母进行排序（先不管如何排序），我们能得到一种排序结果。我们把用来排序的长度叫“有效长度”，现在它为1。

   我们如果能不断扩展这个有效长度，直到它超过n，或者在扩展的过程中发现rank数组已经是0~n-1的一种排列（每个排名都不相同了），这时候我们就排好了

2. 接下来，我们要把有效长度扩展到2，这时候我们要利用后缀的一个很重要的性质：

   后缀s[i:]的前l个字符，恰好是后缀s[i-l:]的[l:2l]部分的字符

   假如从s[l:]到s[n-1:]都已经按有效长度l排好序了，那么后缀s[0:]到s[n-1-l:]的[l:2l]部分的字符的排序就已知了

所以倍增法就是这么一个过程：
当所有后缀已经按有效长度l排好序之后，把所有后缀的[0:l]部分字符当作第一关键字，把[l:2l]部分的字符当作第二关键字，再度进行排序，就能得到有效长度为2l的排序。

而基数排序就特别适合这种模式，按第一关键字的rank值相当于十位，按第二关键字的rank值相当于个位，然后排序。

具体看看如何实现：

初始化

# sa[i]:排名为i的后缀的起始位置
# rk[i]:起始位置为i的后缀的排名
n = len(s)
sa = []
rk = []
for i in xrange(n):
    rk.append(ord(s[i])-ord('a')) # 刚开始时，每个后缀的排名按照它们首字母的排序
    sa.append(i) # 而排名第i的后缀就是从i开始的后缀

倍增循环

l = 0 # l是已经排好序的长度，现在要按2l长度排序
sig = 26 # sig是unique的排名的个数，初始是字符集的大小
while True:
    p = []
    # 对于长度小于l的后缀来说，它们的第二关键字排名肯定是最小的，因为都是空的
    for i in xrange(n-l,n):
        p.append(i)
    # 对于其它长度的后缀来说，起始位置在`sa[i]`的后缀排名第i，而它的前l个字符恰好是起始位置为`sa[i]-l`的后缀的第二关键字
    for i in xrange(n):
        if sa[i]>=l:
            p.append(sa[i]-l)
    # 然后开始基数排序，先对第一关键字进行统计
    # 先统计每个值都有多少
    cnt = [0]*sig
    for i in xrange(n):
        cnt[rk[i]] += 1
    # 做个前缀和，方便基数排序
    for i in xrange(1,sig):
        cnt[i] += cnt[i-1]
    # 然后利用基数排序计算新sa
    for i in xrange(n-1,-1,-1):
        cnt[rk[p[i]]] -= 1
        sa[cnt[rk[p[i]]]] = p[i]
    # 然后利用新sa计算新rk
    def equal(i,j,l):
        if rk[i]!=rk[j]:return False
        if i+l>=n and j+l>=n:
            return True
        if i+l<n and j+l<n:
            return rk[i+l]==rk[j+l]
        return False
    sig = -1
    tmp = [None]*n
    for i in xrange(n):
        # 直接通过判断第一关键字的排名和第二关键字的排名来确定它们的前2l个字符是否相同
        if i==0 or not equal(sa[i],sa[i-1],l):
            sig += 1
        tmp[sa[i]] = sig
    rk = tmp
    sig += 1
    if sig==n:
        break
    # 更新有效长度
    l = l << 1 if l > 0 else 1

然后开始倍增法的循环，直到排好序为止，也就是rank数组里的unique排名数达到n为止
每个循环大概就是这样的思路：

1. 基于老sa数组计算出按第二关键字排序的新sa数组（把它叫做p）
2. 基于p数组和老rank数组，计算出新sa数组
3. 基于新sa数组得到新rank数组，并统计出unique排名数量

p数组的产生

p = []
# 对于长度小于l的后缀来说，它们的第二关键字排名肯定是最小的，因为都是空的
for i in xrange(n-l,n):
    p.append(i)
# 对于其它长度的后缀来说，起始位置在`sa[i]`的后缀排名第i，而它的前l个字符恰好是起始位置为`sa[i]-l`的后缀的第二关键字
for i in xrange(n):
    if sa[i]>=l:
        p.append(sa[i]-l)

基数排序的cnt数组

# 先统计每个值都有多少
cnt = [0]*sig
for i in xrange(n):
    cnt[rk[i]] += 1
# 做个前缀和，方便基数排序
for i in xrange(1,sig):
    cnt[i] += cnt[i-1]

利用基数排序更新排名

# 然后利用基数排序计算新sa
for i in xrange(n-1,-1,-1):
    cnt[rk[p[i]]] -= 1
    sa[cnt[rk[p[i]]]] = p[i]

从后往前更新，按第二关键字排序排在最后的后缀，肯定在它的第一关键字的组里排最后，也就是前缀和cnt[rk[p[i]]]的值
接着每排一个，我们都将对应的cnt数组-1，相当于指针指向了前一位。
现在我们就把新sa数组计算出来了

根据新sa计算新rank

def equal(i,j,l):
    if rk[i]!=rk[j]:return False
    if i+l>=n and j+l>=n:
        return True
    if i+l<n and j+l<n:
        return rk[i+l]==rk[j+l]
    return False
sig = -1
tmp = [None]*n
for i in xrange(n):
    # 直接通过判断第一关键字的排名和第二关键字的排名来确定它们的前2l个字符是否相同
    if i==0 or not equal(sa[i],sa[i-1],l):
        sig += 1
    tmp[sa[i]] = sig
rk = tmp
sig += 1
if sig==n:
    break

这一步需要判断一下排名相邻的两个数组在有效长度范围内是不是相等的，即它们的[0:2l]部分是不是相等的，如果是相等的那么排名不变，这有助于统计unique排名数量sig，目的是为了跳出循环
为了保证判断相等的时间复杂度为O(1)，需要建立一个新的rk的临时数组，同时利用老rk的值来判断是否相等

综上所述，倍增法的所有代码如下：

def doubling(s):
    # sa[i]:排名为i的后缀的起始位置
    # rk[i]:起始位置为i的后缀的排名
    n = len(s)
    sa = []
    rk = []
    for i in xrange(n):
        rk.append(ord(s[i])-ord('a')) # 刚开始时，每个后缀的排名按照它们首字母的排序
        sa.append(i) # 而排名第i的后缀就是从i开始的后缀

    l = 0 # l是已经排好序的长度，现在要按2l长度排序
    sig = 26 # sig是unique的排名的个数，初始是字符集的大小
    while True:
        p = []
        # 对于长度小于l的后缀来说，它们的第二关键字排名肯定是最小的，因为都是空的
        for i in xrange(n-l,n):
            p.append(i)
        # 对于其它长度的后缀来说，起始位置在`sa[i]`的后缀排名第i，而它的前l个字符恰好是起始位置为`sa[i]-l`的后缀的第二关键字
        for i in xrange(n):
            if sa[i]>=l:
                p.append(sa[i]-l)
        # 然后开始基数排序，先对第一关键字进行统计
        # 先统计每个值都有多少
        cnt = [0]*sig
        for i in xrange(n):
            cnt[rk[i]] += 1
        # 做个前缀和，方便基数排序
        for i in xrange(1,sig):
            cnt[i] += cnt[i-1]
        # 然后利用基数排序计算新sa
        for i in xrange(n-1,-1,-1):
            cnt[rk[p[i]]] -= 1
            sa[cnt[rk[p[i]]]] = p[i]
        # 然后利用新sa计算新rk
        def equal(i,j,l):
            if rk[i]!=rk[j]:return False
            if i+l>=n and j+l>=n:
                return True
            if i+l<n and j+l<n:
                return rk[i+l]==rk[j+l]
            return False
        sig = -1
        tmp = [None]*n
        for i in xrange(n):
            # 直接通过判断第一关键字的排名和第二关键字的排名来确定它们的前2l个字符是否相同
            if i==0 or not equal(sa[i],sa[i-1],l):
                sig += 1
            tmp[sa[i]] = sig
        rk = tmp
        sig += 1
        if sig==n:
            break
        # 更新有效长度
        l = l << 1 if l > 0 else 1
    return sa,rk

height数组

我们现在令height[i]为s[sa[i-1]:]和s[sa[i]:]的最长公共前缀长度，即排名相邻的两个后缀的最长公共前缀长度。
比如：

那么如果rank[j]<rank[k]，则后缀S[j:]和S[k:]的最长公共前缀为min(height[rank[j]+1],height[rank[j]+2]...height[rank[k]])。

证明转载自：https://blog.csdn.net/a1035719430/article/details/80217267，补充了加粗部分的话。

同时，我们还有一个结论：height[rank[i]]≥height[rank[i−1]]−1

证明：
设suffix(k)是排在suffix(i−1)前一名的后缀，则它们的最长公共前缀是height[rank[i−1]]
那么suffix(k+1)将排在suffix(i)的前面
并且suffix(k+1)和suffix(i)的最长公共前缀至少是height[rank[i−1]]−1
那么由于suffix(k+1)和suffix(i)的前height[rank[i−1]]−1位都一样，那么排名在它们中间的后缀，这前height[rank[i−1]]−1位也都得一样，不然它们肯定不会排在中间。
所以suffix(i)和在它前一名的后缀的最长公共前缀至少是height[rank[i−1]]−1
证毕。

这样我们按照height[rank[1]],height[rank[2]]...height[rank[n]]的顺序计算，利用height数组的性质，就可以将时间复杂度可以降为O(n)。这是因为height数组的值最多不超过n，每次计算结束我们只会减1，所以总的运算不会超过2n次。

# 计算height数组
k = 0
height = [0]*n
for i in xrange(n):
    if rk[i]>0:
        j = sa[rk[i]-1]
        while i+k<n and j+k<n and s[i+k]==s[j+k]:
            k += 1
        height[rk[i]] = k
        k = max(0,k-1) # 下一个height的值至少从max(0,k-1)开始

总的代码

def doubling(s):
    # sa[i]:排名为i的后缀的起始位置
    # rk[i]:起始位置为i的后缀的排名
    n = len(s)
    sa = []
    rk = []
    for i in xrange(n):
        rk.append(ord(s[i])-ord('a')) # 刚开始时，每个后缀的排名按照它们首字母的排序
        sa.append(i) # 而排名第i的后缀就是从i开始的后缀

    l = 0 # l是已经排好序的长度，现在要按2l长度排序
    sig = 26 # sig是unique的排名的个数，初始是字符集的大小
    while True:
        p = []
        # 对于长度小于l的后缀来说，它们的第二关键字排名肯定是最小的，因为都是空的
        for i in xrange(n-l,n):
            p.append(i)
        # 对于其它长度的后缀来说，起始位置在`sa[i]`的后缀排名第i，而它的前l个字符恰好是起始位置为`sa[i]-l`的后缀的第二关键字
        for i in xrange(n):
            if sa[i]>=l:
                p.append(sa[i]-l)
        # 然后开始基数排序，先对第一关键字进行统计
        # 先统计每个值都有多少
        cnt = [0]*sig
        for i in xrange(n):
            cnt[rk[i]] += 1
        # 做个前缀和，方便基数排序
        for i in xrange(1,sig):
            cnt[i] += cnt[i-1]
        # 然后利用基数排序计算新sa
        for i in xrange(n-1,-1,-1):
            cnt[rk[p[i]]] -= 1
            sa[cnt[rk[p[i]]]] = p[i]
        # 然后利用新sa计算新rk
        def equal(i,j,l):
            if rk[i]!=rk[j]:return False
            if i+l>=n and j+l>=n:
                return True
            if i+l<n and j+l<n:
                return rk[i+l]==rk[j+l]
            return False
        sig = -1
        tmp = [None]*n
        for i in xrange(n):
            # 直接通过判断第一关键字的排名和第二关键字的排名来确定它们的前2l个字符是否相同
            if i==0 or not equal(sa[i],sa[i-1],l):
                sig += 1
            tmp[sa[i]] = sig
        rk = tmp
        sig += 1
        if sig==n:
            break
        # 更新有效长度
        l = l << 1 if l > 0 else 1
    # 计算height数组
    k = 0
    height = [0]*n
    for i in xrange(n):
        if rk[i]>0:
            j = sa[rk[i]-1]
            while i+k<n and j+k<n and s[i+k]==s[j+k]:
                k += 1
            height[rk[i]] = k
            k = max(0,k-1) # 下一个height的值至少从max(0,k-1)开始
    return sa,rk,height