后缀数组 sa 数组计算,倍增算法 - (sunznx) 振翅飞翔
21 September 2019

几年前写的一篇 后缀树 的文章,本文算是填坑之一。

倍增算法计算后缀数组的感受:

脑子:我懂了
手:不,你没懂

接下来,我们先看看如何让脑子懂。以 "aabaaaab" 这个字符串为例。首先,先对,字符串的每个字符做一次 “基数排序/桶排序”,结果如下

然后开始倍增:倍增步数为 1
在这里,我们要计算 "aa ab ba aa aa aa ab b" 的排序关系,因为我们知道 "a a b a a a a b"(设为 x,即第一关键字) 和 "a b a a a a b null" (设为 y,即第二关键字,这里 null 表示空字符串)的排序关系了。我们只要把这两个排序关系 “合并” 起来 然后再执行一次基数排序就可以得出结果

继续倍增:倍增步数为 2
在这里,我们要计算的字符串变成了 "aaba abaa baaa aaaa aaab aab ab b"。同样的,x 为 "aa ab ba aa aa aa ab b",y 为 "ba aa aa aa ab b null null"

继续倍增:倍增步数为 4
在这里,我们要计算的字符串变成了 "aabaaaab abaaaab baaaab aaaab aaab aab ab b"

在开始吹水代码怎么实现之前,有两个知识点是需要了解的

  1. 基数排序
  2. 倍增算法

代码实现:(很难看懂,只能贴代码出来了)

public function __construct($str)
{
    $this->str = $str;
    $this->maxTax = 27;
    $this->prework();
}

private function prework()
{
    $n = \strlen($this->str);
    $this->x = array_fill(0, 2 * $n, 0);
    $this->y = array_fill(0, 2 * $n, 0);
    $this->sa = array_fill(0, $n, 0);

    for ($i = 0; $i < $n; $i++) {
        $this->x[$i] = static::getDigitByRadix($this->str, $i);
    }
    $this->computeSA();
}

private function computeSA()
{
    $n = \strlen($this->str);
    $p = 1;

    while (!$this->isComputeSAok()) {
        for ($i = 0; $i < $n; $i++) {
            $this->y[$i] = $this->x[$i + $p];
        }
        $this->sortByXY();
        $p = $p * 2;
    }

    for ($i = 0; $i < $n; $i++) {
        $this->sa[$i] = $this->x[$i] - 1;
    }
}

// radix 从 0 开始
private static function getDigitByRadix($str, $radix)
{
    return isset($str[$radix]) ? (ord($str[$radix]) - ord('a') + 1) : 0;
}

private function isComputeSAok()
{
    $n = \strlen($this->str);
    $tax = array_fill(0, $this->maxTax + 1, 0);
    for ($i = 0; $i < $n; $i++) {
        $v = $this->x[$i];
        $tax[$v] += 1;
        if ($tax[$v] >= 2) {
            return false;
        }
    }
    return true;
}

prework 函数首先初始化了 x 数组,即上面所说的第一关键字,这部分应该很简单了,接下来就开始计算 sa 数组了。
computeSA 会调用 isComputeSAok 函数来判断 SA 是不是计算完成了,计算完成的依据是 x 数组的每个元素都是不相同的,这意味着已经排好序了。例如字符串 "abcde", x=[1, 2, 3, 4, 5],那么这时候便不需要再计算了。computeSA 会根据 x(第一关键字)p(倍增步长) 来计算出 y(第二关键字) ,代码的关键部分是 sortByXY

 1: private function sortByXY()
 2: {
 3:     $n = \strlen($this->str);
 4:     $taxX = array_fill(0, $this->maxTax + 1, 0);
 5:     $taxY = array_fill(0, $this->maxTax + 1, 0);
 6:     $saY = array_fill(0, $n, 0);
 7:     $saX = array_fill(0, $n, 0);
 8: 
 9:     for ($i = 0; $i < $n; $i++) {
10:         $taxY[$this->y[$i] + 1] += 1;
11:     }
12: 
13:     for ($i = 1; $i < $this->maxTax; $i++) {
14:         $taxY[$i] += $taxY[$i - 1];
15:     }
16: 
17:     for ($i = 0; $i < $n; $i++) {
18:         $digit = $this->y[$i];
19:         $idx = $taxY[$digit]++;
20:         $saY[$idx] = $i;
21:     }
22: 
23:     for ($i = 0; $i < $n; $i++) {
24:         $taxX[$this->x[$i] + 1] += 1;
25:     }
26: 
27:     for ($i = 1; $i < $this->maxTax; $i++) {
28:         $taxX[$i] += $taxX[$i - 1];
29:     }
30: 
31:     for ($i = 0; $i < $n; $i++) {
32:         $r = $saY[$i];
33:         $digit = $this->x[$r];
34:         $idx = $taxX[$digit]++;
35:         $saX[$idx] = $r;
36:     }
37: 
38:     $newX = array_fill(0, 2 * $n, 0);
39: 
40:     $r = 1;
41:     $newX[$saX[0]] = $r++;
42: 
43:     for ($i = 1; $i < $n; $i++) {
44:         $curArrIdx = $saX[$i];
45:         $preArrIdx = $saX[$i - 1];
46:         if (($this->x[$curArrIdx] == $this->x[$preArrIdx]) && $this->y[$curArrIdx] == $this->y[$preArrIdx]) {
47:             $newX[$curArrIdx] = $r - 1;
48:         } else {
49:             $newX[$curArrIdx] = $r++;
50:         }
51:     }
52: 
53:     $this->x = $newX;
54: }

taxXtaxY 是用来做基数排序的桶
saX[i]=j 表示按第一关键字排序后,排名为第 i 名的字符串下标是 j
saY[i]=j 表示按第二关键字排序后,排名为第 i 名的字符串下标是 j

第 9~21 行,对 y 进行基数排序:

for ($i = 0; $i < $n; $i++) {
    $digit = $this->y[$i];
    $idx = $taxY[$digit]++;
    $saY[$idx] = $i;
}

第 23~26 行,对 x 进行基数排序,注意这时候,我们要依赖 y 数组的排序结果(下面的代码为什么这么写,要看明白,没看明白的,说明你没看懂 LSD 基数排序算法,反正我看了好久)

for ($i = 0; $i < $n; $i++) {
    $r = $saY[$i];
    $digit = $this->x[$r];
    $idx = $taxX[$digit]++;
    $saX[$idx] = $r;
}

然后根据 x 和 y 生成新的 x 数组

$newX = array_fill(0, 2 * $n, 0);

$r = 1;
$newX[$saX[0]] = $r++;

for ($i = 1; $i < $n; $i++) {
    $curArrIdx = $saX[$i];
    $preArrIdx = $saX[$i - 1];
    if (($this->x[$curArrIdx] == $this->x[$preArrIdx]) && $this->y[$curArrIdx] == $this->y[$preArrIdx]) {
        $newX[$curArrIdx] = $r - 1;
    } else {
        $newX[$curArrIdx] = $r++;
    }
}

$this->x = $newX;

完整的代码实现:https://github.com/sunznx/php-datastruct/blob/master/src/SuffixArray/SuffixArray.php