Linear-time suffix sorting with SA-IS

By suffix sorting, we mean the problem of creating a suffix array for a string \(T\) of \(n\) characters from an alphabet \(\Sigma\). We write:

• \(T[i]\) for an individual character of \(T\), where \(0≤i < n\),
• \(T[i…j]\) for a substring “\(T[i]\,T[i+1]\,\cdots\,T[j]\)”, where \(0≤i≤j < n\), and
• \(T[i…]\) for a suffix \(T[i…n-1]\), where \(0≤i < n\).

The suffix array is a permutation of the integers \(0\) to \(n-1\) in which \(i\) comes before \(j\) if \(T[i…]\) < \(T[j…]\) when strings are compared lexicographically.

We assume that characters are integers \(0…|\Sigma|-1\), where \(0\) is a sentinel, required to appear at the end of \(T\) and only at the end. A valid \(T\) can be created from an input string by adding one to every character code and appending a zero sentinel like this:

const n = input.length+1;
const T = Array(n);
for (let i = 0; i < n-1; i++) { T[i] = input.charCodeAt(i)+1 }
T[n-1] = 0;

For instance, let \(T=\) “beebybeeybybaby$” where “$” is the sentinel (zero), and letters correspond to positive integers in the order of the alphabet. Then the suffix array is \((15,12,11,0,5,13,9,3,2,1,6,7,14,10,4,8)\) reflecting the lexicographic suffix order:

15: $
12: aby$
11: baby$
 0: beebybeeybybaby$
 5: beeybybaby$
13: by$
 9: bybaby$
 3: bybeeybybaby$
 2: ebybeeybybaby$
 1: eebybeeybybaby$
 6: eeybybaby$
 7: eybybaby$
14: y$
10: ybaby$
 4: ybeeybybaby$
 8: ybybaby$

Several recursive divide-and-conquer algorithms exist for sorting \(T\) in time \(O(n)\), under the assumption that \(|\Sigma|\) is \(O(n)\), without first constructing a suffix tree. The general scheme is:

1. Pick out some positions of \(T\) to create a string \(T'\) such that \(|T'|/|T|\) is (at most) a constant \(d<1\).
2. Recursively suffix sort \(T'\).
3. Induce the order of all the suffixes of \(T\) from the suffix order of \(T'\).

By making sure that steps 1 and 3 take altogether at most \(k·n\) time for some constant \(k\), we get a total time complexity of \(O(n)\), since the recursion formula \(\mathit{time}(n)≤kn+k'+\mathit{time}(dn)\) resolves to \(\mathit{time}(n)≤(k/(1-d))·n+k''\) with constants \(k'\) and \(k''\).

The algorithm by Kärkkäinen, Sanders, and Burkhardt has \(d=2/3\) and is reasonably simple. The algorithm by Kim, Sim, Park, and Park has \(d=1/2\) but a complex step 3. Ko and Aluru came up with a scheme of classifying of each position in \(T\) as either “S” or “L”, and picking the positions with one of these types (the one that appears the smallest number of times) for \(T'\). This makes \(d\) at most \(1/2\), usually smaller. Their idea was refined and simplified by Nong, Zhang, and Chan in two different ways. The variant they named SA-IS is the focus of this text.

We say that a position \(i\) of \(T\) is S-type if either \(T[i…] < T[i+1…]\) (the suffix starting at \(i\) is smaller than suffix starting at the following position) or \(i=n-1\) (the sentinel position). Otherwise (if the suffix starting at \(i\) is larger than suffix starting at the following position), we say that \(i\) is L-type. An S-suffix is a suffix that starts at an S-type position, and an L-suffix starts at an L-type position.

The following function computes the S/L-types of \(T\) in a boolean array isS, right-to-left in time \(O(n)\):

const computeTypeArray = (T) => {
    const isS = Array(T.length);
    isS[T.length-1] = true;                     // sentinel
    for (let i = T.length-1; i > 0; i--) {
        isS[i-1] = T[i-1] < T[i] || T[i-1] === T[i] && isS[i];
    }
    return isS;
}

A key observation is that if two suffixes start with the same character, but one is S and the other is L, then the S-suffix is larger than the L-suffix. This is because if an S-suffix starts with character \(c\), it must, after zero or more repetitions of \(c\), be followed by a character that is larger than \(c\). The opposite is true for an L-suffix: \(c\) must ultimately be followed by a smaller character. One (or both) of these larger/smaller characters determines the order of the two suffixes.

We say that:

• A position \(i\) is an LMS-position (short for leftmost S) if \(i > 0\), \(i\) is S-type, and \(i-1\) is L-type. It is simple to determine if a position is an LMS-position from the type array. Note that there are at most \(⌊n/2⌋\) LMS-positions.
• An LMS-suffix is a suffix starting at an LMS-position.
• An LMS-string is a substring that starts with one LMS-position and ends with the next LMS-position, except for the sentinel which is an LMS-string on its own.

The reduction to a smaller problem (to solve recursively) is done by treating the LMS-strings as basic blocks, as if they were large characters, and creating a shorter string \(T'\) whose characters correspond to the LMS-strings.

\(I\) is an array of length \(n\) which ultimately (when the construction algorithm is finished) holds the suffix array. Its values are suffixes of \(T\), represented by their starting positions \(0, …, n-1\). Saying that a suffix \(α\) is placed at position \(i\) of \(I\) is the same as saying that the integer \(I[i]\) is set to the starting position of \(α\). We also allow \(I\) to contain the value \(-1\) for “none”.

A bucket is a segment of \(I\) dedicated to suffixes with the same first character. There are as many buckets as there are different characters in \(T\), and the size of the \(c\)-bucket is the same as the number of times \(c\) occurs in \(T\).

Let the integer \(R\) be the current size of the alphabet (the maximum number of buckets). At the top recursion level, \(R=|\Sigma|\), but when the algorithm invokes itself recursively, the alphabet is different from the original one, and then \(R\) can be both smaller or larger than \(|\Sigma|\), although never larger than \(n\).

We use an array \(B\) of length \(R\) to hold one position for each bucket. The following function, which is effectively the first part of counting sort (also known as key-indexed counting), sets \(B\) to either the first position (if the end parameter is false) or just after last position (if end is true) of each bucket.

const setBucket = (T, B, end) => {
    B.fill(0);
    for (const c of T) { B[c]++ }
    for (let i = 0, sum = 0; i < B.length; i++) {
        B[i] = (sum += B[i]) - (end ? 0 : B[i]);
    }
}

Since the length of \(B\) is \(R\), both loops are executed at most \(O(n)\) times under the assumption that \(|\Sigma|\) is \(O(n)\), and the running time of the function is \(O(n)\).

The key trick of the SA-IS algorithm is inducing the order of all suffixes from the order of only the LMS-suffixes in the following way. Say that we know the order of the LMS-suffixes. We place them in order in \(I\), spread out at the ends of their buckets, and fill the rest of \(I\) with \(-1\). In our example, the positions of the LMS-suffixes in order are:

15: $
12: aby$
 5: beeybybaby$
 9: bybaby$
 3: bybeeybybaby$

They are placed in \(I\) like this:

Then we do:

1. Set bucket positions in \(B\) to the beginning of each bucket.
2. For each \(I[i]\) from left to right: if \(I[i]-1\) is an L-position, place it at the current position of its bucket, and advance the bucket position by one.
3. Set bucket positions in \(B\) to just past the end of each bucket.
4. For each \(I[i]\) from right to left: if \(I[i]-1\) is an S-position, decrease its bucket position by one, and place \(I[i]-1\) in that position.

In Javascript:

const induceOrder = (isS, I, T, B) => {
    setBucket(T, B, false);
    for (let i = 0; i < T.length; i++) {
        const j = I[i]-1;
        if (isS[j] === false) { I[B[T[j]]++] = j }
    }
    setBucket(T, B, true);
    for (let i = T.length-1; i >= 0; i--) {
        const j = I[i]-1;
        if (isS[j]) { I[--B[T[j]]] = j }
    }
}

Let us look at what effect step 2 has on our example. We start the scan with \(I[0]\) which is \(15\), and find that \(15-1=14\) is an L-position. Suffix \(T[14…]\) should go in the “y”-bucket (because \(T[14]\) is “y”) which starts at position \(12\), so we place \(14\) at \(I[12]\) and advance the “y”-bucket to \(13\). Then we look at \(I[1]\) which is \(12\), and \(12-1=11\) is also an L-position, so we place \(11\) at \(I[2]\), the beginning of the “b”-bucket (because \(T[11]\) is “b”), and advance the “b”-bucket to \(3\). Then we look at \(I[2]\), which we just set to \(11\), and again \(11-1=10\) is an L-position, and since \(T[10]\) is “y”, we place \(10\) at \(I[13]\), the next position of the “y”-bucket, and advance the “y”-bucket to \(14\). Then we look at \(I[3]\), which is \(-1\), and since \(-1-1=-2\) is not an L-position, we place nothing. And so on, until we get to the end of \(I\), which then contains:

Note that in addition to the LMS-positions, \(I\) now contains all the L-positions, in the correct suffix order. Also, they are grouped together at the beginning of their buckets, and since L-suffixes are smaller than S-suffixes that begin with the same character, they are actually in their correct final positions in the suffix array.

This may seem remarkable, but it is no coincidence. We now prove that this is always the case at the end of step 2.

To go from \(T\) to the reduced string \(T'\), we pick out the LMS-positions of \(T\), and translate the LMS-string in each of those positions to a character in a new alphabet. The order of the characters in the new alphabet needs to respect the order of the correpsonding LMS-strings. To find the right alphabet, we sort the LMS-strings, and then assign them numbers in increasing order.

A major simplification in the SA-IS algorithm compared to the algorithm of Ko and Aluru comes from the observation that the same induce procedure that we have already covered can be used for sorting the LMS-strings. We spread out the LMS-positions over the ends of their buckets, and fill the rest of \(I\) with \(-1\). (This is just like when setting up for inducing suffix order, except that now the LMS-strings are not sorted apart for the bucketing on the first characters). Then we run the induce procedure. At the end of step 2, L-suffixes are in place, sorted on their prefixes up to and including the first S-position. At the end of step 4, the S-suffixes have also been put in place, and the order between those of them that start at LMS-positions are in our desired LMS-string order.

We omit the proof, but again it follows the same pattern of reasoning as the correctness proof for inducing suffix order.

Recursive computation of suffix arrays can stop when all characters of \(T\) are distinct, because then suffix sorting only needs to consider the first character of each suffix, and we can find it simply as the inverse permutation of \(T\) in \(O(n)\) time like this:

const inversePermutation = (T) => {
    const I = Array(T.length);
    for (let i = 0; i < T.length; i++) { I[T[i]] = i }
    return I;
}

Now we are ready for the main suffix sorting function, which approximately follows the sample implementation of Nong, Zhang, and Chan. It takes two arguments: an array \(T\) of length \(n\) containing the string, and an integer \(R\) such that \(T[n-1]=0\) and \(1 ≤ T[i] < R\) for \(i=1,…,n-2\).

export const computeSuffixArray = (T, R) => {

First compute the type array and define a function to check whether a given \(i\) is an LMS-position. Allocate \(I\) and the bucket array \(B\).

    const isS = computeTypeArray(T);

    const isLMS = (i) => isS[i] && isS[i-1] === false;

    const I = Array(T.length);                  // position array
    const B = Array(R);                         // bucket array

Spread out the LMS-positions in \(I\), on the right sides of their buckets, with the rest of \(I\) filled with \(–1\). Then use the induce procedure to obtain the order of LMS-strings.

    setBucket(T, B, true);                      // find ends of buckets
    I.fill(-1);                                 // -1 is empty position
    let nʹ = 0;                                 // length of Tʹ
    for (let i = 1; i < T.length; i++) {
        if (isLMS(i)) {
            I[--B[T[i]]] = i;
            nʹ++;
        }
    }

    induceOrder(isS, I, T, B);                  // induce order of LMS-strings

Now the LMS-positions are in \(I\) in LMS-string order. Allocate \(T'\) and use it temporarily to hold the LMS-positions in correct order.

    const Tʹ = Array(nʹ);
    for (let i = 0, j = 0; i < T.length; i++) {
        if (isLMS(I[i])) { Tʹ[j++] = I[i] }
    }

Use \(I\) temporarily to map positions in \(T\) to LMS-string numbers, assigned in increasing order. We scan \(T'\) to find the LMS-positions. For each position, we check whether it points to an LMS-string that is different from its predecessor, in which case we increase the counter for distinct LMS-strings. We then place the character number of the current LMS-string (one less than the counter) in the position of the LMS-string in \(I\). This will let us use \(I\) for mapping the LMS-position to the corresponding character in \(T'\).

    I.fill(-1);                                 // set empty again
    I[Tʹ[0]] = 0;                               // assign 0 to sentinel
    let Rʹ = 1;                                 // counts assigned names
    for (let i = 1, p = Tʹ[0]; i < nʹ; i++) {
        const q = Tʹ[i];                        // the LMS-position to consider

        // Compare q to previous string p, by character and (if equal) S/L.
        for (let d = 0;; d++) {
            if (T[p+d] !== T[q+d] || isS[p+d] !== isS[q+d]) {
                Rʹ++;                           // p different from q
                p = q;                          // compare to this one instead
                break;                          // done with q
            }
            if (d > 0 && isLMS(p+d)) { break }  // LMS-string ends? done with q
        }
        I[q] = Rʹ-1;                            // map position to Tʹ character
    }

Finalize \(T'\) by getting the LMS-string characters from \(I\).

    for (let i = 0, j = 0; j < nʹ; i++) {
        if (I[i] >= 0) { Tʹ[j++] = I[i] }
    }

Compute \(I'\), the suffix array of \(T'\). If the length of \(T'\) is equal to the number of distinct LMS-strings, all characters in \(T'\) are distinct, and we can find \(I'\) as the inverse permutation of \(T'\). Otherwise, we invoke the suffix array computation recursively for \(T'\).

    const Iʹ = Rʹ < nʹ ? computeSuffixArray(Tʹ, Rʹ) : inversePermutation(Tʹ);

Reuse \(T'\) again for holding the LMS-positions of \(T\). Then use \(T'\) as a map to translate the contents of \(I'\) from \(T'\) positions to the corresponding \(T\) positions.

    for (let i = 0, j = 0; j < nʹ; i++) {
        if (isLMS(i)) { Tʹ[j++] = i }
    }
    for (let i = 0; i < nʹ; i++) { Iʹ[i] = Tʹ[Iʹ[i]] }

Spread out the sorted LMS-suffixes at the ends of buckets in \(I\). Then use the induce procedure to obtain the order of all suffixes, and return the result.

    setBucket(T, B, true);
    I.fill(-1);
    for (let i = nʹ-1; i >= 0; i--) {
        const j = Iʹ[i];
        I[--B[T[j]]] = j;
    }

    induceOrder(isS, I, T, B);

    return I;
}

It may be desirable to reduce the \(O(n)\) extra space used in the algorithm (for the S/L-type array, the bucket array, and for \(T'\) and \(I'\)). After all, a suffix tree takes \(O(n)\) space too, and the main reason for preferring a suffix array is that it requires less space.