Suffix arrays

Sort all the endings of a string, and many questions about its pieces become easy: where a pattern occurs, which piece repeats, how many different pieces there are. The sorted list of endings is the suffix array. With a second list, the lcp array, it does most of what a suffix tree does, in a fraction of the memory.

Sorted endings

A suffix of s is an ending of it: s[i:] for some place i. A string of n letters has n non-empty suffixes, one per place. The suffix array lists the places in the alphabetical order of their suffixes.

s = "banana"
# 5  a
# 3  ana
# 1  anana
# 0  banana
# 4  na
# 2  nana
# the suffix array is [5, 3, 1, 0, 4, 2]

In alphabetical order, a string comes before every longer string that starts with it: "a" < "ana". Python compares strings this way.

TRIE-4 mentioned the suffix tree, the compressed trie of all suffixes. The suffix array is the order of that tree’s leaves, read from left to right. It is a plain list of n numbers, while a suffix tree in Python needs many objects per letter.

The one-line way, sorted(range(n), key=lambda i: s[i:]), is correct, but each key is a copy of a suffix: n²/2 letters of memory in all, and comparing two long suffixes can take O(n). For a string like "aaaa…" of 100,000 letters that is 5 billion letters of copies. It is fine for a few thousand letters.

Prefix doubling

The faster method sorts by the first letter, then by the first 2 letters, then 4, 8, and so on. After each round, every suffix gets a rank: its place in the order so far, with equal ranks for suffixes whose first k letters are the same. A suffix shorter than k just uses all its letters.

Both halves already have ranks from the last round. So sorting by the pair (rank[i], rank[i + k]) sorts by the first 2k letters, without looking at a single letter. When i + k is past the end, use −1: nothing comes before everything. Step through the first variant.

# one round: rank[i] orders the suffixes by their first k letters
key = lambda i: (rank[i], rank[i + k] if i + k < n else -1)
sa.sort(key=key)                   # now sorted by their first 2k letters
new = [0] * n
for r in range(1, n):              # equal pairs, equal ranks
    new[sa[r]] = new[sa[r - 1]] + (key(sa[r - 1]) != key(sa[r]))
rank = new                         # then k doubles, until all ranks differ

Cost. After round t the order is right for the first 2^t letters, so after about log₂ n rounds (how many times you can double 1 before passing n) all ranks differ. On ordinary text it stops much sooner, once the first few letters already tell the suffixes apart. Each round here sorts with sort(), O(n log n), so the whole method is O(n log² n) (n times log n times log n). In Python this is the practical choice: the sort runs in C.

The ranks are whole numbers below n, so each round can use radix sort instead (SORT-6): a counting sort by the second number, then a stable counting sort by the first. That makes a round O(n) and the method O(n log n). Faster still are methods that build the array in O(n), such as SA-IS (2009); they are used in DNA tools and in compression.

Finding a pattern

Every place where a pattern P occurs is the start of a suffix that begins with P. In sorted order, those suffixes form one unbroken block. Binary search finds the block’s two ends (BSRC-3), comparing P with the first m letters of a suffix at each step: O(m log n) per search, after one sort.

from bisect import bisect_left, bisect_right

def places(s, sa, p):
    first = lambda i: s[i:i + len(p)]
    lo = bisect_left(sa, p, key=first)     # key= needs Python 3.10+
    hi = bisect_right(sa, p, key=first)
    return sorted(sa[lo:hi])

The lcp array and Kasai’s method

The lcp array stores, for each row of the suffix array after the first, the length of the longest common prefix (the shared start) of that suffix and the one in the row above. For banana: [–, 1, 3, 0, 0, 2]. Comparing each pair of rows from scratch costs O(n²) on "aaaa…".

Kasai’s method (2001) visits the suffixes in the order of the text, i = 0, 1, 2, …, not in sorted order, and keeps the last answer h.

Proof. Let j be the suffix in the row above i, sharing h letters. Drop the first letter of both: suffixes j + 1 and i + 1 share h − 1 letters, and j + 1 still comes first (they had the same first letter, so the rest decided the order). Every suffix between them in sorted order also starts with those h − 1 letters, because sorted order keeps strings with a common start together. In particular the row directly above i + 1 does.

# suffix i, in text order; h is carried over from suffix i - 1
j = sa[rank[i] - 1]                # the suffix in the row above
while i + h < n and j + h < n and s[i + h] == s[j + h]:
    h += 1
lcp[rank[i]] = h
if h > 0:
    h -= 1                         # the next answer is at least h - 1

Cost. h never goes above n and drops by at most 1 per suffix, at most n times in all. So it can rise at most 2n times, and the whole array costs O(n). Step through the second variant of the figure.

What the two arrays answer

• The longest piece that occurs twice: the largest lcp value, at the rows that have it.
• The number of different pieces: suffix sa[r] has n − sa[r] starts, and lcp[r] of them were already counted with the row above. So the count is n(n + 1)/2 minus the sum of the lcp array.
• The common start of any two suffixes: the smallest lcp value in the rows between them. With a sparse table over the lcp array (RNGQ-2), each such question costs O(1).
• The longest piece two strings share: build the arrays for a + "#" + b and take the largest lcp between two neighbouring rows where one suffix starts in a and the other in b.

Now try STRG-P5, Sorted endings.