Abstract
The Suffix Array is a classic text index enabling on-line pattern matching queries via simple binary search. The main drawback of the Suffix Array is that it takes linear space in the text’s length, even when the text itself is extremely compressible. Several works in the literature showed that the Suffix Array can be compressed, but they all rely on complex succinct data structures which in practice tend to exhibit poor cache locality and thus significantly slow down queries. In this paper, we propose a new simple and very efficient solution to this problem by presenting the Suffixient Array: a tiny subset of the Suffix Array sufficient to locate on-line one pattern occurrence (in general, all its Maximal Exact Matches) via binary search, provided that random access to the text is available. We prove that: (i) the Suffixient Array length \(\varvec{\chi }\) is a strong repetitiveness measure, (ii) unlike most existing repetition-aware indexes such as the \(\varvec{r}\)-index, our new index is efficient in the I/O model, and (iii) Suffixient Arrays can be computed in linear time and compressed working space. We show experimentally that, when using well-established compressed random access data structures on repetitive collections, the Suffixient Array \(\varvec{\textrm{sA}}\) is simultaneously (i) faster and orders of magnitude smaller than the Suffix Array and (ii) smaller and one to two orders of magnitude faster than the \(\varvec{r}\)-index. With an average pattern matching query time as low as 3.5 ns per character, our new index gets very close to the ultimate lower bound: the RAM throughput of our workstation (1.18 ns per character).
1 Introduction and Overview of the Paper
In this paper, we propose a new simple and remarkably efficient solution to the well-studied problem of indexing compressed text for fast pattern matching queries [3, 4]. This problem is motivated by the constantly-increasing amount of repetitive data produced in several applicative scenarios. Bioinformatics is an excellent example of such a phenomenon: with the quick advancement of DNA sequencing technologies in the last decade, an enormous amount of genomic data can now be produced efficiently in terms of both time and cost (for instance, since the introduction of Next-Generation Sequencing, the cost of sequencing a Human genome dropped from 1 million to just 600 dollars in roughly a decade [5]). These technological improvements opened the pangenome era, in which typical analyses often involve the comparison of billions of sequence fragments with large collections of reference genomes [6, 7]. This data explosion, however, did not correspond to an information explosion. Since genomes of the same species are nearly-identical (a typical human genome, for instance, differs by no more than 0.2% from the reference genome sequence [9]), a collection of same-species genomes is usually extremely compressible. The goal of a compressed index is to pre-process such a collection into a small (compressed) data structure in order to speed up subsequent pattern matching queries. Importantly, such queries must be supported directly on the compressed data without decompressing it.
Compressed Indexed Pattern Matching: State of the Art
Suffix arrays [10,11,12] are a classical (uncompressed) solution to the indexed pattern matching problem. They use linear space in the text’s length and are extremely simple: the Suffix Array \(\textrm{SA}\) is the lexicographically-sorted list of (the starting positions of) the text’s suffixes. Pattern matching on \(\textrm{SA}\) can be performed via simple binary search combined with text extraction, needed to compare the sought pattern with the text’s suffixes during binary search.
As mentioned above, however, linear space is not acceptable in modern big-data scenarios. Subsequent research therefore focused on techniques for compressing the Suffix Array and for supporting pattern matching on conceptually-different compressed text representations, such as Lempel-Ziv’77 (LZ77) [13] and grammar compression [14]. First research papers on compressed suffix arrays focused on entropy compression and exact pattern matching [15, 16]. Entropy, however, is a weak compressibility measure in the repetitive regime as it is based on symbol frequencies. As a consequence, the empirical entropy of a text T and a collection composed of, say, thousands of copies of T is the same (i.e., entropy is not sensitive to repetitions). Since then, a variety of indexes have been developed to capture and exploit the repetitiveness of the input data in different ways; see [3, 4] for an extensive survey on compressed indexes and repetitiveness measures. Notable works in this research direction include the run-length FM-index of Mäkinen and Navarro [17], which used space proportional to the number of runs r in the Burrows-Wheeler transform (BWT) of the text (where r is known to effectively capture repetitiveness in the text, see [3]) to count the number of occurrences of a pattern. Another key contribution is the r-index of Gagie et al. [18], which extended these functionalities to locate queries (that is, reporting the occurrences of the pattern in the indexed text) in the same asymptotic space. These indexes can be interpreted as techniques to compress the Suffix Array to space proportional to O(r) words. Compared with the Suffix Array, however, they suffer from a big issue: being based on backward searching the Burrows-Wheeler transform and on complex compressed data structures for rank, select, and predecessor queries, they can no longer be considered as simple data structures and are notoriously not cache efficient. As a matter of fact, searching a pattern of length m in those indexes has \(\Omega (m)\) I/O complexity (i.e., cache misses). In practice (as we also investigate in this article), this translates to the fact that these indexes are orders of magnitude slower than suffix arrays. This has recently been mitigated by techniques based on relative Lempel-Ziv (RLZ) compression of the Suffix Array [19], which however incur a larger space usage.
Extensions to Approximate Pattern Matching
While these works focused on exact pattern matching, in real applications one is usually interested in finding approximate occurrences of the pattern under some meaningful string distance (e.g. Hamming or edit distance). For instance, for approximately matching short reads under the edit distance, the tool b-move [20] offers a lossless, repetition-aware (run-length-compressed BWT) solution. Another widely used approach to approximate pattern matching, effective for various types of reads and commonly employed in practice [21], is to locate Maximal Exact Matches (MEMs), and then extend them to full (approximate) occurrences of the pattern. Given a pattern P, a MEM is a maximal (that is, not extendible either to the left nor to the right) substring P[i, j] occurring without errors in the collection. Bannai et al. [22] showed that one can augment the r-index to compute MEMs efficiently within the same O(r) space bound on top of a random-access oracle on the text. Their algorithm has two phases. In the first phase, it scans the pattern from right to left to compute for every suffix of the pattern the position in the text of an occurrence of the longest prefix of that suffix that occurs anywhere in the text. Then it scans the pattern again, but this time from left to right, to compute the length of a longest match of each of its suffixes using a random access oracle and the positions stored in the first phase, from which MEMs can be reported immediately. Subsequently, an efficient construction algorithm for the data structure of Bannai et al. [22] was presented by Rossi et al. [23], who also presented a practical implementation. The two-phases algorithm of Bannai et al. [22] was later simplified into a one-pass (left-to-right scan) algorithm by Boucher et al. [24] which allows for processing patterns in the streaming model.
1.1 Our Contributions
In this paper we present a new Suffix Array compression technique that is repetition-aware, simple, and remarkably efficient in practice. More in detail, we introduce the Suffixient Array (\(\textrm{sA}\)): a sampling scheme of the Prefix ArrayFootnote 1 (the co-lexicographically-sorted list of text prefixes) with the following properties.
-
1.
If random access is available on the text, then a sequence of simple binary searches on \(\textrm{sA}\) suffices to locate on-line one occurrence of the pattern P in the text (i.e., we return an occurrence of every prefix P[1, i] as soon as P[i] arrives) or, more in general, all its MEMs. Using state-of-the-art compressed random access techniques, we obtain a fully-compressed index.
-
2.
\(\textrm{sA}\) is compressed in the sense that its length \(\chi \) is upper-bounded by \(O(\bar{r})\), the number of equal-letter runs in the Burrows-Wheeler transform of the reversed text. Furthermore, \(\chi \) does not depend on the alphabet order (while \(\bar{r}\) does) and there exist infinite text families such that \(\chi \in o(\bar{r})\). This is confirmed in practice: on a repetitive collection containing variants of Human chromosome 19, our index is smaller than the toehold lemma of the r-index (i.e., its subset for locating one pattern occurrence — about half the size of the full r-index).
-
3.
Unlike most existing compressed indexes, pattern matching on our index is also efficient in the I/O model. This is confirmed in practice: our index is one to two orders of magnitude faster than the r-index [18]. We compare our query times with a hard lower bound — the time needed on our workstation to extract |P| contiguous characters from a random position in a large text (which approaches the RAM throughput as |P| increases) — and show that our index approaches this bound.
While our index can locate only one pattern occurrence, standard techniques [18, 25] can be used to list efficiently all other occurrences by enumerating Suffix Array values adjacent to the pattern occurrence. We do not enter in such details in this article and will describe this line of research in a future publication under preparation.
When the uncompressed text is used for random access, our technique can be directly compared with binary searching the Suffix (Prefix) Array. In order to obtain a self-index, however, we need to use a compressed text representation. A vast amount of literature has been devoted to supporting fast random access in compressed space, for example via Straight-Line programs [26], LZ77 compression [27], block trees [28], or relative Lempel-Ziv (RLZ) [29]. See [3] for a complete survey on the topic. Some of these data structures are also known to be very small and fast in practice [28, 29]. As a consequence, they can be used as a black-box in point (2) above, arguably yielding the first simple (and efficient, as we show experimentally) compressed self-index for repetitive sequences.
Suffixient Sets
To achieve this result, we start by introducing a new class of combinatorial objects of independent interest: suffixient sets. Let \(T\in \Sigma ^n\) be a text of length n. A substring \(\alpha = T[i,j-1]\) of T is said to be right-maximal if both \(\alpha \cdot a\) and \(\alpha \cdot b\) appear in T for some characters \(a \ne b\), or if \(\alpha \) is a suffix of T. Consider any set \(S \subseteq [n]\) such that the following holds: for every one-character right-extension T[i, j] of every right-maximal substring \(T[i,j-1]\) of T, there exists \(x\in S\) such that T[i, j] is a suffix of T[1, x]. We call a set S with this property suffixient (being sufficient to “capture” all right-extensions of right-maximal strings or, equivalently, “cover” paths starting at the root of the suffix tree of T and ending in the first character labeling every edge) and denote with \(\chi \) the size of a smallest (not necessarily unique) suffixient set for T.
We proceed by exhibiting a suffixient set of cardinality \(2\bar{r}\), where \(\bar{r}\) is the number of runs in the Burrows-Wheeler transform of T reversed. In particular, this implies \(\chi \le 2\bar{r}\), showing that \(\chi \) is a new meaningful repetitiveness measure. Observe that \(\chi \) is independent of the alphabet ordering, while \(\bar{r}\) is not. Since re-ordering the alphabet may reduce asymptotically \(\bar{r}\) [30], from the bound \(\chi \le 2\bar{r}\) we immediately obtain that there exist string families on which \(\chi = o(\bar{r})\). In other words, \(\chi \) is a better repetitiveness measure than \(\bar{r}\), even though we leave open the problem of determining whether \(\chi \) is reachable (i.e., if \(O(\chi )\) words are always sufficient to store the text). We get however close to this goal: we prove that any suffixient set is a string attractor [31], which implies that \(O(\chi \log (n/\chi ))\) words of space are sufficient to store T.
Suffixient Arrays
We define a Suffixient Array \(\textrm{sA}\) of T to be any (not necessarily unique) suffixient set S of smallest cardinality \(\chi \), sorted according to the co-lexicographic order of the text prefixes represented by S.
A Suffixient Array \(\textrm{sA}\) can be used to locate one occurrence of a pattern P[1, m] in T with the following simple strategy. Assume that P occurs in T and that the prefix P[1, j] (starting with \(j=0\)) is right-maximal in T. (i) Find (by binary search on \(\textrm{sA}\) and random access on T) a prefix T[1, x], for \(x\in \textrm{sA}\), suffixed by \(P[1,j+1]\). (ii) Right-extend the match by comparing \(T[x+1,\dots ]\) and \(P[j+2, \dots ]\). This way, we either match until the end of P or we find a mismatch \(T[x+1+t] \ne P[j+2+t]\), for some \(t\ge 0\). In the former case, we are done. In the latter case, we know that if \(P[1,j+1+t]\) occurs in T, then it is right-maximal in T so we can repeat steps (i) and (ii). Starting with the empty prefix of P, this procedure yields one occurrence of P in T and, as we show in this paper, can be easily extended to return all MEMs of P in T. While this procedure extracts from the text a number of characters proportional to \(|P|^2\) in the worst case, we show that it actually triggers at most \(O((1+m/B)\cdot d\log \chi )\) I/O operations (i.e., cache misses), where B is the I/O block size (i.e., cache line) and \(d \le m\) is the node depth of P in the suffix tree of T. In many practical applications (for example, DNA alignment), \(d \ll m\) holds (on uniform random texts, \(d \in O(\log _\sigma n)\) with high probability).
We furthermore show that the running time can be made linear in |P| also in the worst case by employing z-fast tries [32] and fast longest common prefix/suffix queries on T in compressed space.
Computing the Suffixient Array
We proceed with the following natural question: can the Suffixient Array \(\textrm{sA}\) be computed efficiently? To answer this question, we start by characterizing the smallest suffixient set. Intuitively, we show that a suffixient set is of smallest cardinality when it captures precisely all the \(\chi \) one-character extensions of all right-maximal substrings of T being also left-maximal. We call such substrings of T supermaximal extensions. We then describe four algorithms solving the problem, all based on the Burrows-Wheeler transform (\(\textrm{BWT}\)), the Longest Common Prefix array (\(\textrm{LCP}\)), and the Suffix Array (\(\textrm{SA}\)) of T reversed. The first algorithm runs in quadratic time and is conceptually very simple. The other three algorithms are optimizations of the first one. The second algorithm runs in \(O(n + \bar{r}\log \sigma )\) time and requires only one sequential pass on those three arrays. The third and fourth algorithms require random access on those arrays and run in optimal O(n) time. While being slower than the latter two when \(\bar{r} \in \omega (n/\log \sigma )\), the one-pass property of the second algorithm makes it ideal to be combined with Prefix Free Parsing (PFP) [33], a technique able to stream those three arrays in compressed working space. We show experimentally that this combination of PFP with our one-pass algorithm can quickly compute a smallest suffixient set in compressed working space on massive repetitive text collections.
We conclude the paper with experimental results testing Suffixient Array on pattern matching queries.
2 Preliminaries
We use the following notation: [i, j] for an interval indicating all values \(\{i,i+1,\dots ,j\}\) (if \(j\ell \). If \(\alpha [m]\ne T[x]\) for \(\forall x\in \textrm{sA}\), then return the pair (0, 0).
Observe that Operation \(\textsf{search}(\alpha )\) can be supported in time \(O(m\log \chi )\) by simple binary search on \(\textrm{sA}\) and random access on T, provided that random access on T can be performed in constant time per extracted character (otherwise, this time gets multiplied by time needed to access a single text character). This time can be sped up to \(O(m + \log \chi )\) by using classic longest-common-suffix (LCS) data structures on the Suffixient Array (in [10], the authors show how to achieve this on the Suffix Array using the symmetric LCP array), using \(O(\chi )\) additional words of space.
Algorithm 1 shows how to locate an occurrence of every P[1, i], where P is a query pattern, using the above two operations.
Observe that, when operation \(\textsf{search}()\) is implemented with binary search, Algorithm 1 works by performing a sequence of (at most) m simple binary searches. This is more than the single binary search of pattern matching on the Suffix (Prefix) Array, which however requires linear space in n (unless implemented with more complex compressed data structures such as the rlcsa [39]). Nevertheless, experimentally we will show that in practice binary searching \(\textrm{sA}\) is as fast as binary searching \(\textrm{SA}\).
Lemma 16
Algorithm 1 is correct and complete.
Proof
The claim is equivalent to the following:
Proposition 17
For every \(1\le i\le m\) the following hold. Every time Algorithm 1 reaches Line 8 with values i, j, P[1, i] is a suffix of T[1, j] so, in particular, P[1, i] occurs in T. Conversely, if P[1, i] occurs in T then Algorithm 1 will reach Line 8 with values i, j such that P[1, i] is a suffix of T[1, j].
We prove the claim inductively on i. At the beginning of execution, in Line 3 we set \(i\leftarrow j \leftarrow 1\). Assume P[i] occurs in T. If \(P[i] = T[j]\), then the if statement is not executed, we reach Line 8, and clearly Proposition 17 holds true. If, on the other hand, \(P[i] \ne T[j]\), then the call \(j\leftarrow \textsf{search}(P[1,i])\) at Line 5 will return a position j with \(P[i] = T[j]\) since P[i] is a one-character extension of a right maximal string (the empty string) and \(\textrm{sA}\) is a suffixient set. We conclude that, again, we reach Line 8 and Proposition 17 holds true. Assume now that P[i] does not occur in T. Then, it must be the case \(P[i] \ne T[j]\) so the if statement will be executed and the call \(j\leftarrow \textsf{search}(P[1,i])\) will set j to 0, causing the termination of the algorithm with message NOT_FOUND. Also in this case, Proposition 17 holds true.
Assume inductively that \(P[1,i-1]\), with \(i>1\), occurs in T (otherwise, the algorithm halted in a previous step) and that the claim holds for all prefixes of \(P[1,i-1]\). Let \((i-1,j-1)\) be the pair output by the algorithm on prefix \(P[1,i-1]\). We consider two cases. (a) P[1, i] occurs in T. Then, either (a.1) \(P[i] = T[j]\), in which case the if statement is not executed, we reach Line 8, and Proposition 17 holds true, or (a.2) \(P[i] \ne T[j]\). In this case, observe that \(P[1,i-1]\) is right-maximal in T since both P[1, i] and \(P[1,i-1]\cdot T[j]\) occur in T. Then, the call \(j\leftarrow \textsf{search}(P[1,i])\) at Line 5 will return a position j with \(P[1,i] = T[j-i+1,j]\) since P[1, i] is a one-character extension of a right maximal string and \(\textrm{sA}\) is a suffixient set. We conclude that, again, we reach Line 8 and Proposition 17 holds true. (b) P[1, i] does not occur in T. Then, necessarily \(P[i] \ne T[j]\), so the if statement will be executed and the call \(j\leftarrow \textsf{search}(P[1,i])\) will set j to 0, causing the termination of the algorithm with message NOT_FOUND. Also in this case, Proposition 17 holds true. \(\square \)
Algorithm 1 calls \(\textsf{search}(P[1,i])\) at most m times. We can say more: it is actually not hard to see that Algorithm 1 calls \(\textsf{search}(P[1,i])\) at most \(d \le m\) times, where d is the node depth of the locus of P in the suffix tree of T. If Algorithm 1 outputs NOT_FOUND, then d is 1 more than the node depth of the locus of the longest prefix of P that occurs in T. We conclude that, if \(\textsf{search}(P[1,i])\) is implemented via binary search on \(\textrm{sA}\) and if random access on T costs constant time per extracted character, then Algorithm 1 runs in \(O(m\cdot d\log \chi ) \subseteq O(m^2\log \chi )\) time. More in general we obtain:
Theorem 18
Let T[1, n] be a text stored in a data structure supporting the extraction of T[i], for any \(i\in [n]\), in time \(t_a\). Then, given a pattern P[1, m], using the Suffixient Array of T (\(\chi \) words) we can locate one occurrence of every prefix P[1, j] (\(j\in [m]\)) in total \(O(t_a\cdot m\cdot d \log \chi ) \subseteq O(t_a\cdot m^2 \log \chi )\) time, where \(d\le m\) is the node depth of the locus of P in the suffix tree of T.
The above query time can be reduced to \(O(t_a\cdot m\cdot d +d\log \chi )\) by using classic longest-common-suffix (LCS) data structures taking \(O(\chi )\) additional words of space (see [10]). For random strings the length of the longest repeated substring is \(O(\log _\sigma n)\) with high probability, which implies \(d\in O(\log _\sigma n)\). We leave to future experiments the study of d in practical applications. However, the main feature of our solution is that, differently from the FM-index [15] and the r-index [18], our index is cache-efficient, provided that the text is stored contiguously in memory. While this seems a hard requirement for a compressed random access oracle, in Section 8 we will show such an efficient oracle. We analyze this phenomenon in the I/O model where a block (e.g. a cache line) fits B characters:
Theorem 19
Let T[1, n] be a text stored contiguously in main memory. Then, given a pattern P[1, m], using the Suffixient Array of T (\(\chi \) words) we can locate one occurrence of P with \(O( (1+m/B) \cdot d \log \chi )\) I/O complexity, where \(d\le m\) is the node depth of the locus of P in the suffix tree of T.
Proof
As observed above, Algorithm 1 runs at most d binary searches, totaling \(O(d\log \chi )\) binary search steps. At each step we need to compare at most m consecutive characters of P and T, an operation having \(O(1+ m/B)\) I/O complexity. \(\square \)
In comparison, the FM-index [15] and the r-index [18] have \(\Omega (m)\) I/O complexity (i.e., every pattern character triggers a I/O operation). In Section 7 we will discuss heuristics that in practice tend to drastically reduce the impact of the term \(d \log \chi \) of Theorem 19, effectively bringing the I/O complexity very close to the optimum.
Remark 20
(Locating all occurrences) As mentioned in the introduction, with our mechanism it is actually possible to count and locate all pattern occurrences (rather than just one) efficiently in the I/O model using \(O(\bar{r})\) space on top of the random access oracle. This line of research, however, is out of the scope of this article and will be described in a future publication under preparation.
4.2 Finding Maximal Exact Matches
Algorithm 1 can be extended so that it returns all MEMs of P in T within the same time complexity. We present this procedure in Algorithm 2.
If P[i] does not occur in T, then Algorithm 2 simply does not output any triple of the form \((i,\cdot ,\cdot )\). Assuming some oracle for functions \(\textsf{search}(\cdot )\) and \(LCP(\cdot , \cdot )\), we prove:
Lemma 21
Given a pattern P[1, m], Algorithm 2 reports all the MEMs of P in T.
Proof
The claim is implied by the following invariant. Every time the condition of the while statement (Line 3) is evaluated, the following Properties hold:
-
1.
\(\ell = LCS(T[1,j-1],P[1,i-1])\),
-
2.
if \(\ell >0\) and \(i\le m\), then \(P[i] \ne T[j]\),
-
3.
\(LCS(T[1,j'-1],P[1,i-1]) \le \ell \) for every \(j' \in [n]\), and
-
4.
All MEMs \((i',j',\ell ')\) such that \(i' < i-1\) have been reported.
The invariant trivially holds the first time we enter in the while loop.
We proceed by induction. Assume that the invariant holds at Line 3 for some values of \(j,\ell \), and \(i\le m\). We want to show that it still holds after executing the operations in Lines 4-7. First observe that, by Property (1) of the invariant, \(\textsf{search}(P[i-\ell ,i]) = \textsf{search}(P[1,i])\) holds in Line 4. This is a first difference with Algorithm 1, which instead calls \(\textsf{search}(P[1,i])\) (we will need the variant of Line 4 in the next subsection).
Line does not modify \(i,j,\ell \), so it maintains the invariant valid. After running Line 4, by definition of \(\textsf{search}()\) we have that \(P[i-\ell '+1,i] = T[j'-\ell '+1,j']\) is the longest string (of length \(\ell '\)) suffixing a text prefix T[1, x] ending in a position \(x\in \textrm{sA}\) (in particular, \(x=j'\) holds in this case). Observe also that \(\ell '\le \ell +1\) because \(\ell '\) cannot exceed the length of \(P[i-\ell ,i]\), which is \(\ell +1\). Moreover, in Line 5 we do not report any MEM with \(\ell =0\), so we can assume that \(\ell >0\) to analyze which MEMs are reported in Line 5. We have therefore to only distinguish two cases in Line 5.
(Case A) If \(\ell ' = \ell +1\), then \((i-1,j-1,\ell )\) is not a MEM (because \(P[i-\ell ,i] = T[j'-\ell ,j']\), i.e., \(P[i-\ell ,i-1]\) can be extended to the right), and we correctly do not output it.
(Case B) If, on the other hand, \(\ell ' < \ell +1\) then we claim that \((i-1,j-1,\ell )\) is a MEM (and we correctly report it). To see this, observe that (i) \(P[i-\ell ,i-1]\) cannot be extended to the left (i.e., \(P[i-\ell -1,i-1]\) does not occur in T). This is implied by Property (3) of the invariant. (ii) \(P[i-\ell ,i-1]\) cannot be extended to the right (i.e., \(P[i-\ell ,i]\) does not occur in T). To see this, assume for a contradiction that \(P[i-\ell ,i]\) occurs in T. Let \(\alpha = P[i-\ell ,i-1]\). Then, property (2) of the invariant implies that \(\alpha \) is right maximal, since both \(\alpha \cdot P[i]\) and \(\alpha \cdot T[j]\) occur in the text (with \(P[i]\ne T[j]\)). But then, by definition of suffixient set, there must exist \(j''\in \textrm{sA}\) such that \(\alpha \cdot P[i] = P[i-\ell ,i]\) (of length \(\ell +1\)) suffixes \(T[1,j'']\). This contradicts the fact that \(\textsf{search}(P[i-\ell ,i])\) returned \((j',\ell ')\) with \(\ell '<\ell +1\). From (i) and (ii), we conclude that \((i-1,j-1,\ell )\) is a MEM so Line 5 correctly reports it.
Note that the above discussion implies that the algorithm correctly reported all MEMs ending before pattern position \(i-1\) included.
At this point, we prove that \(P[i-\ell '+1,i]\) cannot be extended to the left, i.e., that \(P[i-\ell ',i]\) does not occur in T. Assume, for a contradiction, that \(P[i-\ell ',i]\) (of length \(\ell '+1\)) occurs in T and let \(j''\) be such that \(P[i-\ell ',i] = T[j''-\ell ',j'']\). Recall that \(\ell ' \le \ell +1\) must hold. If \(\ell ' = \ell +1\), then \(LCS(T[1,j''-1], P[1,i-1]) \ge \ell ' = \ell +1\), contradicting Property (3) of the invariant. If \(\ell ' \le \ell \), then \(T[j''-\ell ',j''-1]\) (of length \(\ell '\)) is a suffix of \(T[j-\ell ,j-1]\), i.e., \(T[j''-\ell ',j''-1] = T[j-\ell ',j-1] = P[i-\ell ',i-1]\). But then \(\alpha = P[i-\ell ',i-1]\) is right-maximal: both \(\alpha \cdot P[i]\) and \(\alpha \cdot T[j]\) (with \(P[i]\ne T[j]\) by Property (2) of the invariant), of length \(\ell '+1\), occur in T. Again, by definition of suffixient set, this contradicts the fact that \(\textsf{search}(i-\ell , i)\) returned \((\ell ',j')\).
To conclude, let \(\Delta = LCP(P[i+1,m], T[j'+1,n])\) as per Line 6. By definition of LCP(): (i) Since above we proved that \(P[i-\ell '+1,i]\) cannot be extended to the left, it follows that the same holds for \(P[i-\ell '+1,i+\Delta ]\), (ii) \(P[i+\Delta +1] \ne T[j'+\Delta +1]\) (importantly, note that \(j'+\Delta +1 \le n\) because of the terminator $ at the end of T; \(i+\Delta +1\), on the other hand, can be equal to \(m+1\): in that case, we exit the while loop), (iii) \(P[i-\ell '+1,i+\Delta ] = T[j'-\ell '+1,j'+\Delta ]\), and (iv) \(P[i-\ell '+1,i+\delta ]\) does not correspond to a MEM for all \(0 \le \delta < \Delta \). Conditions (i–iv) immediately imply the invariant is still valid on \(i,j,\ell \) after the assignment at Line 7: Properties (1) and (3) of the invariant follow from (i) and (iii), Property (2) is equivalent to (ii), and property (4) follows from (iv) and from the fact that, as proved above, the algorithm correctly reported all MEMs ending before pattern position \(i-1\) included. \(\square \)
Implementing \(\textsf{search}(\cdot )\) by simple binary search on \(\textrm{sA}\) and random access on T, and \(LCP(\alpha ,T[i,n]) = \Delta \) with a sequence of \(\Delta +1\) random access queries on T, we obtain:
Theorem 22
Let T[1, n] be a text stored in a data structure supporting the extraction of T[i], for any \(i\in [n]\), in time \(t_a\). Then, given a pattern P[1, m], using the Suffixient Array of T (\(\chi \) words) we can find all MEMs of P in T in total \(O(t_a\cdot m^2\log \chi )\) time.
Proof
Variable i in Algorithm 2 increases at least by one unit at every iteration of the while loop, and the loop terminates when \(i > m\). In every iteration, we perform one call to \(\textsf{search}()\) (\(O(t_a\cdot m\log \chi )\) time with binary search on \(\textrm{sA}\) and random access on T), so overall all the calls to this operation cost \(O(t_a\cdot m^2\log \chi )\) time. Finally, operation \(LCP(P[i+1,m], T[j'+1,n]) = \Delta \) implemented via random access costs \(O(t_a\cdot \Delta )\) time. Since i is incremented by \(\Delta +1\) in each step (until surpassing value m), it follows that the total time spent on LCP() is \(O(t_a\cdot m)\). \(\square \)
4.3 Faster Queries
Using state-of-the-art techniques, our query times can be reduced from quadratic to linearithmic. Since a solution for MEMs also implies a solution for exact pattern matching, in this subsection we only discuss solutions for MEMs. We need two ingredients: z-fast tries [32] (in the revisited version of [40]Footnote 2) and Straight-Line Programs.
Lemma 23
(z-fast trie, [40, Thm. 6-7]) Given a set \(Z=\{\alpha _1,\cdots ,\alpha _q\}\) of q strings sorted in co-lexicographic order (\(i m\). In every iteration, we perform one call to \(\textsf{search}()\) (\(O(\log n)\) time with Lemma 26), so overall all the calls to this operation cost \(O(m\log n)\) time. Finally, operation \(LCP(P[i+1,m], T[j'+1,n]) = \Delta \) implemented via random access costs \(O(\Delta \log n)\) time using the structure of Lemma 24. Since i is incremented by \(\Delta +1\) in each step (until surpassing value m), it follows that the total time spent on LCP() is also \(O(m\log n)\).
We now prove correctness. If the structure of Lemma 26 always returns the correct result, then the algorithm is correct by Lemma 21. Note that condition (1) of the invariant used in the proof of Lemma 21 implies that, in Line 4 of Algorithm 2, \(P[i-\ell ,i-1]\) is always a substring of T. But then, letting \(\beta = P[i-\ell ,i]\) and \(k = |\beta |\), this means that all calls to \(\textsf{search}\) are of the form \(\textsf{search}(\beta )\), where \(\beta [1,k-1]\) appears in T. This is precisely the condition required by Lemma 26 in order for the structure to return correct results. This concludes the proof. \(\square \)
Remark 28
As far as the construction time is concerned, one can build the data structures described in this section in \(O(n \log n)\) expected time where n is the text length, provided that an SLP with g rules is given. Computing a smallest suffixient set of size \(\chi \) takes O(n) time as we will see later in Section 6. The z-fast trie in Lemma 23 for the prefixes ending at positions in the suffixient set can be constructed in O(n) expected time [40]. More specifically, one can build the enhanced suffix tree of the reversed text and Karp-Rabin fingerprints in O(n) time, then insert the subset of the suffixes corresponding to the suffixient set in \(O(\chi )\subseteq O(n)\) time with the computation of relevant information according to [40, Section 4]; if a collision is found, the procedure is repeated with a different hash function, and the expected number of repeats is O(1). The random-access data structure in Lemma 24 can be built in O(n) time [26, Theorem 1.1]. The construction time for the data structure in Lemma 25 is \(O(n\log n)\) expected time [41, Lemma 1] whose bottleneck is to find a collision-free Karp-Rabin hash function relying on [42, Section 7.2.1] and the rest of this data structure can be constructed in \(O(g)\subseteq O(n)\) time.
4.4 Even Faster Queries
We provide a faster solution in the case where pattern P is not chosen adversarially. With this term, we mean that substrings of P do not collide with substrings of T through the hash functions used by the structures of Lemmas 25 and 26 with high probability; in other words, we assume that an adversary did not craft P so that it generates such collisions. This is a realistic scenario in practice since a string chosen independently from (the random choices in) those data structures will collide with a text substring of T through those hash functions with low (inverse polynomial) probability only (see [40, 41]):
Remark 29
([40, 41]) The structures of Lemma 25 and Lemma 26 return the correct result with high probability, assuming that their inputs (respectively, P[i, j] and \(\beta \)) are not chosen by an adversary on the basis of the random seeds in the hash functions used by those data structures. To abbreviate, in such a case we will say that the inputs are “ not chosen adversarially ”.
Observe that the non-adversarial hypothesis is the same required in other randomized data structures such as classic hash tables and Bloom filters [43] (where the random choices precede the arrival of the input) and is common in other works using z-fast tries to index text [40].
Algorithm 2 can be understood as a traversal of \(d\le m\) suffix tree edges, where d is the number of times we would fully or partially descend edges in the suffix tree of T while finding MEMs (a generalization to the context of MEMs of parameter d used in Theorem 18). Figure 2 shows an example where d is much smaller than m.
By following the same analysis carried out in the proof of Theorem 27 but using the data structure in Lemma 25 for computing LCP(), we immediately obtain:
Theorem 30
Let T[1, n] be a text with Suffixient Array size \(\chi \) for which there exists a SLP of size g. Then Algorithm 2, when implemented with the structure of Lemma 26 to support operation \(\textsf{search}()\) and with the structure of Lemma 25 to support function LCP(), finds all MEMs of any pattern P[1, m] in T in \(O(m + d\log n)\) time using \(O(\chi + g)\) words of space, where \(d\le m\) is the number of times we would fully or partially descend edges in the suffix tree of T while finding those MEMs. The result is correct with high probability, provided that P is not chosen adversarially.
4.5 Offline Solutions with only Random Access and z-Fast Tries
Interestingly, fast random access to T and a z-fast trie on the set of strings \(Z = \{T[1,j]\,\ j\in \textrm{sA}\}\) are by themselves sufficient for fast offline pattern matching. To see why, assume P does occur in T and we already know the location of an occurrence \(T [j - i + 1, j]\) of P[1, i] in T. If \(T [j + 1] = P [i + 1]\) then \(T [j - i + 1, j + 1] = P[1, i + 1]\). Otherwise, querying the z-fast trie with \(P [1, i + 1]\) will return in \(O (\log (i + 1)) \subseteq O (\log m)\) time a position \(j' + 1\) such that \(T [j' - i + 1, j' + 1] = P [1, i + 1]\).
It follows that in \(O (t_a \cdot m + d \log m)\) time we can find a position \(j^*\) such that if P occurs in T then \(T [j^* - m + 1, j^*] = P\), which we can then check in \(O (t_a \cdot m)\) time. Here, \(d \le m\) is again the number of times we would fully or partially descend edges in the suffix tree of T while finding the MEMs of P with respect to T (so the node depth of the locus of P if P occurs in T). If P occurs in T then our time bound holds in the worst case; otherwise, it holds with high probability assuming an adversary has not chosen P on the basis of the hash function used in our z-fast trie. In any case, our solution works in \(O (m (t_a + \log m))\) time.
This solution is offline in the sense that we may not know until we have finished whether any prefix P[1, i] of P occurs in T. Moreover, if P does not occur in T then we do not learn which is the longest prefix of P that occurs in T.
Theorem 31
Suppose we already have \(t_a\)-time random access to a text T[1, n] with Suffixient Array size \(\chi \). Then we can store in \(O (\chi )\) space a z-fast trie such that in \(O (m (t_a + \log m))\) time we can check whether P[1, m] occurs in T and, if so, return the location of one occurrence.
We can adapt this idea to offline MEM-finding. Suppose T has a Suffixient Array \(\textrm{sA}\) of size \(\chi \), \(T^{\text {rev}}\) has a Suffixient Array \(\textrm{sA}'\) of size \(\chi '\), and we have fast random access to T (and thus also to \(T^{\text {rev}}\)) and z-fast tries for the sets of strings \(Z = \{T[1,j]\,\ j\in \textrm{sA}\}\) and \(Z' = \{T^{\text {rev}} [1, j]\,\ j \in \textrm{sA}'\}\). These z-fast tries take \(O (\chi + \chi ')\) total space.
Assume we already know the ending position j of a prefix T[1, j] of T with the longest common suffix with P[1, i] of any prefix of T, possibly without knowing the length of that common suffix. If \(T [j + 1] = P [i + 1]\) then \(T [1, j + 1]\) has the longest common suffix with \(P [1, i + 1]\) of any prefix of T. Otherwise, querying the z-fast trie for Z with \(P [1, i + 1]\) will return in \(O (\log (i + 1)) \subseteq O (\log m)\) time a position \(j' + 1\) such that \(T [1, j' + 1]\) has the longest common suffix with \(P [1, i + 1]\) of any prefix of T.
It follows that in \(O (t_a \cdot m + d \log m))\) time we can build an array J[1, m] such that T[1, J[i]] has the longest common suffix with P[1, i] of any prefix of T, for \(i \le m\), where \(d \le m\) is again the number of times we would fully or partially descend edges in the suffix tree of T while finding the MEMs of P with respect to T.
Because the prefixes of \(T^{\text {rev}}\) are the suffixes of T, in \(O (t_a \cdot m + d' \log m)\) time we can also build an array \(J' [1, m]\) such that \(T [J' [i], n]\) has the longest common prefix with P[i, m] of any suffix of T, where \(d'\) is the number of times we would fully or partially descend edges in the suffix tree of \(T^{\text {rev}}\) while finding the MEMs of \(P^{\text {rev}}\) with respect to \(T^{\text {rev}}\).
We can find the ending position \(e_1\) of the leftmost MEM \(P [1, e_1]\) in P in \(O (t_a \cdot e_1)\) time, by comparing the characters in P[1, m] and \(T [J' [1], n]\) until we find a mismatch at \(P [e_1 + 1]\). Since \(P [e_1 + 1]\) is in the second MEM \(P [s_2, e_2]\) in left-to-right order, we can find that MEM’s starting position \(s_2\) in \(O (t_a (e_1 - s_2 + 1))\) time by comparing the characters in \(P [1, e_1 + 1]\) and \(T [1, J [e_1 + 1]]\) from right to left until we find a mismatch at \(P [s_2 - 1]\).
Since the second MEM in P is the first MEM in \(P [s_2, m]\), we can repeat this whole process to find the ending position of the second MEM and the starting position of the third MEM, and so on. This step is similar to Li’s [44] forward-backward algorithm and takes time proportional to \(t_a\) times the total length of the MEMs (we could also use as a base an algorithm from [45] and [34, Section 2.6], which is faster in practice).
With high probability, our algorithm works correctly and in \(O (t_a \cdot m + (d + d') \log m)\) time plus time proportional to \(t_a\) times the total length of the MEMs, assuming an adversary has not chosen P on the basis of the hash function used in our z-fast tries.
Theorem 32
Suppose we already have \(t_a\)-time random access to a text T[1, n] with Suffixient Array size \(\chi \) whose reverse has Suffixient Array size \(\chi '\). Then we can store in \(O (\chi + \chi ')\) space two z-fast tries such such that with high probability we can find the MEMs of P[1, m] with respect to T in \(O (t_a \cdot m + (d + d') \log m))\) time plus time proportional to \(t_a\) times the total length of the MEMs, assuming an adversary has not chosen P on the basis of the hash function used in our z-fast trie.
5 Characterization of Smallest Suffixient Sets
We now move to the problems of computing Suffixient Arrays. We start in this section by characterizing smallest suffixient sets. We first define the concept of supermaximal extensions (Definition 33). Then, in Definition 34 we define a set (actually, a family of sets) \(\mathcal {S}\) containing all positions in T “capturing” all (and only the) supermaximal extensions of T. Finally, in Lemma 35 we prove that \(\mathcal {S}\) is a suffixient set of minimum cardinality.
Definition 33
We say that T[i, j] (with \(j\ge i\)) is a supermaximal extension if \(T[i,j-1]\) is right-maximal and, for each right-maximal \(T[i',j'-1] \ne T[i,j-1]\) (with \(i' \le j' \le n\)), T[i, j] is not a suffix of \(T[i',j']\).
The following definition depends on a particular total order \(<_t\) on [n]. We use \(<_t\) to break ties between equivalent candidate positions of the suffixient set, choosing the position of maximum rank according to \(<_t\) in case of a tie. In order to make notation lighter, we do not parameterize the set on the particular tie-breaking strategy (\(<_t\) will always be clear from the context).
Definition 34
Let \(<_t\) be any total order on [n]. We define a set \(\mathcal {S}\subseteq [n]\) as follows: \(x \in \mathcal {S}\) if and only if there exists a supermaximal extension T[i, j] such that (i) T[i, j] is a suffix of T[1, x], and (ii) for all prefixes T[1, y] suffixed by T[i, j], if \(y\ne x\) then \(y <_t x\).
We prove that \(\mathcal {S}\) is a suffixient set of minimum cardinality:
Lemma 35
For any tie-breaking strategy (total order) \(<_t\), \(\mathcal {S}\) is a suffixient set of minimum cardinality \(\chi \) for T.
Proof
To see that \(\mathcal {S}\) is suffixient, consider any right-maximal substring \(T[i,j-1]\) (\(i \le j \le n\)). We want to prove that there exists \(x \in \mathcal {S}\) such that T[i, j] is a suffix of T[1, x]. Let T[1, x] be a prefix of T being suffixed by T[i, j], breaking ties (on its endpoint x) by \(<_t\). If T[i, j] is a supermaximal extension, then by Definition 34 it holds \(x\in \mathcal {S}\) and we are done. Otherwise, let \(T[i',j'] \ne T[i,j]\) be a one-character extension of a right-maximal string \(T[i',j'-1]\) such that T[i, j] is a suffix of \(T[i',j']\).
Without loss of generality, let \(T[i',j']\) be the string of maximum length \(|T[i',j']| = j'-i' + 1\) with this property. Observe that there cannot be a right-extension \(T[i'',j''] \ne T[i',j']\) of a right-maximal string \(T[i'',j''-1]\) such that \(T[i',j']\) is a suffix of \(T[i'',j'']\), otherwise T[i, j] would be a suffix of \(T[i'',j'']\) and \(|T[i'',j'']|>|T[i',j']|\), contradicting the fact that \(T[i',j']\) is the longest such string. It follows that \(T[i',j']\) is a supermaximal extension. Let T[1, y] be a prefix of T being suffixed by \(T[i',j']\), breaking ties (on y) by \(<_t\). Then, by Definition 34, \(y\in \mathcal {S}\). Also in this case we are done, since T[i, j] suffixes \(T[i',j']\) and \(T[i',j']\) suffixes T[1, y], so by transitivity of the suffix relation, T[i, j] suffixes T[1, y].
To see that \(\mathcal {S}\) is of minimum cardinality, let \(\mathcal {S'}\) be a suffixient set. We prove \(|\mathcal {S}| \le |\mathcal {S'}|\) by exhibiting an injective function from \(\mathcal {S}\) to \(\mathcal {S'}\).
Let \(\textit{SE}\subseteq \Sigma ^+\) be the set of all supermaximal extensions. We preliminarily show that any prefix T[1, j] can be suffixed by at most one supermaximal extension: assume, for a contradiction, that T[1, j] is suffixed by two distinct supermaximal extensions \(T[i,j] \ne T[i',j]\). Without loss of generality, \(i' \ell \) and such that \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }-1]\) is right-maximal. Let \(\hat{\imath }= bwt(\hat{\jmath })\). Since \(T[j'-\ell ,j']\) is a suffix of \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }]\) we have that \(T[j'-\ell ,j'-1]\) is a suffix of \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }-1]\), hence it must be the case that \(\hat{\imath }\in box(i)\). Since \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }-1]\) is right-maximal and \(T[\hat{\jmath }] = c\), without loss of generality we can assume that \(\hat{\imath }\in \{k-1,k\}\), where k is a c-run break (such a c-run break must exist in box(i)) and \(\textrm{LCP}[k] \ge \hat{\ell }\). In particular, \(k \in box(i)\). Then, \(\textrm{LCP}[k] \ge \hat{\ell }> \ell \) and \(k\in box(i)\) yield \(\max \textrm{LCP}[B_{i,c}] \ge \hat{\ell }> \ell \), a contradiction.
(2) Let \(T[j'-\ell ,j']\) be a supermaximal extension and let \(c = T[j']\). Let \(i' = bwt(j')\). Since \(T[j'-\ell ,j'-1]\) is right-maximal, we can assume without loss of generality that either \(\textrm{BWT}[i'] \ne \textrm{BWT}[i'-1]\) or \(\textrm{BWT}[i'] \ne \textrm{BWT}[i'+1]\). Let therefore \(i' \in \{i,i-1\}\), where i is a c-run break. By definition of \(T[j'-\ell ,j']\), it holds \(\ell = \textrm{LCP}[i]\) (because \(T[j'-\ell ,j'-1]\) is the longest right-maximal string suffixing \(T[1,j'-1]\)). We need to prove that \(\ell = \max \textrm{LCP}[B_{i,c}]\). Assume, for a contradiction, that \(\ell < \hat{\ell }= \max \textrm{LCP}[B_{i,c}]\). This means that there exists a supermaximal extension \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }]\) with (i) \(T[\hat{\jmath }] = c\), (ii) \(\hat{\imath }\in \{k-1,k\}\), where \(\hat{\imath }= bwt(\hat{\jmath })\) and \(k \in B_{i,c}\) is a c-run break in box(i), and (iii) \(\textrm{LCP}[k] = \hat{\ell }\). But then, since \(k\in box(i)\) and \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }-1]\) is a right-maximal string of length \(\hat{\ell }> \ell \), we have that \(T[j'-\ell ,j'-1]\) is a (proper) suffix of \(T[\hat{\jmath }-\hat{\ell },\hat{\jmath }-1]\). This fact and \(T[\hat{\jmath }]=c\) contradict the fact that \(T[j'-\ell ,j']\) is a supermaximal extension. \(\square \)
We immediately obtain a simple quadratic algorithm computing a suffixient set of smallest cardinality: see Algorithm 3. Note that, in Lemma 41, if \(i''\in B_{i,c}\) and \(\textrm{LCP}[i''] = \textrm{LCP}[i]\), then \(box(i) = box(i'')\). Since there could be multiple c-run breaks \(i'' \in B_{i,c}\) such that \(\textrm{LCP}[i''] = \textrm{LCP}[i]\), if \(\textrm{LCP}[i] = \textrm{LCP}[i''] = \max \textrm{LCP}[B_{i,c}]\) then we have to break ties and insert only one of the corresponding text positions in the suffixient set. In Line 7 of Algorithm 3 we choose the largest such c-run break: as a matter of fact, this defines a particular tie-breaking strategy \(<_t\) between positions [n] (as required by Definition 34); i.e., for \(y,x\in [n]\) \(y <_t x\) iff \(T[1,y]<_{\text {colex}}T[1,x]\).
Computing \(\textrm{BWT}\), \(\textrm{LCP}\), and \(\textrm{SA}\) in Line 1 takes O(n) time. Additionally, for each \(i = 2, \dots , n\), Algorithm 3 computes the set \(B_{i,c}\) in line 7 by scanning \(\textrm{LCP}[\dots ,i-1, i, i+1, \dots ]\) in order to identify \(\textrm{LCP}[box(i)]\) (O(n) time for each such scan). It follows that the total running time is bounded by \(O(n^2)\). Correctness follows immediately from Lemma 41, Definition 34, and Lemma 35.
6.2 A One-Pass Algorithm
In this section, we speed up the simple quadratic algorithm provided in the previous section. We will refer to this version as “one-pass” since it only requires one scan of the \(\textrm{BWT}\), \(\textrm{SA}\), and \(\textrm{LCP}\) arrays for \(T^{\text {rev}}\). The algorithm is summarized in Algorithms 4 and 5. See below for a description.
Observe that Algorithm 3 runs in quadratic time due to the fact that it needs to scan \(\textrm{LCP}[box(i)]\) for each run break i, and those boxes overlap. The main intuition of Algorithm 4 is that instead of scanning box() for every position, we simply detect the end of the corresponding box for each suffixient set position. Suppose we already know a position \(j\in [n]\) such that j corresponds to a suffixient set position \(n-SA[j]+1\) to be reported. Note that by definition, \(\textrm{LCP}[j]\) is the maximum LCP value in \(B_{j,c}\) where \(c\in \Sigma \) is such that j is a c-break. Let \(\ell ,r\) be such that \(box(j)=[\ell ,r]\). Then by definition, it must hold that \(r=n\) or \(\textrm{LCP}[r+1]<\textrm{LCP}[j]\). Therefore, if we keep track of such a position j with a locally maximal LCP value for each \(c\in \Sigma \), we can determine whether its corresponding box box(j) ends at a particular position by simply comparing the LCP value. Note that it is not necessary to check this at every position. Instead, it is sufficient to check at the end of every BWT-run whether the corresponding box has already ended, by performing the comparison with the minimum LCP value within each BWT-run.
Let us describe Algorithm 4 in detail. Assume we scanned the arrays \(\textrm{BWT}\), \(\textrm{SA}\), and \(\textrm{LCP}\) up to position i. R is a map associating each \(c \in \Sigma \) with information related with some c-run break \(j \le i\) corresponding to a candidate supermaximal extension \(\alpha \cdot c\). Each entry \(R[c] = (len, pos, active)\) is composed of three values: \(R[c].len = \textrm{LCP}[j]\) (i.e., the length of the right-maximal string \(\alpha \)), \(R[c].pos = text(j') = n - \textrm{SA}[j'] +1\), where \(j'\in \{j-1,j\}\) is such that \(\textrm{BWT}[j']=c\) (i.e., the position of the last character of the supermaximal extension \(\alpha \cdot c\) in T), and a boolean flag R[c].active which is set to \(\textit{true}\) if and only if R[c].len is the maximum value in \(\textrm{LCP}[B_{j,c}\cap [i]]\) (breaking ties by smaller j: we keep the first encountered such maximum). Depending on the value of flag active, we will distinguish between active and inactive candidates.
Example 42
See Fig. 3 and assume we scanned \(\textrm{BWT}\) up to position \(i=4\). The candidate associated with letter \(\mathtt T\) is \(R[\texttt {T}] = (2,3,\textit{true})\), where \(R[\texttt {T}].len = \textrm{LCP}[4]=2\) and \(R[\texttt {T}].pos=n-SA[3]+1=3\). This candidate is active since \(\textrm{LCP}[4]\) is the largest value in \(\textrm{LCP}[B_{4,\texttt {T}}\) \(\cap [4]] = \{2\}\). When processing position \(i=6\), we find another \(\mathtt T\)-run break. Here we update \(R[\texttt {T}] = (\textrm{LCP}[i],n-SA[i']+1,\textit{true})=(5,6,\textit{true})\) where \(i'=i-1=5\) since \(6 \in B_{4,\texttt {T}} \cap [6]\), and \(5=\textrm{LCP}[6] > \textrm{LCP}[4]=2\).
Algorithm 4 works as follows. We scan the \(\textrm{BWT}\) left-to-right: \(\textrm{BWT}[2]\), \(\dots \), \(\textrm{BWT}[i]\) (for \(i=2, \dots , n\)). For \(c\in \Sigma \), let \(j R[c].len\), then \(i \in B_{j,c} \cap [i]\) and R[c].len is not a local maximum in \(\textrm{LCP}[B_{j,c}]\) (line 10). In this case, we have to replace the candidate stored in R[c] with the information associated with the new candidate i: letting \(i'\in \{i-1,i\}\) be such that \(\textrm{BWT}[i'] = c\), we replace \(R[c] \leftarrow (\textrm{LCP}[i], n-\textrm{SA}[i']+1, \textit{true})\) (line 11).
Example 43
Continuing Example 42. When processing the \(\mathtt T\)-run break \(i=7\), we see that \(\textrm{LCP}[7] < R[\texttt {T}].len\). Thus, we reached the end of \(box(6) = [6,6]\), where max \(\textrm{LCP}[B_{6,\texttt {T}}] = R[\texttt {T}].len = 5\), and insert \(R[\texttt {T}].pos\) in \(\mathcal {S}\). Next, we update \(R[\texttt {T}] = (4,0,\textit{false})\). On \(i=8\), we again update \(R[\texttt {T}] = (1,0,\textit{false})\). Finally, on \(i=9\), we get \(\textrm{LCP}[9] > R[\texttt {T}].len\) since box(6) and box(9) are now disjoint. We update \(R[\texttt {T}] = (\textrm{LCP}[9],11,\textit{true})\) to the active state.
We obtain:
Lemma 44
OnepassProof Given a text T[1, n] over alphabet of size \(\sigma \), Algorithm 4 computes a smallest suffixient set \(\mathcal {S}\) in \(O(n + \bar{r}\cdot \sigma )\) time and O(n) words of space.
Proof
Algorithm 4 scans the \(\textrm{BWT}\), \(\textrm{LCP}\), and \(\textrm{SA}\) exactly once and, for each run break, checks if any of the \(\sigma \) candidates is a supermaximal extension using Algorithm 5. Since Algorithm 5 performs a constant number of operations for each \(c \in \Sigma \), the whole algorithm runs in \(O(n + \bar{r}\cdot \sigma )\) time. In addition, the only supplementary data structure we need is R, which takes \(O(\sigma )\) words of space; thus if \(\sigma < n\), altogether, we take O(n) words of space.
We prove the correctness of Algorithm 4 by showing it computes the same output of Algorithm 3. In particular, given \(\mathcal {S}'\), the output of Algorithm 4, and \(\mathcal {S}\), the output of Algorithm 3, we show that given any value \(s = n - \textrm{SA}[i] + 1\), where \(i \in [n]\), (1) \(s \in \mathcal {S}\implies s \in \mathcal {S}'\) and (2) \(s \notin \mathcal {S}\implies s \notin \mathcal {S}'\).
- (1):
-
Let \((n - \textrm{SA}[i'] + 1) \in \mathcal {S}\), where i is a c-run break, such that \(i' \in \{i-1,i\}\) and \(\textrm{BWT}[i'] = c\). By Lemma 41, it follows that \(\max \textrm{LCP}[B_{i,c}] = \textrm{LCP}[i]\). Assume we scanned the \(\textrm{BWT}\) up to position i, there are three cases: \(i'\) is the position of the first occurrence of c in the \(\textrm{BWT}\), or there exists another c-run break at position \(j < i\), such that either \(i \in B_{j,c}\) and \(\textrm{LCP}[i] > \textrm{LCP}[j]\), or \(box(i) \cap box(j) = \emptyset \). For all these three cases, \(\textrm{LCP}[i] > R[c].len\); thus, i is set as the new active c candidate in R (lines 10-12). Now, if \(i'\) is the position of the last occurrence of c in the \(\textrm{BWT}\), R[c].pos is inserted in \(\mathcal {S}\) at the end of the algorithm (line 17); otherwise, let \(box(i) = [b,e]\), where \(b \le i \le e\), we have \(\forall s \in \textrm{LCP}[b,e], s \ge \textrm{LCP}[i]\) and no \(j' \in [i+1,b]\) is a c-run break such that \(\textrm{LCP}[j'] > \textrm{LCP}[i]\). Thus, R[c] is not updated until we scan \(\textrm{LCP}[e+1]\). At this point, \(R[c].len > \textrm{LCP}[e+1]\) and since R[c] is an active candidate; thus, Algorithm 4 insert \(R[c].pos = (n - \textrm{SA}[i] + 1)\) in \(\mathcal {S}'\) (Algorithm 5, line 4).
- (2):
-
Now let \((n - \textrm{SA}[i'] + 1) \notin \mathcal {S}\), it follows by Lemma 41, that \(\exists j \in B_{i,c}\) such that \(\textrm{LCP}[j] > \textrm{LCP}[i]\) (for simplicity we consider the smallest j). Assume we scanned the \(\textrm{BWT}\) up to position i; we need to consider two cases: (i) \(j > i\) and (ii) \(j < i\). (i) If \(j > i\), let \(R[c].len = \textrm{LCP}[i]\), then \(\forall l \in \textrm{LCP}[i+1,j], l \ge R[c].len\); thus, \((n - \textrm{SA}[i'] + 1)\) is not inserted in \(\mathcal {S}\) (lines 2-6, Algorithm 5) until position j. Here, we get \(\textrm{LCP}[j] > R[c].len\) and, update \(R[c].pos = (n - \textrm{SA}[j'] + 1)\) (lines 10-12), where \(j'\in \{j-1,j\}\) and \(\textrm{BWT}[j']=c\). Due to this, \((n - \textrm{SA}[i'] + 1)\) is dropped and not inserted in \(\mathcal {S}'\). (ii) If \(j < i\), then \(\exists l \in \textrm{LCP}[j+1,i], l < \textrm{LCP}[j]\); thus when we read l, R[c] is updated to the inactive state (line 6, Algorithm 5). However, since \(\forall s \in B_{i,c} \cap [j+1,i], \textrm{LCP}[s] < \textrm{LCP}[j]\) it means that R[c] is never updated to the active state for any position in \([j+1,i]\) (we skip lines 10-12), so also in this case \((n - \textrm{SA}[i'] + 1)\) is not inserted in \(\mathcal {S}'\).\(\square \)
Observe that in Algorithm 5 (procedure \(eval(\cdot )\)), only entries such that \(l < R[c].len\) are possibly modified. This suggests that R could be sorted in order to speed up operations. Indeed, it turns out that R can be replaced with a data structure based on balanced search trees, such that operation \(eval(\Sigma ,m,R,\mathcal {S})\) at Line 8 of Algorithm 4 costs \(O(\log \sigma )\) amortized time. With this modification, Algorithm 4 runs in \(O(n + \bar{r} \log \sigma )\) time.
Lemma 45
Algorithm 4 can be implemented with \(O(n + \bar{r}\log \sigma )\) running time and O(n) words of space where n is the text length, \(\bar{r}\) is the number of runs in \(\textrm{BWT}(T^{\text {rev}})\), and \(\sigma \) is the alphabet size.
Proof
We will prove that \(eval(\cdot )\) can be implemented with \(O(\log \sigma )\) amortized running time. We maintain a balanced binary search tree (BST) in which each node has key \(l=R[c].len\) for some \(c\in \Sigma \). Note that the BST has at most \(\sigma \) nodes, thus the height is \(O(\log \sigma )\). Each node stores such corresponding alphabet symbols c’s, maintaining them using two doubly linked lists, one for active symbols (i.e., those with \(R[c].active=\textit{true}\)), and the other for inactive symbols (i.e., those with \(R[c].active=\textit{false}\)). For each symbol \(c\in \Sigma \), we store R[c].pos and R[c].active, as well as the pointer to the BST node in which c is stored so that R[c].len can be retrieved by accessing c, and the pointer to the corresponding node of the linked list so that we can move symbols between different linked lists in O(1) time.
We implement Line 11 of Algorithm 4 as follows. We move the linked-list node corresponding to symbol \(\textrm{BWT}[i']\) into the linked list for active symbols of the BST node associated with l (after creating a new node if it does not exist) and update R[c].pos and R[c].active for \(c=\textrm{BWT}[i']\) accordingly. This takes \(O(\log \sigma )\) time.
Our optimization of Algorithm 5 then works as follows. We retrieve all nodes associated with lengths greater than l, which can be performed in \(O(\log \sigma + k_{ret})\) time where \(k_{ret}\) is the number of retrieved nodes. We iterate all active symbols \(c\in \Sigma \) stored in the retrieved nodes to add R[c].pos to \(\mathcal {S}\) in \(O(k_{act})\) time where \(k_{act}\) is the total number of active symbols to report. Then we concatenate all the linked lists in the retrieved nodes in \(O(k_{ret})\) time, and append it to the linked list for inactive symbols in the node with key l (we create a new node if it does not exist), which takes \(O(\log \sigma )\) time. Then we delete all the retrieved BST nodes (since they no longer contain any symbol) in \(O(k_{ret}\log \sigma )\) time. As a consequence, Algorithm 5 runs in \(O(k_{act}+(1+k_{ret})\log \sigma )\) time. Observe that at each c-run break, at most two symbols can become active and only O(1) new BST nodes can be created. Therefore, the sum of \(k_{ret}+k_{act}\) over all executions of \(eval(\cdot )\) is \(O(\bar{r})\). Therefore, the amortized running time is \(O(\log \sigma )\). \(\square \)
6.2.1 One-Pass Algorithm in Compressed Space
One important feature of our algorithm is that we only need one scan of the \(\textrm{BWT}\), \(\textrm{SA}\), and \(\textrm{LCP}\) to compute \(\mathcal {S}\). This means that our algorithm also works in the streaming scenario where these arrays are provided one element at a time, from first to last. This feature can be exploited to run the algorithm in compressed working space using Prefix-free Parsing [33]. This optimization does not provide strong theoretical bounds, but in practice (see Section 8) it drastically reduces the amount of time and working space needed to compute the Suffixient Array.
Prefix-free parsing (PFP) is a technique introduced by Boucher et al. [33] to ease the computational burden of computing the \(\textrm{BWT}\) of large and repetitive texts. Briefly, with one linear-time scan of the text T, PFP divides T into overlapping segments, called phrases, of variable length, which are then used to construct what is referred to as the dictionary D (i.e., the set of distinct phrases) and parse P of the text (i.e., the text T encoded as a sequence of phrases, represented as indexes in D). Then, with a separate linear-time algorithm, the \(\textrm{BWT}\) of T is directly computed from D and P; thus using space proportional to the combined size \(|D|+|P|\) of the two data structures. Very repetitive texts will tend to generate very small D, since text repetitions translate to repeated phrases. The size of P (number of phrases in which T is parsed), is instead controlled by a user-defined parameter and is a small fraction of T’s length n. In [23, 46], it was shown how to modify PFP in order to compute also the \(\textrm{SA}\) and the \(\textrm{LCP}\) arrays. This version streams the three arrays \(\textrm{BWT}\), \(\textrm{LCP}\), and \(\textrm{SA}\) of T in \(O(|D| + |P|)\) compressed space, from their first to last entry and in parallel (i.e., the triples \((\textrm{BWT}[i], \textrm{LCP}[i], \textrm{SA}[i])\) are output for \(i=1, 2, \dots , n\)); this is sufficient for running Algorithm 4 in compressed space, without affecting its running time. The only change consists in reading the text T backwards, in order to compute \(\textrm{BWT}(T^{\text {rev}})\), \(\textrm{LCP}(T^{\text {rev}})\), and \(\textrm{SA}(T^{\text {rev}})\). If the input T resides on disk, this can be achieved very easily with a backward scan of T (without increasing the number of I/Os - page swaps - with respect to reading T forward). If, instead, the text is read from a forward stream, then PFP yields \(\textrm{BWT}(T)\), \(\textrm{LCP}(T)\), and \(\textrm{SA}(T)\) and our algorithm computes a suffixient set for \(T^{\text {rev}}\). Our index still works, except that we have to reverse the query pattern before searching for it in the index. It is worth mentioning that it is actually also possible to modify PFP so that it streams \(\textrm{BWT}(T^{\text {rev}})\), \(\textrm{LCP}(T^{\text {rev}})\), and \(\textrm{SA}(T^{\text {rev}})\) while reading a forward stream for T. However, in view of the simpler solutions described above, we do not enter into such details.
In Section 8 we will show that our PFP-based optimization makes it possible to compute the Suffixient Array efficiently on massive repetitive datasets.
6.3 A Linear-Time Algorithm using LF-Mapping
In this section, we further speed up the one-pass algorithm provided in the previous section and achieve linear time. Our new algorithm is summarized in Algorithm 6.
As previously discussed, in Algorithm 4, for the one-character extension of every \(\textrm{BWT}\) equal-letter run \(\textrm{BWT}[i^*,i] = cc\dots ccx\) (with \(x\ne c\)) we run procedure \(eval \) (Algorithm 5) to check if \(\min \textrm{LCP}[i^*+1,i]\) drops below the current \(\textrm{LCP}\) maxima associated to y-run breaks, for all \(y\in \Sigma \). In the end, this step charges an additional \(O(\bar{r} \sigma )\) term (which we mentioned can be reduced to \(O(\bar{r} \log \sigma )\) by using opportune data structures). Intuitively, Algorithm 6 avoids this cost by calling the eval procedure only on the two candidates \(R[c']\), where \(c' \in \textrm{BWT}[i-1,i]\). We achieve this by introducing a new array \(\textrm{LF}[1,\sigma ]\) updated on-the-fly, implementing the LF-mapping property of the Burrows-Wheeler transform. More formally: assume we have scanned the \(\textrm{BWT}\) up to the c-run break i, and let \(i' \in \{i-1,i\}\) be such that \(\textrm{BWT}[i']=c\). Then, \(\textrm{LF}\) is such that \(\textrm{LF}[\textrm{BWT}[i']] = \textrm{LF}[c] = j\), with \(\textrm{SA}[j] = \textrm{SA}[i'] + 1\). In other words, \(\textrm{BWT}[j]\) is the character preceding \(\textrm{BWT}[i']\) in \(T^{\text {rev}}\). Our linear-time algorithm is based on the following idea. Assume for simplicity that i is not the first c-run break. Let \(i^* < i\) be the largest integer such that \(i^*\) is a c-run break as well (i.e., \(i^*\) is the c-run break immediately preceding i). As in Algorithm 4, our new Algorithm 6 maintains \(R[c].len = \textrm{LCP}[i^*]\). By Lemma 41, we have to discover if \(i\in B_{i^*,c}\); if this is the case, then we compare \(\textrm{LCP}[i]\) and \(\textrm{LCP}[i^*]\) and decide if \(\textrm{LCP}[i]\) is the new maximum in \(B_{i^*,c} \cap [i]\) (i.e., \(\textrm{LCP}[i] > \textrm{LCP}[i^*]\)) or if the local maximum among c-run breaks in \(B_{i^*,c} \cap [i]\) had already been found (i.e., \(\textrm{LCP}[i] \le \textrm{LCP}[i^*]\)). If, on the other hand, \(i\notin B_{i^*,c}\), then we insert in \(\mathcal {S}\) the local maximum R[c].pos relative to \(B_{i^*,c}\) if and only if \(R[c].active = \textit{true}\) (if \(R[c].active = \textit{false}\), then the suffixient position corresponding to the local maximum of \(B_{i^*,c}\) had already been stored in \(\mathcal {S}\)).
In order to discover if \(i\in B_{i^*,c}\), we distinguish two cases.
-
(i)
If \(\textrm{BWT}[i-1]=c\) then, since \(i^*\) is the previous c-run break, it must be the case that \(\textrm{BWT}[i^*-1, i^*, \dots , i-1, i] = y c c \dots c c x\), for some \(x\ne c\) and \(y\ne c\). It follows that \(i\in B_{i^*,c}\) if and only if \(\min \textrm{LCP}[i^*,i] \ge R[c].len = \textrm{LCP}[i^*]\). Algorithm 6 computes \(\min \textrm{LCP}[i^*,i]\) analogously to Algorithm 4, i.e., by updating a variable m storing the minimum \(\textrm{LCP}\) inside intervals corresponding to \(\textrm{BWT}\) equal-letter runs.
-
(ii)
If, on the other hand, \(\textrm{BWT}[i-1]=x \ne c\), then since \(i^*\) is the previous c-run break, it must be that \(\textrm{BWT}[i^*-1,i^*,\ldots ,i-1,i] = cy...xc\), for some \(x \ne c\) and \(y \ne c\) (and there are no other occurrences of c between \(\textrm{BWT}[i^*-1]\) and \(\textrm{BWT}[i]\)). As in the previous case, the goal is to compute \(\min \textrm{LCP}[i^*,i]\). We achieve this by using the \(\textrm{LF}\) array. Due to the way the array \(\textrm{LF}\) is defined and constructed, we know that \(\min \textrm{LCP}[i^*,i] = \textrm{LCP}[\textrm{LF}[c]]-1\). As a result, we have that \(i\in B_{i^*,c}\) if and only if \(\textrm{LCP}[\textrm{LF}[c]] - 1 \ge R[c].len = \textrm{LCP}[i^*]\).
From the above intuition, we obtain:
Lemma 46
Given a text T[1, n] , Algorithm 6 computes a smallest suffixient set \(\mathcal {S}\) in O(n) time and O(n) words of space.
Proof
Algorithm 6 scans the \(\textrm{BWT}\) exactly once and, for each position of it, performs O(1) accesses to the \(\textrm{SA}\) and the \(\textrm{LCP}\) and calls Algorithm 5 at most once for each run break. Since Algorithm 5 performs a constant number of operations and the LF-mapping array can be computed in O(n) time, the whole algorithm runs in O(n) time. In addition, the only supplementary data structure we need is R, which consumes \(O(\sigma )\) words of space; thus, altogether, we take \(O(n + \sigma )\) words of space. Since we assumed \(\sigma \le n\), the space consumption of Algorithm 6 is O(n) words.
We prove the correctness of Algorithm 6 by showing that it computes the same output as Algorithm 3. In particular, given \(\mathcal {S}'\), the output of Algorithm 6, and \(\mathcal {S}\), the output of Algorithm 3,
we show that for any value \(s = n - \textrm{SA}[i] + 1\), where \(i \in [n]\), (1) \(s \in \mathcal {S}\implies s \in \mathcal {S}'\) and (2) \(s \notin \mathcal {S}\implies s \notin \mathcal {S}'\).
- (1):
-
Let \((n - \textrm{SA}[i'] + 1) \in \mathcal {S}\), where i is a c-run break, such that \(i' \in \{i-1,i\}\) and \(\textrm{BWT}[i'] = c\). By Lemma 41, we have \(\max \textrm{LCP}[B_{i,c}] = \textrm{LCP}[i]\). Assume we have scanned the \(\textrm{BWT}\) up to position i. There are three cases: \(i'\) is the position of the first occurrence of c in the \(\textrm{BWT}\), or there exists another c-run break at position \(j < i\), such that either \(i \in B_{j,c}\) or \(box(i) \cap box(j) = \emptyset \). In all three cases, since \((n - \textrm{SA}[i'] + 1) \in \mathcal {S}\), we have \(\textrm{LCP}[i] > R[c].len\); thus, i is set as the new active c candidate in R (lines 17-19). If \(i'\) is the position of the last occurrence of \(\textrm{BWT}[i']\), since \(R[c].active = \textit{true}\), after traversing the whole \(\textrm{BWT}\), Algorithm 6 inserts \(R[c].pos = (n - \textrm{SA}[i] + 1)\) in \(\mathcal {S}'\) (line 4, Algorithm 5). Otherwise, let \(i < k'\) be the position of the next occurrence of c in the \(\textrm{BWT}\), such that \(k' \in \{k - 1, k\}\). Due to \(\max \textrm{LCP}[B_{i,c}] = \textrm{LCP}[i]\), it must hold \(\textrm{LCP}[k] < \textrm{LCP}[i]\) or \(k \not \in box(i)\). Here either, \(\textrm{BWT}[i-1,k'] = xcc...ccy\) where \(x \ne c\) and \(y \ne c\); thus, \(\exists l \in \textrm{LCP}[i+1,k']\) such that \(l < R[c].len\), or \(\textrm{LCP}[\textrm{LF}[c]] - 1 < R[c].len\). Again, since \(R[c].active = \textit{true}\), then Algorithm 6 insert \(R[c].pos = (n - \textrm{SA}[i] + 1)\) in \(\mathcal {S}'\) (line 4, Algorithm 5).
- (2):
-
Now let \((n - \textrm{SA}[i'] + 1) \notin \mathcal {S}\), this means that \(\exists j \in B_{i,c}\) such that \(\textrm{LCP}[j] > \textrm{LCP}[i]\) and j is a c-run break. We need to consider two cases: (i) \(j > i\) and (ii) \(j < i\). (i) If \(j > i\), \(\forall l \in \textrm{LCP}[i + 1, j]\), \( l \ge \textrm{LCP}[i]\), which implies \(\textrm{LCP}[\textrm{LF}[c]] - 1 \ge \textrm{LCP}[i]\) for all c-run breaks in \([i+1,j]\); thus, we never update R[c] to the inactive state (lines 2-6, Algorithm 5) until we scan position j. Here, we get \(\textrm{LCP}[j] > R[c].len\) and update \(R[c].pos = (n - \textrm{SA}[j'] + 1)\) (lines 17-19), where \(j' \in \{j - 1, j\}\) and \(\textrm{BWT}[j'] = c\). Due to this, \((n - \textrm{SA}[i'] + 1)\) is dropped and not inserted in \(\mathcal {S}'\). (ii) If \(j < i\), then \(\exists l \in [j+1,i]\) such that \(\textrm{LCP}[\textrm{LF}[c]]-1 < R[c].len\) or \(\textrm{LCP}[l] < \textrm{LCP}[i]\); thus, when we read l, R[c] is updated to the inactive state (line 6, Algorithm 5). However, since \(\forall s \in B_{i,c} \cap [j + 1, i]\), \( \textrm{LCP}[s] < \textrm{LCP}[j]\) it means that R[c] is never updated to the active state for any position in \([j + 1, i]\) (we skip lines 17-19), so also in this case \((n - \textrm{SA}[i'] + 1)\) is not inserted in \(\mathcal {S}'\).
\(\square \)
In Fig. 4 we show an example of how Algorithm 6 works.
6.4 Another Linear-Time Algorithm via Precomputing Boxes
In this section we propose a second linear algorithm which in practice is faster than Algorithm 6, at the cost of increasing by an additive term O(n) the space consumption. In Section 8 we show our experimental results comparing both algorithms.
For a given c-run break \(i_f\) for \(f\in [k]\), let \(i_f' \in \{i_f - 1, i_f\}\) be such that \(\textrm{BWT}[i_f'] = c\). Algorithms 4 and 6 will add to the suffixient set the position \(n - \textrm{SA}[i_f'] + 1\) if and only if
In the following, we expose a method to evaluate (1) in O(1) time for each \(i \in [n]\), which will lead to a linear-time algorithm. As a first step, we define first c-maximum positions (Definition 47). Later, we prove that those positions are exactly the positions satisfying (1) (Proposition 48). Finally, we show an algorithm computing a suffixient set of smallest cardinality by finding first c-maximum positions (Algorithm 7).
Definition 47
(first c-maximum) Let \(i_f\) the position of a c-run break. We say \(i_f\) is first c-candidate if one of the following conditions hold:
-
\(f = 1\), or
-
\(i_{f - 1} \le l_f - 1\)
Moreover, we say that \(i_f\) is first c-maximum if it is first c-candidate and one of the following conditions hold:
-
\(f = k\), or
-
for each \(f < j \le k\), if \(i_f \le l_{j} - 1\) then \(r_f + 1 < i_{j}\).
Intuitively, first c-candidate positions are those such that its \(\textrm{LCP}\) value is smaller than the \(\textrm{LCP}\) of the previous c-run break. On the other hand, first c-maximum positions are first c-candidate positions whose box only intersects c-run break positions not having larger \(\textrm{LCP}\) values. The following proposition establishes the equivalency between finding c-run breaks which are the leftmost local maxima within a box and finding first c-maximum positions.
Proposition 48
Let \(i_f\) be a c-run break. We have
if and only if \(i_f\) is first c-maximum.
Proof
Let \(B_{i_f, c} = \{i_{f - a}, \dots , i_f, \dots , i_{f + b}\}\). By definition it holds for every \(j\in [f-a,f+b]\) that \(\textrm{LCP}[i_j]\ge \textrm{LCP}[i_f]\), and it also holds that
where \([l_f,r_f]=box(i_f)\) and we define \(i_0:=0\) and \(i_{k+1}:=n+1\) for brevity.
\((\Rightarrow )\) Suppose
Then by the maximality of \(\textrm{LCP}[i_f]\), it holds that \(\textrm{LCP}[i_f]=\textrm{LCP}[i']\) for all \(i'\in {B_{i_f,c}}\). This implies that \(a = 0\) (due to the minimality of \(i_f\)) and \(box(i_f) = \cdots = box(i_{f + b})=[l_{f},r_{f}]\). Since it holds \(i_{f-1}\le l_f-1< l_f\le i_{f}\) by (2), \(i_{f}\) is first c-candidate. To prove that \(i_f\) is also first c-maximum, we need to prove that it holds for every \(f \textrm{LCP}[i_{j}]\) for every \(h\in [l_f,r_f]\), which implies that \(l_{j}\le l_{f}\le i_f\). Since \(l_{j}-1i_m\), it follows that \(i_f\) cannot be first c-maximum. On the other hand, if \(m < f\), then \(a \ge 1\). From \(l_f\le i_{f - a}\) in (2), it follows that \(l_f-1 < l_f\le i_{f-a}\le i_{f-1}\). Therefore, \(i_f\) cannot be first c-candidate. \(\square \)
As shown in the proof of Proposition 48, if \(i_f\) is first c-candidate and there exists a position \(i_f < i_e\) such that \(l_{f + 1} - 1< i_f, \dots , l_{e - 1} - 1 < i_f\), \(i_f \le l_e - 1\) and \(r_f + 1 < i_e\), then for any other position \(i_e < i_g\) such that \(i_f \le l_g - 1\) (if any exist) it holds \(r_f + 1 < i_g\). Indeed, \(i_e\) is also first c-candidate, which means that, to find first c-maximum positions, we just need to sequentially scan the \(\textrm{BWT}\) looking for first c-candidate positions and compare their boxes’ boundaries. Consider the two arrays defined based on the LCP array as the following.
Definition 49
(PSV/NSV arrays) For a given length-n integer array A[1, n], the previous smaller value array of A is an integer array of length n defined as \(\textrm{PSV}(A)[i] = \max (\{j \,|\, j< i, A[j] < A[i]\}\cup \{0\})\) for all \(i\in [n]\). In a similar way, the next smaller value array of A is defined as \(\textrm{NSV}(A)[i] = \min (\{j \,|\, j > i,A[j] < A[i]\}\cup \{n+1\})\) for all \(i\in [n]\).
Interestingly, if we have access to the \(\textrm{PSV}(\textrm{LCP})\) and \(\textrm{NSV}(\textrm{LCP})\) arrays, we can compute \(box(i) = [l_i, r_i] = [\textrm{PSV}(\textrm{LCP})[i] + 1, \textrm{NSV}(\textrm{LCP})[i] - 1]\) in O(1) time. These arrays can be computed in O(n) time using a stack-based algorithm [48] and one can also use space-efficient data structures (e.g., [47]) if a smaller working space is of interest. For simplicity, in the rest of this section \(\textrm{PSV}\), \(\textrm{NSV}\) denote \(\textrm{PSV}(\textrm{LCP})\), and \(\textrm{NSV}(\textrm{LCP})\), respectively.
Following the above ideas, we show in Algorithm 7 a procedure adding to the suffixient set first c-maximum positions. In the same spirit as Algorithms 4 and 6, R[c] keeps the information of the last c-run break found up to now which is first c-candidate. The algorithm sequentially scans the \(\textrm{BWT}\) looking for run breaks. For each c-run break at position i, if i is first c-candidate, R[c] is updated as \((sa\_pos \leftarrow i, text\_pos \leftarrow n - \textrm{SA}[i'] + 1, active \leftarrow \textit{true}, nsv \leftarrow \textrm{NSV}[i])\). Whenever R[c] is being updated, we check if the last first c-candidate is first c-maximum, and report it if so. We conclude this section with the following.
Lemma 50
Given a text T[1, n] over alphabet of size \(\sigma \), Algorithm 7 computes a smallest suffixient set \(\mathcal {S}\) in O(n) time and O(n) words of space.
Proof
We prove the correctness by showing the following invariant: After we scanned the \(\textrm{BWT}\) up to position i, it holds for every \(c\in \Sigma \) that (i) R[c] has the information of the last c-run break found up to position i which is first c-candidate unless c has not appeared yet, and (ii) all first c-maximum positions \(i'\) with \(\textrm{NSV}[i']0\)), and (iii) optimizes the extraction of long substrings by accessing the predecessor structure (storing phrase borders) only when needed, i.e., a number of times equal to the number of phrases overlapping with the string to be extracted.
Suffixient Array-based Indexes ( \(\textrm{sA}\) and opt- \(\textrm{sA}\) )
We implemented two variants of our suffixient-based index and combined them with the three random access oracles described above. The two variants implement the Optimizations described in Section 7. The first variant, labeled \(\textrm{sA}\) in the plots below, implements only Optimization I. The second variant, deemed opt-\(\textrm{sA}\), implements both Optimizations I and II. Note that we implemented these optimizations assuming DNA alphabet, i.e., \(\sigma = 4\), and we always set \(k = 14\).
Prefix Array index (PA)
This index implements classic binary search on the Prefix Array, combining it with the three random access oracles described above.
Toehold Lemma of the r -Index (toehold)
We used the subset of the r-index [18] (we employ the original C++ implementationFootnote 11) sufficient to locate one pattern occurrence (the so-called Toehold lemma); from the original implementation, we removed all data structures required to locate all other pattern occurrences. This made the index substantially smaller (compared to the full r-index). To ease reproducibility of our experiments, we included this subset of the r-index in the repository containing our Suffixient Array index.
8.4.2 Datasets and Patterns Extraction
In this experiment, we used the three genomic datasets from Table 1: Chromosome 19, Salmonella, and SARS-CoV-2. Since our opt-\(\textrm{sA}\) implementation is designed for DNA alphabets of four characters, we preprocess the datasets to remove all non-DNA characters. This ensures that all four competing methods use the same three input texts. For each dataset, we generated three sets of 100,000 patterns of lengths 10, 100, and 1000 each. Specifically, given a text T[1, n] and a pattern length m, we selected patterns by drawing uniform random positions \(i \in [1,n-m+1]\).
8.4.3 Results
The plots below show a comparison between the query times of the four indexes and the RAM throughput of our workstation when extracting substrings of length \(m = \) 10, 100, and 1000 from uniformly-chosen positions in a uniform text of 1 billion characters. On our workstation, we obtained the following RAM throughput benchmarks: 6.23341 ns/character for \(m=10\), 2.31762 ns/character for \(m=100\), and 1.41736 ns/character for \(m=1000\). The RAM throughput for reading the entire text of 1 billion characters was of 1.17591 ns/character. Note that the RAM throughput indicates a hard lower bound for processing a pattern of a certain length since it provides the maximum rate at which a block of length m can be retrieved from our system’s internal memory. In turn, this is a hard lower bound for any deterministic pattern matching algorithm guaranteeing a correct answer.
In all experiments, the lz77 oracle was always orders of magnitude slower than the rlz oracle, while always using essentially the same space. Symmetrically, the rlz oracle was always as fast as the bitpacked text, while using orders of magnitude less space. Since rlz always dominated the other two oracles, in Fig. 7 we only show results using the rlz oracle.
As the plots show, the optimization opt-\(\textrm{sA}\) only slightly increases the space usage of \(\textrm{sA}\) while dramatically speeding up query times (even by one order of magnitude in some cases). Our results show that opt-\(\textrm{sA}\) always dominates the Prefix Array by a wide margin: opt-\(\textrm{sA}\) was always faster than the Prefix Array while using orders of magnitude less space. As expected from our analysis of opt-\(\textrm{sA}\) in the I/O model (Theorem 19), opt-\(\textrm{sA}\) dominates the toehold lemma of the r-index by a wide margin on the moderately-repetitive Chr 19 dataset. On this dataset, opt-\(\textrm{sA}\) was always smaller than the toehold lemma and about half the size of the full r-index. At the same time, opt-\(\textrm{sA}\) was always from one to two orders of magnitude faster than the toehold lemma. On the more repetitive Salmonella and SarsCov2 datasets, opt-\(\textrm{sA}\) was only slightly larger than the toehold lemma and always smaller than the full r-index while at the same time being always from one to two orders of magnitude faster.
The plots show that opt-\(\textrm{sA}\) was always at most 1 order of magnitude slower than the RAM throughput. This is the case of Chrom. 19 and patterns of length 100. On the other hand, for patterns of length 10 and 1000 the results were much more favorable. In the best case, for SarsCoV2 and length 1000, opt-\(\textrm{sA}\) (about 3.5 ns per character) was only 2.5 times slower than the RAM throughput (about 1.4 ns per character). This aligns with our analysis in the I/O model, as for long patterns and repetitive collections made up of variations of a nearly-uniform string, fewer binary search steps are expected. We remark that our implementation is not yet heavily optimized and we expect further improvements with more optimization work.
9 Future Work
We list some open problems on which we are currently working.
Queries in \(O(\chi )\) Space
We have already mentioned the problem of determining whether \(\chi \) is a reachable measure, i.e., \(O(\chi )\) words are sufficient to store the text and to perform efficient random access. A possible way to attack this problem is to observe that our random access data structure based on string attractors (using \(O(\chi \log (n/\chi ))\) space) does not exploit any structural property of suffixient sets other than them being a string attractor. Suffixient sets are, however, very particular string attractors. Can we use other combinatorial properties of suffixient sets to support random access in \(o(\chi \log (n/\chi ))\) space?
Another interesting problem is to support counting and locating all pattern occurrences (rather than just one) in space \(O(\chi )\) on top of the random access oracle. Eventually, we would also like to compute the BWT range of the pattern with our mechanism, therefore completely replacing the Burrows-Wheeler transform. Locating could be achieved, for example, by encoding the \(\phi \) function of the r-index [18] in \(O(\chi )\) space. As we briefly mentioned in the paper, this is possible in space \(O(\bar{r})\) on top of the random access oracle (while retaining I/O efficiency), but \(\bar{r}\) can be asymptotically larger than \(\chi \).
The Repetitiveness Measure \(\chi \)
The relation \(\chi \le \bar{r}\) was always true in our experiments, even though our theory only predicts \(\chi \le 2\bar{r}\) (Lemma 10). Can we prove better upper-bounds for \(\chi \) as a function of \(\bar{r}\) or, more in general, as a function of other repetitiveness measures [3]?
Recently, [50] has shown that \(\chi \le 2r\) where r is the number of equal-letter runs in the BWT of the text. Since it is known that \(\bar{r}=O(r \lg ^2 n)\) [51], it would be interesting to see if there is such a relationship between the size of a smallest suffixient set of a string and that of the reversed.
Necessary Samples
While Suffixient Arrays are the optimal solution to the problem of minimizing the cardinality of a suffixient set (Definition 9), they are not necessarily the smallest sampling of the Prefix Array guaranteeing the correctness of Algorithm 1. In a sense, Suffixient Arrays guarantee a sufficient but not necessary condition for Algorithm 1 to work correctly. While we have determined what this necessary condition is, efficient algorithms to optimize it are still under investigation.
Optimizations
Our implementation is not yet heavily optimized: by using more cache-efficient data structures (in particular, Elias-Fano dictionaries), we believe we can get even closer to the RAM throughput of our machine. We will explore this direction in a future practical implementation of our index in the context of a usable DNA aligner for repetitive collections.
Approximate Pattern Matching
We are interested in whether suffixient sets can be adapted for use with the technique of finding all substrings of T within edit distance k of P by cutting P into \(k + 1\) pieces and, for each piece, trying to match that piece exactly and then match the other pieces allowing for a total of k edits to them [52]. More generally, this technique is referred to as search schemes [53]. If we can overcome the challenges then we believe suffixient sets are well-suited to this task since, if P is reasonably long and k is reasonably small, after we have exactly matched one piece we will be fairly deep in the suffix tree of T, where edge labels tend to be long and it is advantageous to descend edges without checking all the characters in their labels (which we do via \(\textrm{LCP}\) queries between suffixes of P and T).
The problem we face is determining all the ways to approximately match the other pieces. This can be done with backtracking in the suffix tree, but with our current implementation we do not see how to perform this backtracking: when we compute the \(\textrm{LCP}\) of a suffix of P and any suffix of T, we may descend several edges at once without realizing it and, for backtracking, later we should go back and partly explore the subtrees hanging off the path we descended. We are currently looking for ways to detect when we descend past a node in the suffix tree and find its string depth, so we can return to it later and visit its other children.
Data Availability
We provide all the necessary information to access the datasets used in our experimental analysis in the GitHub repository of our project.
Notes
Equivalently, the Suffix Array of the reversed text; for a formal definition of Prefix Array, see Definition 4. Of course, one can define the same concepts symmetrically as a sampling of the Suffix Array of the original text; we decided to sample the Prefix (rather than Suffix) Array in order to support on-line pattern matching (read next).
Although [40, Thm. 6] states a query time of \(O(\log m + \sigma )\), below the theorem’s statement the authors briefly observe that, within the same asymptotic space, query time can be reduced to \(O(\log m)\).
While in [41, Lemma 1] they only discuss LCP(), the symmetric function LCS() can be supported by simply reversing the text, an operation that does not change the SLP size g.
As a matter of fact, only one binary search is really useful: the one for the longest pattern prefix P[1, i] suffixing T[1, x] for some \(x \in \textrm{sA}\). Since, however, we do not know i in advance, our algorithm needs to run several binary searches on different pattern prefixes before finding i.
References
Cenzato, D., Olivares, F., Prezza, N.: On computing the smallest suffixient set. In: Proceedings of the 31st International Symposium on String Processing and Information Retrieval, pp. 73–87 (2024). https://doi.org/10.1007/978-3-031-72200-4_6
Depuydt, L., Gagie, T., Langmead, B., Manzini, G., Prezza, N.: Suffixient Sets. arxiv:2312.01359 (2024)
Navarro, G.: Indexing highly repetitive string collections, part I: Repetitiveness measures. ACM Comput. Surv. 54(2), 29–12931 (2022). https://doi.org/10.1145/3434399
Navarro, G.: Indexing highly repetitive string collections, Part II: Compressed indexes. ACM Comput. Surv. 54(2), 26–12631 (2022). https://doi.org/10.1145/3432999
DNA Sequencing Costs: Data. https://www.genome.gov/about-genomics/fact-sheets/DNA-Sequencing-Costs-Data Accessed: (07 Jan 2025)
Medini, D., Donati, C., Tettelin, H., Masignani, V., Rappuoli, R.: The microbial pan-genome. Curr. Opin. Genetics Dev. 15(6), 589–594 (2005). https://doi.org/10.1016/j.gde.2005.09.006
Wang, T., Antonacci-Fulton, L., Howe, K., Lawson, H.A., Lucas, J.K., Phillippy, A.M., Popejoy, A.B., Asri, M., Carson, C., Chaisson, M.J., et al.: The human pangenome project: A global resource to map genomic diversity. Nature 604(7906), 437–446 (2022). https://doi.org/10.1038/s41586-022-04601-8
Altshuler, D., Pollara, V.J., Cowles, C.R., Etten, W.J.V., Baldwin, J., Linton, L., Lander, E.S.: An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature 407, 513–516 (2000). https://doi.org/10.1038/35035083
Consortium, T.G.P.: A global reference for human genetic variation. Nature 526(7571), 68–74 (2015). https://doi.org/10.1038/nature15393
Manber, U., Myers, G.: Suffix arrays: A new method for on-line string searches. In: First Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 319–327. ACM (1990)
Baeza-Yates, R.A., Gonnet, G.H.: A new approach to text searching. In: Proceedings of the 12th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval. SIGIR ’89, pp. 168–175. ACM, New York, NY, USA (1989). https://doi.org/10.1145/75334.75352
Gonnet, G.H., Baeza-Yates, R.A., Snider, T.: New indices for text: Pat trees and pat arrays. Inf. Retrieval: Data Struc. Algorithms 66, 82 (1992)
Ziv, J., Lempel, A.: A universal algorithm for sequential data compression. IEEE Trans. Inf. Theory 23(3), 337–343 (1977)
Larsson, N.J., Moffat, A.: Off-line dictionary-based compression. Proc. IEEE 88(11), 1722–1732 (2000). https://doi.org/10.1109/5.892708
Ferragina, P., Manzini, G.: Opportunistic data structures with applications. In: Proceedings 41st Annual Symposium on Foundations of Computer Science, pp. 390–398 (2000). https://doi.org/10.1109/SFCS.2000.892127
Grossi, R., Vitter, J.S.: Compressed suffix arrays and suffix trees with applications to text indexing and string matching. In: Proceedings of the Thirty-second Annual ACM Symposium on Theory of Computing (STOC), pp. 397–406 (2000). https://doi.org/10.1145/335305.33535
Mäkinen, V., Navarro, G.: Succinct suffix arrays based on run-length encoding. Nordic J. Comput. 12(1), 40–66 (2005)
Gagie, T., Navarro, G., Prezza, N.: Fully functional suffix trees and optimal text searching in BWT-runs bounded space. J. ACM. 67(1) (2020). https://doi.org/10.1145/3375890
Puglisi, S.J., Zhukova, B.: Smaller RLZ-compressed suffix arrays. In: 2021 Data Compression Conference (DCC), pp. 213–222. IEEE (2021). https://doi.org/10.1109/DCC50243.2021.00029
Depuydt, L., Renders, L., Vyver, S., Veys, L., Gagie, T., Fostier, J.: b-move: Faster Bidirectional Character Extensions in a Run-Length Compressed Index. In: Pissis, S.P., Sung, W.-K. (eds.) 24th International Workshop on Algorithms in Bioinformatics (WABI 2024). Leibniz International Proceedings in Informatics (LIPIcs), vol. 312, pp. 10–11018. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany (2024). https://doi.org/10.4230/LIPIcs.WABI.2024.10
Li, H.: Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arxiv:1303.3997 (2013)
Bannai, H., Gagie, T., I., T.: Refining r-index. Theoret. Comput. Sci. 812, 92–108 (2020). https://doi.org/10.1016/j.tcs.2019.08.005
Rossi, M., Oliva, M., Langmead, B., Gagie, T., Boucher, C.: MONI: A pangenomic index for finding maximal exact matches. J. Comput. Biol. 29(2), 169–187 (2022). https://doi.org/10.1089/CMB.2021.0290
Boucher, C., Gagie, T., I, T., Köppl, D., Langmead, B., Manzini, G., Navarro, G., Pacheco, A., Rossi, M.: PHONI: Streamed Matching Statistics with Multi-Genome References. In: Data Compression Conference (DCC), pp. 193–202 (2021). https://doi.org/10.1109/DCC50243.2021.00027
Nishimoto, T., Tabei, Y.: Optimal-time queries on BWT-runs compressed indexes. In: 48th International Colloquium on Automata, Languages, and Programming (ICALP 2021), pp. 101–1. Schloss Dagstuhl–Leibniz-Zentrum für Informatik, (2021). https://doi.org/10.4230/LIPIcs.ICALP.2021.101
Bille, P., Landau, G.M., Raman, R., Sadakane, K., Satti, S.R., Weimann, O.: Random access to grammar-compressed strings and trees. SIAM J. Comput. 44, 513–539 (2015). https://doi.org/10.1137/130936889
Kreft, S., Navarro, G.: On compressing and indexing repetitive sequences. Theoret. Comput. Sci. 483, 115–133 (2013). https://doi.org/10.1016/j.tcs.2012.02.006
Belazzougui, D., Cáceres, M., Gagie, T., Gawrychowski, P., Kärkkäinen, J., Navarro, G., Ordóñez, A., Puglisi, S.J., Tabei, Y.: Block trees. J. Comput. Syst. Sci. 117, 1–22 (2021). https://doi.org/10.1016/j.jcss.2020.11.002
Kuruppu, S., Puglisi, S.J., Zobel, J.: Relative Lempel-Ziv compression of genomes for large-scale storage and retrieval. In: International Symposium on String Processing and Information Retrieval, pp. 201–206. Springer, (2010). https://doi.org/10.1007/978-3-642-16321-0_20
Bentley, J.W., Gibney, D., Thankachan, S.V.: On the complexity of BWT-runs minimization via alphabet reordering. In: Proceedings of the 28th Annual European Symposium on Algorithms (ESA 2020). LIPIcs, vol. 173, pp. 15–11513 (2020). https://doi.org/10.4230/LIPICS.ESA.2020.15
Kempa, D., Prezza, N.: At the roots of dictionary compression: string attractors. In: Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing (STOC 2018). STOC 2018, pp. 827–840. Association for Computing Machinery, New York, NY, USA (2018). https://doi.org/10.1145/3188745.3188814 . https://doi.org/10.1145/3188745.3188814
Belazzougui, D., Boldi, P., Pagh, R., Vigna, S.: Monotone minimal perfect hashing: searching a sorted table with \(o(1)\) accesses. In: Proceedings of the Twentieth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 785–794 (2009). https://dl.acm.org/doi/10.5555/1496770.1496856
Boucher, C., Gagie, T., Kuhnle, A., Langmead, B., Manzini, G., Mun, T.: Prefix-free parsing for building big BWTs. Algorithms Mol. Biol. 14(1), 13–11315 (2019). https://doi.org/10.1186/S13015-019-0148-5
Li, H.: Exploring single-sample snp and indel calling with whole-genome de novo assembly. Bioinformatics 28(14), 1838–1844 (2012). https://doi.org/10.1093/bioinformatics/bts280
Manber, U., Myers, G.: Suffix arrays: A new method for on-line string searches. SIAM J. Comput. 22(5), 935–948 (1993). https://doi.org/10.1137/0222058
Weiner, P.: Linear pattern matching algorithms. In: Proceedings of the 14th Annual Symposium on Switching and Automata Theory, pp. 1–11 (1973). https://doi.org/10.1109/SWAT.1973.13
Burrows, M., Wheeler, D.J.: A block sorting lossless data compression algorithm, Digital Equipment Corporation (1994). (124 Technical Report)
Bentley, J.W., Gibney, D., Thankachan, S.V.: On the Complexity of BWT-Runs Minimization via Alphabet Reordering. In: 28th Annual European Symposium on Algorithms (ESA 2020). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, (2020). https://doi.org/10.4230/LIPIcs.ESA.2020.15
Sirén, J., Välimäki, N., Mäkinen, V., Navarro, G.: Run-length compressed indexes are superior for highly repetitive sequence collections. In: Amir, A., Turpin, A., Moffat, A. (eds.) String Processing and Information Retrieval, 15th International Symposium, SPIRE 2008, Melbourne, Australia, November 10-12, 2008. Proceedings. Lecture Notes in Computer Science, vol. 5280, pp. 164–175. Springer, (2008). https://doi.org/10.1007/978-3-540-89097-3_17
Boldi, P., Vigna, S.: Kings, Name Days, Lazy Servants and Magic. In: Ito, H., Leonardi, S., Pagli, L., Prencipe, G. (eds.) 9th International Conference on Fun with Algorithms (FUN 2018). Leibniz International Proceedings in Informatics (LIPIcs), vol. 100, pp. 10–11013. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany (2018). https://doi.org/10.4230/LIPIcs.FUN.2018.10
Baláž, A., Gagie, T., Goga, A., Heumos, S., Navarro, G., Petescia, A., Sirén, J.: Wheeler maps. In: Soto, J.A., Wiese, A. (eds.) LATIN 2024: Theoretical Informatics, pp. 178–192. Springer, Cham (2024). https://doi.org/10.1007/978-3-031-55598-5_12
Bille, P., Gørtz, I.L., Cording, P.H., Sach, B., Vildhøj, H.W., Vind, S.: Fingerprints in compressed strings. J. Comput. Syst. Sci. 86, 171–180 (2017). https://doi.org/10.1016/j.jcss.2017.01.002
Bloom, B.H.: Space/time trade-offs in hash coding with allowable errors. Commun. ACM 13(7), 422–426 (1970). https://doi.org/10.1145/362686.362692
Li, H.: Exploring single-sample SNP and INDEL calling with whole-genome de novo assembly. Bioinformatics 28(14), 1838–1844 (2012). https://doi.org/10.1093/bioinformatics/bts280
Gagie, T.: How to find long maximal exact matches and ignore short ones. In: International Conference on Developments in Language Theory, pp. 131–140 (2024). https://doi.org/10.1007/978-3-031-66159-4_10
Kuhnle, A., Mun, T., Boucher, C., Gagie, T., Langmead, B., Manzini, G.: Efficient construction of a complete index for pan-genomics read alignment. J. Comput. Biol. 27(4), 500–513 (2020). https://doi.org/10.1089/CMB.2019.0309
Fischer, J., Mäkinen, V., Navarro, G.: Faster entropy-bounded compressed suffix trees. Theoret. Comput. Sci. 410(51), 5354–5364 (2009)
Berkman, O., Schieber, B., Vishkin, U.: Optimal doubly logarithmic parallel algorithms based on finding all nearest smaller values. J. Algorithms 14(3), 344–370 (1993). https://doi.org/10.1006/jagm.1993.1018
Belazzougui, D., Navarro, G.: Optimal lower and upper bounds for representing sequences. ACM Trans. Algorithms 11(4), 31–13121 (2015). https://doi.org/10.1145/2629339
Navarro, G., Romana, G., Urbina, C.: Smallest suffixient sets as a repetitiveness measure. In: Proceedings of the 32nd International Symposium on String Processing and Information Retrieval, pp. 217–232 (2025). https://doi.org/10.1007/978-3-032-05228-5_18
Kempa, D., Kociumaka, T.: Resolution of the Burrows-Wheeler transform conjecture. Commun. ACM 65(6), 91–98 (2022). https://doi.org/10.1145/3531445
Lam, T.W., Li, R., Tam, A., Wong, S., Wu, E., Yiu, S.: High throughput short read alignment via bi-directional BWT. In: IEEE International Conference on Bioinformatics and Biomedicine, pp. 31–36 (2009). https://doi.org/10.1109/BIBM.2009.42
Kucherov, G., Salikhov, K., Tsur, D.: Approximate String Matching Using a Bidirectional Index. In: Kulikov, A.S., Kuznetsov, S.O., Pevzner, P. (eds.) Combinatorial Pattern Matching, pp. 222–231. Springer, Cham (2014). https://doi.org/10.1007/978-3-319-07566-2_23
Acknowledgements
We thank Ragnar Groot Koerkamp and Giulio Ermanno Pibiri for fruitful discussions on the topic.
Funding
Davide Cenzato, Sung-Hwan Kim, Nicola Prezza: Funded by the European Union (ERC, REGINDEX, 101039208). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.
Lore Depuydt: Funded by a PhD Fellowship FR (1117322N), Research Foundation – Flanders (FWO).
Travis Gagie: Funded by NSERC Discovery Grant RGPIN-07185-2020.
Giovanni Manzini: Funded by the NextGenerationEU programme PNRR ECS00000017 Tuscany Health Ecosystem (CUP: I53C22000780001) and by the project PAN-HUB, Italian Ministry of Health (CUP: I53C22001300001).
Francisco Olivares: Funded by scholarship ANID-Subdirección de Capital Humano/Doctorado Nacional/2021-21210579, ANID, Chile and by Basal Funds FB0001, ANID, Chile.
Author information
Authors and Affiliations
Contributions
D.C., L.D., T.G., S.K., G.M., F.O., N.P. contributed equally to the paper.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/
About this article
Cite this article
Cenzato, D., Depuydt, L., Gagie, T. et al. Suffixient Arrays: A New Efficient Suffix Array Compression Technique. Theory Comput Syst 70, 50 (2026). https://doi.org/10.1007/s00224-026-10287-6
Received:
Accepted:
Published:
Version of record:
DOI: https://doi.org/10.1007/s00224-026-10287-6