a2d6a1e0cd95200d1f3a5464d8359d5736b14cbe
[wimlib] / src / lcpit_matchfinder.c
1 /*
2  * lcpit_matchfinder.c
3  *
4  * A match-finder for Lempel-Ziv compression based on bottom-up construction and
5  * traversal of the Longest Common Prefix (LCP) interval tree.
6  *
7  * The following copying information applies to this specific source code file:
8  *
9  * Written in 2014-2015 by Eric Biggers <ebiggers3@gmail.com>
10  *
11  * To the extent possible under law, the author(s) have dedicated all copyright
12  * and related and neighboring rights to this software to the public domain
13  * worldwide via the Creative Commons Zero 1.0 Universal Public Domain
14  * Dedication (the "CC0").
15  *
16  * This software is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18  * FOR A PARTICULAR PURPOSE. See the CC0 for more details.
19  *
20  * You should have received a copy of the CC0 along with this software; if not
21  * see <http://creativecommons.org/publicdomain/zero/1.0/>.
22  */
23
24 #ifdef HAVE_CONFIG_H
25 #  include "config.h"
26 #endif
27
28 #include <limits.h>
29
30 #include "wimlib/divsufsort.h"
31 #include "wimlib/lcpit_matchfinder.h"
32 #include "wimlib/util.h"
33
34 #define LCP_BITS                6
35 #define LCP_MAX                 (((u32)1 << LCP_BITS) - 1)
36 #define LCP_SHIFT               (32 - LCP_BITS)
37 #define LCP_MASK                (LCP_MAX << LCP_SHIFT)
38 #define POS_MASK                (((u32)1 << (32 - LCP_BITS)) - 1)
39 #define MAX_NORMAL_BUFSIZE      (POS_MASK + 1)
40
41 #define HUGE_LCP_BITS           7
42 #define HUGE_LCP_MAX            (((u32)1 << HUGE_LCP_BITS) - 1)
43 #define HUGE_LCP_SHIFT          (64 - HUGE_LCP_BITS)
44 #define HUGE_LCP_MASK           ((u64)HUGE_LCP_MAX << HUGE_LCP_SHIFT)
45 #define HUGE_POS_MASK           0xFFFFFFFF
46 #define MAX_HUGE_BUFSIZE        ((u64)HUGE_POS_MASK + 1)
47 #define HUGE_UNVISITED_TAG      0x100000000
48
49 #define PREFETCH_SAFETY         5
50
51 /*
52  * Build the LCP (Longest Common Prefix) array in linear time.
53  *
54  * LCP[r] will be the length of the longest common prefix between the suffixes
55  * with positions SA[r - 1] and SA[r].  LCP[0] will be undefined.
56  *
57  * Algorithm taken from Kasai et al. (2001), but modified slightly:
58  *
59  *  - With bytes there is no realistic way to reserve a unique symbol for
60  *    end-of-buffer, so use explicit checks for end-of-buffer.
61  *
62  *  - For decreased memory usage and improved memory locality, pack the two
63  *    logically distinct SA and LCP arrays into a single array SA_and_LCP.
64  *
65  *  - Since SA_and_LCP is accessed randomly, improve the cache behavior by
66  *    reading several entries ahead in ISA and prefetching the upcoming
67  *    SA_and_LCP entry.
68  *
69  *  - If an LCP value is less than the minimum match length, then store 0.  This
70  *    avoids having to do comparisons against the minimum match length later.
71  *
72  *  - If an LCP value is greater than the "nice match length", then store the
73  *    "nice match length".  This caps the number of bits needed to store each
74  *    LCP value, and this caps the depth of the LCP-interval tree, without
75  *    usually hurting the compression ratio too much.
76  *
77  * References:
78  *
79  *      Kasai et al.  2001.  Linear-Time Longest-Common-Prefix Computation in
80  *      Suffix Arrays and Its Applications.  CPM '01 Proceedings of the 12th
81  *      Annual Symposium on Combinatorial Pattern Matching pp. 181-192.
82  */
83 static void
84 build_LCP(u32 SA_and_LCP[restrict], const u32 ISA[restrict],
85           const u8 T[restrict], const u32 n,
86           const u32 min_lcp, const u32 max_lcp)
87 {
88         u32 h = 0;
89         for (u32 i = 0; i < n; i++) {
90                 const u32 r = ISA[i];
91                 prefetchw(&SA_and_LCP[ISA[i + PREFETCH_SAFETY]]);
92                 if (r > 0) {
93                         const u32 j = SA_and_LCP[r - 1] & POS_MASK;
94                         const u32 lim = min(n - i, n - j);
95                         while (h < lim && T[i + h] == T[j + h])
96                                 h++;
97                         u32 stored_lcp = h;
98                         if (stored_lcp < min_lcp)
99                                 stored_lcp = 0;
100                         else if (stored_lcp > max_lcp)
101                                 stored_lcp = max_lcp;
102                         SA_and_LCP[r] |= stored_lcp << LCP_SHIFT;
103                         if (h > 0)
104                                 h--;
105                 }
106         }
107 }
108
109 /*
110  * Use the suffix array accompanied with the longest-common-prefix array ---
111  * which in combination can be called the "enhanced suffix array" --- to
112  * simulate a bottom-up traversal of the corresponding suffix tree, or
113  * equivalently the lcp-interval tree.  Do so in suffix rank order, but save the
114  * superinterval references needed for later bottom-up traversal of the tree in
115  * suffix position order.
116  *
117  * To enumerate the lcp-intervals, this algorithm scans the suffix array and its
118  * corresponding LCP array linearly.  While doing so, it maintains a stack of
119  * lcp-intervals that are currently open, meaning that their left boundaries
120  * have been seen but their right boundaries have not.  The bottom of the stack
121  * is the interval which covers the entire suffix array (this has lcp=0), and
122  * the top of the stack is the deepest interval that is currently open (this has
123  * the greatest lcp of any interval on the stack).  When this algorithm opens an
124  * lcp-interval, it assigns it a unique index in intervals[] and pushes it onto
125  * the stack.  When this algorithm closes an interval, it pops it from the stack
126  * and sets the intervals[] entry of that interval to the index and lcp of that
127  * interval's superinterval, which is the new top of the stack.
128  *
129  * This algorithm also set pos_data[pos] for each suffix position 'pos' to the
130  * index and lcp of the deepest lcp-interval containing it.  Alternatively, we
131  * can interpret each suffix as being associated with a singleton lcp-interval,
132  * or leaf of the suffix tree.  With this interpretation, an entry in pos_data[]
133  * is the superinterval reference for one of these singleton lcp-intervals and
134  * therefore is not fundamentally different from an entry in intervals[].
135  *
136  * To reduce memory usage, this algorithm re-uses the suffix array's memory to
137  * store the generated intervals[] array.  This is possible because SA and LCP
138  * are accessed linearly, and no more than one interval is generated per suffix.
139  *
140  * The techniques used in this algorithm are described in various published
141  * papers.  The generation of lcp-intervals from the suffix array (SA) and the
142  * longest-common-prefix array (LCP) is given as Algorithm BottomUpTraverse in
143  * Kasai et al. (2001) and Algorithm 4.1 ("Computation of lcp-intervals") in
144  * Abouelhoda et al. (2004).  Both these papers note the equivalence between
145  * lcp-intervals (including the singleton lcp-interval for each suffix) and
146  * nodes of the suffix tree.  Abouelhoda et al. (2004) furthermore applies
147  * bottom-up traversal of the lcp-interval tree to Lempel-Ziv factorization, as
148  * does Chen at al. (2008).  Algorithm CPS1b of Chen et al. (2008) dynamically
149  * re-uses the suffix array during bottom-up traversal of the lcp-interval tree.
150  *
151  * References:
152  *
153  *      Kasai et al. Linear-Time Longest-Common-Prefix Computation in Suffix
154  *      Arrays and Its Applications.  2001.  CPM '01 Proceedings of the 12th
155  *      Annual Symposium on Combinatorial Pattern Matching pp. 181-192.
156  *
157  *      M.I. Abouelhoda, S. Kurtz, E. Ohlebusch.  2004.  Replacing Suffix Trees
158  *      With Enhanced Suffix Arrays.  Journal of Discrete Algorithms Volume 2
159  *      Issue 1, March 2004, pp. 53-86.
160  *
161  *      G. Chen, S.J. Puglisi, W.F. Smyth.  2008.  Lempel-Ziv Factorization
162  *      Using Less Time & Space.  Mathematics in Computer Science June 2008,
163  *      Volume 1, Issue 4, pp. 605-623.
164  */
165 static void
166 build_LCPIT(u32 intervals[restrict], u32 pos_data[restrict], const u32 n)
167 {
168         u32 * const SA_and_LCP = intervals;
169         u32 next_interval_idx;
170         u32 open_intervals[LCP_MAX + 1];
171         u32 *top = open_intervals;
172         u32 prev_pos = SA_and_LCP[0] & POS_MASK;
173
174         *top = 0;
175         intervals[0] = 0;
176         next_interval_idx = 1;
177
178         for (u32 r = 1; r < n; r++) {
179                 const u32 next_pos = SA_and_LCP[r] & POS_MASK;
180                 const u32 next_lcp = SA_and_LCP[r] & LCP_MASK;
181                 const u32 top_lcp = *top & LCP_MASK;
182
183                 prefetchw(&pos_data[SA_and_LCP[r + PREFETCH_SAFETY] & POS_MASK]);
184
185                 if (next_lcp == top_lcp) {
186                         /* Continuing the deepest open interval  */
187                         pos_data[prev_pos] = *top;
188                 } else if (next_lcp > top_lcp) {
189                         /* Opening a new interval  */
190                         *++top = next_lcp | next_interval_idx++;
191                         pos_data[prev_pos] = *top;
192                 } else {
193                         /* Closing the deepest open interval  */
194                         pos_data[prev_pos] = *top;
195                         for (;;) {
196                                 const u32 closed_interval_idx = *top-- & POS_MASK;
197                                 const u32 superinterval_lcp = *top & LCP_MASK;
198
199                                 if (next_lcp == superinterval_lcp) {
200                                         /* Continuing the superinterval */
201                                         intervals[closed_interval_idx] = *top;
202                                         break;
203                                 } else if (next_lcp > superinterval_lcp) {
204                                         /* Creating a new interval that is a
205                                          * superinterval of the one being
206                                          * closed, but still a subinterval of
207                                          * its superinterval  */
208                                         *++top = next_lcp | next_interval_idx++;
209                                         intervals[closed_interval_idx] = *top;
210                                         break;
211                                 } else {
212                                         /* Also closing the superinterval  */
213                                         intervals[closed_interval_idx] = *top;
214                                 }
215                         }
216                 }
217                 prev_pos = next_pos;
218         }
219
220         /* Close any still-open intervals.  */
221         pos_data[prev_pos] = *top;
222         for (; top > open_intervals; top--)
223                 intervals[*top & POS_MASK] = *(top - 1);
224 }
225
226 /*
227  * Advance the LCP-interval tree matchfinder by one byte.
228  *
229  * If @record_matches is true, then matches are written to the @matches array
230  * sorted by strictly decreasing length and strictly decreasing offset, and the
231  * return value is the number of matches found.  Otherwise, @matches is ignored
232  * and the return value is always 0.
233  *
234  * How this algorithm works:
235  *
236  * 'cur_pos' is the position of the current suffix, which is the suffix being
237  * matched against.  'cur_pos' starts at 0 and is incremented each time this
238  * function is called.  This function finds each suffix with position less than
239  * 'cur_pos' that shares a prefix with the current suffix, but for each distinct
240  * prefix length it finds only the suffix with the greatest position (i.e. the
241  * most recently seen in the linear traversal by position).  This function
242  * accomplishes this using the lcp-interval tree data structure that was built
243  * by build_LCPIT() and is updated dynamically by this function.
244  *
245  * The first step is to read 'pos_data[cur_pos]', which gives the index and lcp
246  * value of the deepest lcp-interval containing the current suffix --- or,
247  * equivalently, the parent of the conceptual singleton lcp-interval that
248  * contains the current suffix.
249  *
250  * The second step is to ascend the lcp-interval tree until reaching an interval
251  * that has not yet been visited, and link the intervals to the current suffix
252  * along the way.   An lcp-interval has been visited if and only if it has been
253  * linked to a suffix.  Initially, no lcp-intervals are linked to suffixes.
254  *
255  * The third step is to continue ascending the lcp-interval tree, but indirectly
256  * through suffix links rather than through the original superinterval
257  * references, and continuing to form links with the current suffix along the
258  * way.  Each suffix visited during this step, except in a special case to
259  * handle outdated suffixes, is a match which can be written to matches[].  Each
260  * intervals[] entry contains the position of the next suffix to visit, which we
261  * shall call 'match_pos'; this is the most recently seen suffix that belongs to
262  * that lcp-interval.  'pos_data[match_pos]' then contains the lcp and interval
263  * index of the next lcp-interval that should be visited.
264  *
265  * We can view these arrays as describing a new set of links that gets overlaid
266  * on top of the original superinterval references of the lcp-interval tree.
267  * Each such link can connect a node of the lcp-interval tree to an ancestor
268  * more than one generation removed.
269  *
270  * For each one-byte advance, the current position becomes the most recently
271  * seen suffix for a continuous sequence of lcp-intervals from a leaf interval
272  * to the root interval.  Conceptually, this algorithm needs to update all these
273  * nodes to link to 'cur_pos', and then update 'pos_data[cur_pos]' to a "null"
274  * link.  But actually, if some of these nodes have the same most recently seen
275  * suffix, then this algorithm just visits the pos_data[] entry for that suffix
276  * and skips over all these nodes in one step.  Updating the extra nodes is
277  * accomplished by creating a redirection from the previous suffix to the
278  * current suffix.
279  *
280  * Using this shortcutting scheme, it's possible for a suffix to become out of
281  * date, which means that it is no longer the most recently seen suffix for the
282  * lcp-interval under consideration.  This case is detected by noticing when the
283  * next lcp-interval link actually points deeper in the tree, and it is worked
284  * around by just continuing until we get to a link that actually takes us
285  * higher in the tree.  This can be described as a lazy-update scheme.
286  */
287 static forceinline u32
288 lcpit_advance_one_byte(const u32 cur_pos,
289                        u32 pos_data[restrict],
290                        u32 intervals[restrict],
291                        u32 next[restrict],
292                        struct lz_match matches[restrict],
293                        const bool record_matches)
294 {
295         u32 ref;
296         u32 super_ref;
297         u32 match_pos;
298         struct lz_match *matchptr;
299
300         /* Get the deepest lcp-interval containing the current suffix. */
301         ref = pos_data[cur_pos];
302
303         /* Prefetch upcoming data, up to 3 positions ahead.  Assume the
304          * intervals are already visited.  */
305
306         /* Prefetch the superinterval via a suffix link for the deepest
307          * lcp-interval containing the suffix starting 1 position from now.  */
308         prefetchw(&intervals[pos_data[next[0]] & POS_MASK]);
309
310         /* Prefetch suffix link for the deepest lcp-interval containing the
311          * suffix starting 2 positions from now.  */
312         next[0] = intervals[next[1]] & POS_MASK;
313         prefetchw(&pos_data[next[0]]);
314
315         /* Prefetch the deepest lcp-interval containing the suffix starting 3
316          * positions from now.  */
317         next[1] = pos_data[cur_pos + 3] & POS_MASK;
318         prefetchw(&intervals[next[1]]);
319
320         /* There is no "next suffix" after the current one.  */
321         pos_data[cur_pos] = 0;
322
323         /* Ascend until we reach a visited interval, the root, or a child of the
324          * root.  Link unvisited intervals to the current suffix as we go.  */
325         while ((super_ref = intervals[ref & POS_MASK]) & LCP_MASK) {
326                 intervals[ref & POS_MASK] = cur_pos;
327                 ref = super_ref;
328         }
329
330         if (super_ref == 0) {
331                 /* In this case, the current interval may be any of:
332                  * (1) the root;
333                  * (2) an unvisited child of the root;
334                  * (3) an interval last visited by suffix 0
335                  *
336                  * We could avoid the ambiguity with (3) by using an lcp
337                  * placeholder value other than 0 to represent "visited", but
338                  * it's fastest to use 0.  So we just don't allow matches with
339                  * position 0.  */
340
341                 if (ref != 0)  /* Not the root?  */
342                         intervals[ref & POS_MASK] = cur_pos;
343                 return 0;
344         }
345
346         /* Ascend indirectly via pos_data[] links.  */
347         match_pos = super_ref;
348         matchptr = matches;
349         for (;;) {
350                 while ((super_ref = pos_data[match_pos]) > ref)
351                         match_pos = intervals[super_ref & POS_MASK];
352                 intervals[ref & POS_MASK] = cur_pos;
353                 pos_data[match_pos] = ref;
354                 if (record_matches) {
355                         matchptr->length = ref >> LCP_SHIFT;
356                         matchptr->offset = cur_pos - match_pos;
357                         matchptr++;
358                 }
359                 if (super_ref == 0)
360                         break;
361                 ref = super_ref;
362                 match_pos = intervals[ref & POS_MASK];
363         }
364         return matchptr - matches;
365 }
366
367 /* Expand SA from 32 bits to 64 bits.  */
368 static void
369 expand_SA(void *p, u32 n)
370 {
371         typedef u32 _may_alias_attribute aliased_u32_t;
372         typedef u64 _may_alias_attribute aliased_u64_t;
373
374         aliased_u32_t *SA = p;
375         aliased_u64_t *SA64 = p;
376
377         u32 r = n - 1;
378         do {
379                 SA64[r] = SA[r];
380         } while (r--);
381 }
382
383 /* Like build_LCP(), but for buffers larger than MAX_NORMAL_BUFSIZE.  */
384 static void
385 build_LCP_huge(u64 SA_and_LCP64[restrict], const u32 ISA[restrict],
386                const u8 T[restrict], const u32 n,
387                const u32 min_lcp, const u32 max_lcp)
388 {
389         u32 h = 0;
390         for (u32 i = 0; i < n; i++) {
391                 const u32 r = ISA[i];
392                 prefetchw(&SA_and_LCP64[ISA[i + PREFETCH_SAFETY]]);
393                 if (r > 0) {
394                         const u32 j = SA_and_LCP64[r - 1] & HUGE_POS_MASK;
395                         const u32 lim = min(n - i, n - j);
396                         while (h < lim && T[i + h] == T[j + h])
397                                 h++;
398                         u32 stored_lcp = h;
399                         if (stored_lcp < min_lcp)
400                                 stored_lcp = 0;
401                         else if (stored_lcp > max_lcp)
402                                 stored_lcp = max_lcp;
403                         SA_and_LCP64[r] |= (u64)stored_lcp << HUGE_LCP_SHIFT;
404                         if (h > 0)
405                                 h--;
406                 }
407         }
408 }
409
410 /*
411  * Like build_LCPIT(), but for buffers larger than MAX_NORMAL_BUFSIZE.
412  *
413  * This "huge" version is also slightly different in that the lcp value stored
414  * in each intervals[] entry is the lcp value for that interval, not its
415  * superinterval.  This lcp value stays put in intervals[] and doesn't get moved
416  * to pos_data[] during lcpit_advance_one_byte_huge().  One consequence of this
417  * is that we have to use a special flag to distinguish visited from unvisited
418  * intervals.  But overall, this scheme keeps the memory usage at 12n instead of
419  * 16n.  (The non-huge version is 8n.)
420  */
421 static void
422 build_LCPIT_huge(u64 intervals64[restrict], u32 pos_data[restrict], const u32 n)
423 {
424         u64 * const SA_and_LCP64 = intervals64;
425         u32 next_interval_idx;
426         u32 open_intervals[HUGE_LCP_MAX + 1];
427         u32 *top = open_intervals;
428         u32 prev_pos = SA_and_LCP64[0] & HUGE_POS_MASK;
429
430         *top = 0;
431         intervals64[0] = 0;
432         next_interval_idx = 1;
433
434         for (u32 r = 1; r < n; r++) {
435                 const u32 next_pos = SA_and_LCP64[r] & HUGE_POS_MASK;
436                 const u64 next_lcp = SA_and_LCP64[r] & HUGE_LCP_MASK;
437                 const u64 top_lcp = intervals64[*top];
438
439                 prefetchw(&pos_data[SA_and_LCP64[r + PREFETCH_SAFETY] & HUGE_POS_MASK]);
440
441                 if (next_lcp == top_lcp) {
442                         /* Continuing the deepest open interval  */
443                         pos_data[prev_pos] = *top;
444                 } else if (next_lcp > top_lcp) {
445                         /* Opening a new interval  */
446                         intervals64[next_interval_idx] = next_lcp;
447                         pos_data[prev_pos] = next_interval_idx;
448                         *++top = next_interval_idx++;
449                 } else {
450                         /* Closing the deepest open interval  */
451                         pos_data[prev_pos] = *top;
452                         for (;;) {
453                                 const u32 closed_interval_idx = *top--;
454                                 const u64 superinterval_lcp = intervals64[*top];
455
456                                 if (next_lcp == superinterval_lcp) {
457                                         /* Continuing the superinterval */
458                                         intervals64[closed_interval_idx] |=
459                                                 HUGE_UNVISITED_TAG | *top;
460                                         break;
461                                 } else if (next_lcp > superinterval_lcp) {
462                                         /* Creating a new interval that is a
463                                          * superinterval of the one being
464                                          * closed, but still a subinterval of
465                                          * its superinterval  */
466                                         intervals64[next_interval_idx] = next_lcp;
467                                         intervals64[closed_interval_idx] |=
468                                                 HUGE_UNVISITED_TAG | next_interval_idx;
469                                         *++top = next_interval_idx++;
470                                         break;
471                                 } else {
472                                         /* Also closing the superinterval  */
473                                         intervals64[closed_interval_idx] |=
474                                                 HUGE_UNVISITED_TAG | *top;
475                                 }
476                         }
477                 }
478                 prev_pos = next_pos;
479         }
480
481         /* Close any still-open intervals.  */
482         pos_data[prev_pos] = *top;
483         for (; top > open_intervals; top--)
484                 intervals64[*top] |= HUGE_UNVISITED_TAG | *(top - 1);
485 }
486
487 /* Like lcpit_advance_one_byte(), but for buffers larger than
488  * MAX_NORMAL_BUFSIZE.  */
489 static forceinline u32
490 lcpit_advance_one_byte_huge(const u32 cur_pos,
491                             u32 pos_data[restrict],
492                             u64 intervals64[restrict],
493                             u32 prefetch_next[restrict],
494                             struct lz_match matches[restrict],
495                             const bool record_matches)
496 {
497         u32 interval_idx;
498         u32 next_interval_idx;
499         u64 cur;
500         u64 next;
501         u32 match_pos;
502         struct lz_match *matchptr;
503
504         interval_idx = pos_data[cur_pos];
505
506         prefetchw(&intervals64[pos_data[prefetch_next[0]] & HUGE_POS_MASK]);
507
508         prefetch_next[0] = intervals64[prefetch_next[1]] & HUGE_POS_MASK;
509         prefetchw(&pos_data[prefetch_next[0]]);
510
511         prefetch_next[1] = pos_data[cur_pos + 3] & HUGE_POS_MASK;
512         prefetchw(&intervals64[prefetch_next[1]]);
513
514         pos_data[cur_pos] = 0;
515
516         while ((next = intervals64[interval_idx]) & HUGE_UNVISITED_TAG) {
517                 intervals64[interval_idx] = (next & HUGE_LCP_MASK) | cur_pos;
518                 interval_idx = next & HUGE_POS_MASK;
519         }
520
521         matchptr = matches;
522         while (next & HUGE_LCP_MASK) {
523                 cur = next;
524                 do {
525                         match_pos = next & HUGE_POS_MASK;
526                         next_interval_idx = pos_data[match_pos];
527                         next = intervals64[next_interval_idx];
528                 } while (next > cur);
529                 intervals64[interval_idx] = (cur & HUGE_LCP_MASK) | cur_pos;
530                 pos_data[match_pos] = interval_idx;
531                 if (record_matches) {
532                         matchptr->length = cur >> HUGE_LCP_SHIFT;
533                         matchptr->offset = cur_pos - match_pos;
534                         matchptr++;
535                 }
536                 interval_idx = next_interval_idx;
537         }
538         return matchptr - matches;
539 }
540
541 static forceinline u64
542 get_pos_data_size(size_t max_bufsize)
543 {
544         return (u64)max((u64)max_bufsize + PREFETCH_SAFETY,
545                         DIVSUFSORT_TMP_LEN) * sizeof(u32);
546 }
547
548 static forceinline u64
549 get_intervals_size(size_t max_bufsize)
550 {
551         return ((u64)max_bufsize + PREFETCH_SAFETY) *
552                 (max_bufsize <= MAX_NORMAL_BUFSIZE ? sizeof(u32) : sizeof(u64));
553 }
554
555 /*
556  * Calculate the number of bytes of memory needed for the LCP-interval tree
557  * matchfinder.
558  *
559  * @max_bufsize - maximum buffer size to support
560  *
561  * Returns the number of bytes required.
562  */
563 u64
564 lcpit_matchfinder_get_needed_memory(size_t max_bufsize)
565 {
566         return get_pos_data_size(max_bufsize) + get_intervals_size(max_bufsize);
567 }
568
569 /*
570  * Initialize the LCP-interval tree matchfinder.
571  *
572  * @mf - the matchfinder structure to initialize
573  * @max_bufsize - maximum buffer size to support
574  * @min_match_len - minimum match length in bytes
575  * @nice_match_len - only consider this many bytes of each match
576  *
577  * Returns true if successfully initialized; false if out of memory.
578  */
579 bool
580 lcpit_matchfinder_init(struct lcpit_matchfinder *mf, size_t max_bufsize,
581                        u32 min_match_len, u32 nice_match_len)
582 {
583         if (lcpit_matchfinder_get_needed_memory(max_bufsize) > SIZE_MAX)
584                 return false;
585         if (max_bufsize > MAX_HUGE_BUFSIZE - PREFETCH_SAFETY)
586                 return false;
587
588         mf->pos_data = MALLOC(get_pos_data_size(max_bufsize));
589         mf->intervals = MALLOC(get_intervals_size(max_bufsize));
590         if (!mf->pos_data || !mf->intervals) {
591                 lcpit_matchfinder_destroy(mf);
592                 return false;
593         }
594
595         mf->min_match_len = min_match_len;
596         mf->nice_match_len = min(nice_match_len,
597                                  (max_bufsize <= MAX_NORMAL_BUFSIZE) ?
598                                  LCP_MAX : HUGE_LCP_MAX);
599         return true;
600 }
601
602 /*
603  * Build the suffix array SA for the specified byte array T of length n.
604  *
605  * The suffix array is a sorted array of the byte array's suffixes, represented
606  * by indices into the byte array.  It can equivalently be viewed as a mapping
607  * from suffix rank to suffix position.
608  *
609  * To build the suffix array, we use libdivsufsort, which uses an
610  * induced-sorting-based algorithm.  In practice, this seems to be the fastest
611  * suffix array construction algorithm currently available.
612  *
613  * References:
614  *
615  *      Y. Mori.  libdivsufsort, a lightweight suffix-sorting library.
616  *      https://github.com/y-256/libdivsufsort
617  *
618  *      G. Nong, S. Zhang, and W.H. Chan.  2009.  Linear Suffix Array
619  *      Construction by Almost Pure Induced-Sorting.  Data Compression
620  *      Conference, 2009.  DCC '09.  pp. 193 - 202.
621  *
622  *      S.J. Puglisi, W.F. Smyth, and A. Turpin.  2007.  A Taxonomy of Suffix
623  *      Array Construction Algorithms.  ACM Computing Surveys (CSUR) Volume 39
624  *      Issue 2, 2007 Article No. 4.
625  */
626 static void
627 build_SA(u32 SA[], const u8 T[], u32 n, u32 *tmp)
628 {
629         /* Note: divsufsort() requires a fixed amount of temporary space.  The
630          * implementation of divsufsort() has been modified from the original to
631          * use the provided temporary space instead of allocating its own, since
632          * we don't want to have to deal with malloc() failures here.  */
633         divsufsort(T, SA, n, tmp);
634 }
635
636 /*
637  * Build the inverse suffix array ISA from the suffix array SA.
638  *
639  * Whereas the suffix array is a mapping from suffix rank to suffix position,
640  * the inverse suffix array is a mapping from suffix position to suffix rank.
641  */
642 static void
643 build_ISA(u32 ISA[restrict], const u32 SA[restrict], u32 n)
644 {
645         for (u32 r = 0; r < n; r++)
646                 ISA[SA[r]] = r;
647 }
648
649 /*
650  * Prepare the LCP-interval tree matchfinder for a new input buffer.
651  *
652  * @mf - the initialized matchfinder structure
653  * @T - the input buffer
654  * @n - size of the input buffer in bytes.  This must be nonzero and can be at
655  *      most the max_bufsize with which lcpit_matchfinder_init() was called.
656  */
657 void
658 lcpit_matchfinder_load_buffer(struct lcpit_matchfinder *mf, const u8 *T, u32 n)
659 {
660         /* intervals[] temporarily stores SA and LCP packed together.
661          * pos_data[] temporarily stores ISA.
662          * pos_data[] is also used as the temporary space for divsufsort().  */
663
664         build_SA(mf->intervals, T, n, mf->pos_data);
665         build_ISA(mf->pos_data, mf->intervals, n);
666         if (n <= MAX_NORMAL_BUFSIZE) {
667                 for (u32 i = 0; i < PREFETCH_SAFETY; i++) {
668                         mf->intervals[n + i] = 0;
669                         mf->pos_data[n + i] = 0;
670                 }
671                 build_LCP(mf->intervals, mf->pos_data, T, n,
672                           mf->min_match_len, mf->nice_match_len);
673                 build_LCPIT(mf->intervals, mf->pos_data, n);
674                 mf->huge_mode = false;
675         } else {
676                 for (u32 i = 0; i < PREFETCH_SAFETY; i++) {
677                         mf->intervals64[n + i] = 0;
678                         mf->pos_data[n + i] = 0;
679                 }
680                 expand_SA(mf->intervals, n);
681                 build_LCP_huge(mf->intervals64, mf->pos_data, T, n,
682                                mf->min_match_len, mf->nice_match_len);
683                 build_LCPIT_huge(mf->intervals64, mf->pos_data, n);
684                 mf->huge_mode = true;
685         }
686         mf->cur_pos = 0; /* starting at beginning of input buffer  */
687         for (u32 i = 0; i < ARRAY_LEN(mf->next); i++)
688                 mf->next[i] = 0;
689 }
690
691 /*
692  * Retrieve a list of matches with the next position.
693  *
694  * The matches will be recorded in the @matches array, ordered by strictly
695  * decreasing length and strictly decreasing offset.
696  *
697  * The return value is the number of matches found and written to @matches.
698  * This can be any value in [0, nice_match_len - min_match_len + 1].
699  */
700 u32
701 lcpit_matchfinder_get_matches(struct lcpit_matchfinder *mf,
702                               struct lz_match *matches)
703 {
704         if (mf->huge_mode)
705                 return lcpit_advance_one_byte_huge(mf->cur_pos++, mf->pos_data,
706                                                    mf->intervals64, mf->next,
707                                                    matches, true);
708         else
709                 return lcpit_advance_one_byte(mf->cur_pos++, mf->pos_data,
710                                               mf->intervals, mf->next,
711                                               matches, true);
712 }
713
714 /*
715  * Skip the next @count bytes (don't search for matches at them).  @count is
716  * assumed to be > 0.
717  */
718 void
719 lcpit_matchfinder_skip_bytes(struct lcpit_matchfinder *mf, u32 count)
720 {
721         if (mf->huge_mode) {
722                 do {
723                         lcpit_advance_one_byte_huge(mf->cur_pos++, mf->pos_data,
724                                                     mf->intervals64, mf->next,
725                                                     NULL, false);
726                 } while (--count);
727         } else {
728                 do {
729                         lcpit_advance_one_byte(mf->cur_pos++, mf->pos_data,
730                                                mf->intervals, mf->next,
731                                                NULL, false);
732                 } while (--count);
733         }
734 }
735
736 /*
737  * Destroy an LCP-interval tree matchfinder that was previously initialized with
738  * lcpit_matchfinder_init().
739  */
740 void
741 lcpit_matchfinder_destroy(struct lcpit_matchfinder *mf)
742 {
743         FREE(mf->pos_data);
744         FREE(mf->intervals);
745 }