Suffix arrays are powerful tools for string processing and analysis. A suffix array is a data structure that allows efficient solution of various string-related problems.
There are multiple implementations of suffix arrays, with the doubling method achieving a time complexity of \(O(n \log n)\). While it's possible to implement suffix arrays in \(O(n)\) time, such implementations are considerably more complex.
Throughout this article (including code examples), we'll use 1-based indexing for strings. We'll denote suffix \(i\) as the substring of string \(S\) of length \(n\) starting at position \(i\) and extending to the end of the string.
Algorithmic Approach
Fundamental Concept
The core idea behind suffix arrays is to preprocess the ranking of each suffix and then perform at most \(O(\log n)\) sorting operations to arrange all suffixes in lexicographical order. This is achieved using radix sort, resulting in an overall time complexity of \(O(n \log n)\).
Let's define \(rk_i\) as the rank of suffix \(i\) and \(sa_i\) as the starting position of the suffix with rank \(i\). An important property is that initially, \(rk_{sa_i} = sa_{rk_i} = i\).
Our primary focus is updating the \(rk\) values using a doubling approach. Let \(rk_{w,i}\) represent the lexicographical rank of the substring of length \(w\) starting at position \(i\) among all substrings of length \(w\) in string \(S\).
Given that we've already computed \(rk_{w,i}\), we can determine the values for \(rk_{2w,i}\). For the substring of length \(2w\) starting at position \(i\), the first half has rank \(rk_{w,i}\) and the second half has rank \(rk_{w,i+w}\).
We sort all pairs \((rk_{w,i}, rk_{w,i+w})\) for \(1 \leq i \leq n\) (with \(rk_{i+w} = -\infty\) when \(i+w > n\)). The rank \(rk_{2w,i}\) is then simply the position of the pair \((rk_{w,i}, rk_{w,i+w})\) in this sorted list—essentially its discrete index. This allows us to sort all suffixes of \(S\) in lexicographical order in \(O(\log n)\) sorting steps. Initially, \(rk_i = s_i\).
Intuitively, setting \(rk_i = s_i\) initially sorts all suffixes by their first character. When multiple suffixes share the same first character, we need to sort them by their second character. The second character of suffix \(i\) is equivalent to the first character of suffix \(i+1\), whose rank we've already computed as \(rk_{1,i+1}\). Therefore, we can sort by \(rk_{1,i}\) as the primary key and \(rk_{1,i+1}\) as the secondary key, treating each suffix's first two characters as two separate components and calculating their combined rank.
This process continues for longer substrings. After the second sorting step, we sort all sufffixes by their first 4 characters. For suffix \(i\), we can split its first 4 characters into two substrings of length 2 starting at positions \(i\) and \(i+2\), then use the ranks obtained in the first sorting step to determine their combined rank.
By treating suffixes of length \(i\) as substrings of length \(2^n-1\) in an extended string \(S'\) formed by concatenating another copy of \(S\), we ensure that after \(i\) doubling steps, we can distinguish suffixes with identical first \(2^i\) characters. Since suffixes have different lengths, their lexicographical orders will differ, and thus \(O(\log n)\) sorting steps suffice to sort all suffixes.
Radix Sort
Radix sort is a stable sorting algorithm whose approach is similar to bucket sort. When sorting elements with \(k\) keys, we first sort by the \(k\)-th key, then by the \(k-1\)-th key, and so on, until finally sorting by the first key. Due to its stability, elements with the same current key value maintain they relative order from the previous sort.
Each sorting step resembles bucket sort. We first count the occurrences of each key value, compute prefix sums, then scan through the array, placing each element at the position indicated by the prefix sum of its key value, and decrementing that count.
Radix sort has a time complexity of \(O(n)\) and a space complexity proportional to the range of key values.
Reference Implementation (Basic Version)
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
const int MAXN = 1000005;
const int MAXV = 2000005;
int n, m;
vector<int> sa, id, rk, oldrk, cnt;
string s;
int main() {
// Read input string (1-indexed)
cin >> s;
s = " " + s; // Make string 1-indexed
n = s.length() - 1;
m = max(n, 300);
// Initialize vectors
sa.resize(n + 1);
id.resize(n + 1);
rk.resize(MAXV);
oldrk.resize(MAXV);
cnt.resize(m + 1);
// Initial ranking based on character ASCII values
for (int i = 1; i <= n; i++)
cnt[rk[i] = s[i]]++;
// Prefix sum for radix sort
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
// First suffix sort
for (int i = n; i >= 1; i--)
sa[cnt[rk[i]]--] = i;
// Doubling process
for (int j = 1; j < n; j <<= 1) {
// Second key radix sort
fill(cnt.begin(), cnt.end(), 0);
for (int i = 1; i <= n; i++)
id[i] = sa[i];
for (int i = 1; i <= n; i++)
cnt[rk[id[i] + j]]++;
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; i--)
sa[cnt[rk[id[i] + j]]--] = id[i];
// First key radix sort
fill(cnt.begin(), cnt.end(), 0);
for (int i = 1; i <= n; i++)
id[i] = sa[i];
for (int i = 1; i <= n; i++)
cnt[rk[id[i]]]++;
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; i--)
sa[cnt[rk[id[i]]]--] = id[i];
// Update ranks with discretization
oldrk = rk;
int p = 0;
for (int i = 1; i <= n; i++) {
if (oldrk[sa[i]] == oldrk[sa[i - 1]] &&
oldrk[sa[i] + j] == oldrk[sa[i - 1] + j])
rk[sa[i]] = p;
else
rk[sa[i]] = ++p;
}
}
// Output suffix array
for (int i = 1; i <= n; i++)
cout << sa[i] << " ";
cout << endl;
return 0;
}
Optimization Techniques
The basic doubling implementation of suffix arrays has a large constant factor. While it might pass on platforms like Luogu (problem P3809), it's likely to time out on more demanding platforms like LOJ#111. The following optimizations can significantly improve performance, reducing runtime by a factor of 4-5 in practice.
Radix Sort Optimization
We can optimize the second key radix sort since we only use it to count occurrences of each rank. The actual order of elements with the same second key doesn't matter. We can first record all positions \(i\) where \(i + j > n\), then iterate through the array and record \(sa_i - j > 0\). This can be implemented as:
int p = 0;
for (int i = n; i > n - j; i--)
id[++p] = i;
for (int i = 1; i <= n; i++)
if (sa[i] > j)
id[++p] = sa[i] - j;
Value Range Optimization
During discretization, we maintain a variable \(p\) representing the current discrete index. The rank values will always be in the range \([1, p]\). After each sort, we set \(m = p\) to optimize the value range for subsequent operations.
Non-contiguous Memory Access Optimization
The discretization process involves accessing non-contiguous memory. We can replace direct comparisons with a function:
bool compare_suffixes(int x, int y, int w) {
return oldrk[x] == oldrk[y] && oldrk[x + w] == oldrk[y + w];
}
Early Termination
If all suffixes have distinct ranks, the sorting is complete. At this point, the rank values will be in the range \([1, n]\). We can check if the maximum discrete index equals \(n\) and terminate early if true.
Reference Implementation (Optimized Version)
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstring>
using namespace std;
const int MAXN = 1000005;
const int MAXM = 2000005;
int n, m;
vector<int> rk, oldrk, sa, id, cnt;
string s;
bool compare_suffixes(int x, int y, int w) {
return oldrk[x] == oldrk[y] && oldrk[x + w] == oldrk[y + w];
}
int main() {
int p = 0;
cin >> s;
s = " " + s; // Make string 1-indexed
n = s.length() - 1;
m = 300; // Initial value range (ASCII values)
// Initialize vectors
rk.resize(MAXM);
oldrk.resize(MAXM);
sa.resize(n + 1);
id.resize(n + 1);
cnt.resize(m + 1);
// Initial ranking based on character ASCII values
for (int i = 1; i <= n; i++)
cnt[rk[i] = s[i]]++;
// Prefix sum for radix sort
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
// First suffix sort
for (int i = n; i >= 1; i--)
sa[cnt[rk[i]]--] = i;
// Doubling process with optimizations
for (int j = 1; ; j <<= 1) {
// Optimized second key sort
p = 0;
for (int i = n; i > n - j; i--)
id[++p] = i;
for (int i = 1; i <= n; i++)
if (sa[i] > j)
id[++p] = sa[i] - j;
// First key sort
fill(cnt.begin(), cnt.end(), 0);
for (int i = 1; i <= n; i++)
cnt[rk[id[i]]]++;
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; i--)
sa[cnt[rk[id[i]]]--] = id[i];
// Update ranks with discretization
oldrk = rk;
p = 0;
for (int i = 1; i <= n; i++) {
if (compare_suffixes(sa[i], sa[i - 1], j))
rk[sa[i]] = p;
else
rk[sa[i]] = ++p;
}
// Value range optimization and early termination
if (p == n) {
for (int i = 1; i <= n; i++)
sa[rk[i]] = i;
break;
}
m = p; // Update value range
}
// Output suffix array
for (int i = 1; i <= n; i++)
cout << sa[i] << " ";
cout << endl;
return 0;
}
Longest Common Prefix (LCP)
Basic Concepts
LCP (Longest Common Prefix) is an important concept in suffix array applications. In suffix array problems, \(LCP(i, j)\) is defined as the length of the longest common prefix between the suffixes \(sa_i\) and \(sa_j\). The LCP array and its properties enable efficient solutions to many string problems.
Here are some fundamental LCP properties (proofs are not included but can be found in specialized literature):
- \(LCP(i, j) = LCP(j, i)\)
- \(LCP(i, i) = \text{len}(sa_i) = n - i + 1\)
- \(LCP Lemma: \forall 1 \leq i \leq j \leq k \leq n, LCP(i, k) = \min\{LCP(i, j), LCP(j, k)\}\)
- \(LCP Theorem: \forall 1 \leq i < j \leq k, LCP(i, j) = \min\{LCP(j, j - 1)\}\)
The LCP Theorem can be derived from the LCP Lemma.
Example Problem and Solution
Problem: Given a string \(S\), find the length of the longest common prefix between any two substrings of \(S\). \(1 \leq |S| \leq 10^5\).
This problem can be solved using suffix arrays. First, consider a simplified version: finding the LCP between any two suffixes of \(S\).
For two suffixes starting at positions \(i\) and \(j\), we need to find \(LCP(rk_i, rk_j)\). Let's define \(height_i = LCP(i, i - 1)\). According to the LCP Theorem, the LCP between suffixes \(i\) and \(j\) is \(\min\{height_k | i < k \leq j\}\). This can be implemented efficiently using a Sparse Table (ST) in \(O(\log n)\) time per query.
To compute the \(height\) array, we can use the following property: \(\forall 1 < i \leq n, height_i \geq height_{i-1} - 1\). This allows us to compute the \(height\) array in \(O(n)\) time, similar to the approach used in the KMP algorithm or Manacher's algorithm. We start with \(height_{sa_1} = 0\), and for each subsequent position, we begin comparing from \(height_{i-1}\) and increment until we find a mismatch.
void compute_height() {
int k = 0;
for (int i = 1; i <= n; i++)
rk[sa[i]] = i;
for (int i = 1; i <= n; i++) {
if (rk[i] == 1) {
height[1] = 0;
continue;
}
if (k) k--;
int j = sa[rk[i] - 1];
while (i + k <= n && j + k <= n && s[i + k] == s[j + k])
k++;
height[rk[i]] = k;
}
}
Returning to the original problem, if we have substrings \([x, y]\) and \([l, r]\), their LCP length equals the LCP betweeen suffixes \(x\) and \(l\), but we must also consider the length constraints of the substrings. The final answer is the minimum of the LCP length and the lengths of both substrings.
Reference Implementation for LCP Problem
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 100005;
const int MAXM = 200005;
int n, q;
string s;
vector<int> sa, rk, oldrk, height, st, log_table;
// Suffix array computation (optimized version)
void build_suffix_array() {
int m = 300;
vector<int> cnt(m + 1), id(n + 1);
// Initial ranking
for (int i = 0; i < n; i++)
cnt[rk[i] = s[i]]++;
// Prefix sum
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
// First sort
for (int i = n - 1; i >= 0; i--)
sa[--cnt[rk[i]]] = i;
// Doubling
for (int j = 1, p = 0; j < n && p < n; j <<= 1, m = p, p = 0) {
// Second key sort
for (int i = n - j; i < n; i++)
id[p++] = i;
for (int i = 0; i < n; i++)
if (sa[i] >= j)
id[p++] = sa[i] - j;
// First key sort
fill(cnt.begin(), cnt.begin() + m + 1, 0);
for (int i = 0; i < n; i++)
cnt[rk[id[i]]]++;
for (int i = 1; i <= m; i++)
cnt[i] += cnt[i - 1];
for (int i = n - 1; i >= 0; i--)
sa[--cnt[rk[id[i]]]] = id[i];
// Update ranks
oldrk = rk;
p = 0;
rk[sa[0]] = p++;
for (int i = 1; i < n; i++) {
bool equal = (oldrk[sa[i]] == oldrk[sa[i - 1]] &&
oldrk[sa[i] + j] == oldrk[sa[i - 1] + j]);
rk[sa[i]] = equal ? p - 1 : p++;
}
}
}
// Compute height array
void compute_height() {
for (int i = 0; i < n; i++)
rk[sa[i]] = i;
height.resize(n);
int k = 0;
for (int i = 0; i < n; i++) {
if (rk[i] == 0) {
height[0] = 0;
continue;
}
if (k) k--;
int j = sa[rk[i] - 1];
while (i + k < n && j + k < n && s[i + k] == s[j + k])
k++;
height[rk[i]] = k;
}
}
// Build sparse table for RMQ
void build_sparse_table() {
int logn = log2(n) + 1;
st.resize(n * (logn + 1));
// Initialize for length 1
for (int i = 0; i < n; i++)
st[i] = height[i];
// Build table
for (int j = 1; j <= logn; j++)
for (int i = 0; i + (1 << j) <= n; i++)
st[i + (n * j)] = min(st[i + (n * (j - 1))],
st[i + (1 << (j - 1)) + (n * (j - 1))]);
// Precompute log2 values
log_table.resize(n + 1);
for (int i = 2; i <= n; i++)
log_table[i] = log_table[i / 2] + 1;
}
// Query LCP in range [l, r]
int query_lcp(int l, int r) {
if (l > r) swap(l, r);
if (l == r) return height[l];
int j = log_table[r - l];
return min(st[l + n * j], st[r - (1 << j) + 1 + n * j]);
}
int main() {
// Read input
cin >> s;
n = s.length();
// Build suffix array
build_suffix_array();
// Compute height array
compute_height();
// Build sparse table for RMQ
build_sparse_table();
// Process queries
int q;
cin >> q;
while (q--) {
int x, y, l, r;
cin >> x >> y >> l >> r;
x--; y--; l--; r--; // Convert to 0-based
// Special case when substrings start at same position
if (x == l) {
cout << min(y - x + 1, r - l + 1) << endl;
continue;
}
// Get ranks of starting positions
int rkx = rk[x], rkl = rk[l];
if (rkx > rkl) {
swap(x, l);
swap(y, r);
swap(rkx, rkl);
}
// Compute LCP
int lcp = query_lcp(rkx + 1, rkl);
int ans = min({lcp, y - x + 1, r - l + 1});
cout << ans << endl;
}
return 0;
}