33acd61237022b56699f5f1bc07e3ffa81af3c18
[wimlib] / src / lz_sarray.c
1 /*
2  * lz_sarray.c
3  *
4  * Suffix array match-finder for Lempel-Ziv compression.
5  */
6
7 /*
8  * Copyright (c) 2013, 2014 Eric Biggers.  All rights reserved.
9  *
10  * Redistribution and use in source and binary forms, with or without
11  * modification, are permitted provided that the following conditions
12  * are met:
13  *
14  * 1. Redistributions of source code must retain the above copyright
15  *    notice, this list of conditions and the following disclaimer.
16  *
17  * 2. Redistributions in binary form must reproduce the above copyright
18  *    notice, this list of conditions and the following disclaimer in the
19  *    documentation and/or other materials provided with the distribution.
20  *
21  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
24  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE
25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
28  * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
29  * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
30  * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
31  * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32  */
33
34 #ifdef HAVE_CONFIG_H
35 #  include "config.h"
36 #endif
37
38 #include "wimlib/lz_sarray.h"
39 #include "wimlib/util.h"
40 #include "divsufsort/divsufsort.h"
41 #include <string.h>
42
43 #define DIVSUFSORT_TMP1_SIZE (256 * sizeof(saidx_t))
44 #define DIVSUFSORT_TMP2_SIZE (256 * 256 * sizeof(saidx_t))
45
46 /* If ENABLE_LZ_DEBUG is defined, verify that the suffix array satisfies its
47  * definition.
48  *
49  * @SA          The constructed suffix array.
50  * @T           The original data.
51  * @found       Temporary 'bool' array of length @n.
52  * @n           Length of the data (length of @SA, @T, and @found arrays).
53  *
54  * WARNING: this is for debug use only as it does not necessarily run in linear
55  * time!!!  */
56 static void
57 verify_suffix_array(const lz_sarray_pos_t SA[restrict],
58                     const u8 T[restrict],
59                     bool found[restrict],
60                     lz_sarray_pos_t n)
61 {
62 #ifdef ENABLE_LZ_DEBUG
63         /* Ensure the SA contains exactly one of each i in [0, n - 1].  */
64         for (lz_sarray_pos_t i = 0; i < n; i++)
65                 found[i] = false;
66         for (lz_sarray_pos_t r = 0; r < n; r++) {
67                 lz_sarray_pos_t i = SA[r];
68                 LZ_ASSERT(i < n);
69                 LZ_ASSERT(!found[i]);
70                 found[i] = true;
71         }
72
73         /* Ensure the suffix with rank r is lexicographically lesser than the
74          * suffix with rank (r + 1) for all r in [0, n - 2].  */
75         for (lz_sarray_pos_t r = 0; r < n - 1; r++) {
76
77                 lz_sarray_pos_t i1 = SA[r];
78                 lz_sarray_pos_t i2 = SA[r + 1];
79
80                 lz_sarray_pos_t n1 = n - i1;
81                 lz_sarray_pos_t n2 = n - i2;
82
83                 int res = memcmp(&T[i1], &T[i2], min(n1, n2));
84                 LZ_ASSERT(res < 0 || (res == 0 && n1 < n2));
85         }
86 #endif /* ENABLE_LZ_DEBUG  */
87 }
88
89 /* Compute the inverse suffix array @ISA from the suffix array @SA in linear
90  * time.
91  *
92  * Whereas the suffix array is a mapping from suffix rank to suffix position,
93  * the inverse suffix array is a mapping from suffix position to suffix rank.
94  */
95 static void
96 compute_inverse_suffix_array(lz_sarray_pos_t ISA[restrict],
97                              const lz_sarray_pos_t SA[restrict],
98                              lz_sarray_pos_t n)
99 {
100         lz_sarray_pos_t r;
101
102         for (r = 0; r < n; r++)
103                 ISA[SA[r]] = r;
104 }
105
106
107 /* Compute the LCP (Longest Common Prefix) array in linear time.
108  *
109  * LCP[r] will be the length of the longest common prefix between the suffixes
110  * with positions SA[r - 1] and  SA[r].  LCP[0] will be undefined.
111  *
112  * Algorithm adapted from Kasai et al. 2001: "Linear-Time Longest-Common-Prefix
113  * Computation in Suffix Arrays and Its Applications".  Modified slightly to
114  * take into account that with bytes in the real world, there is no unique
115  * symbol at the end of the string.  */
116 static void
117 compute_lcp_array(lz_sarray_pos_t LCP[restrict],
118                   const lz_sarray_pos_t SA[restrict],
119                   const lz_sarray_pos_t ISA[restrict],
120                   const u8 T[restrict],
121                   lz_sarray_pos_t n)
122 {
123         lz_sarray_pos_t h, i, r, j, lim;
124
125         h = 0;
126         for (i = 0; i < n; i++) {
127                 r = ISA[i];
128                 if (r > 0) {
129                         j = SA[r - 1];
130                         lim = min(n - i, n - j);
131
132                         while (h < lim && T[i + h] == T[j + h])
133                                 h++;
134                         LCP[r] = h;
135                         if (h > 0)
136                                 h--;
137                 }
138         }
139 }
140
141 /* If ENABLE_LZ_DEBUG is defined, verify that the LCP (Longest Common Prefix)
142  * array satisfies its definition.
143  *
144  * WARNING: this is for debug use only as it does not necessarily run in linear
145  * time!!!  */
146 static void
147 verify_lcp_array(lz_sarray_pos_t LCP[restrict],
148                  const lz_sarray_pos_t SA[restrict],
149                  const u8 T[restrict],
150                  lz_sarray_pos_t n)
151 {
152 #ifdef ENABLE_LZ_DEBUG
153         for (lz_sarray_pos_t r = 0; r < n - 1; r++) {
154                 lz_sarray_pos_t i1 = SA[r];
155                 lz_sarray_pos_t i2 = SA[r + 1];
156                 lz_sarray_pos_t lcp = LCP[r + 1];
157
158                 lz_sarray_pos_t n1 = n - i1;
159                 lz_sarray_pos_t n2 = n - i2;
160
161                 LZ_ASSERT(lcp <= min(n1, n2));
162
163                 LZ_ASSERT(memcmp(&T[i1], &T[i2], lcp) == 0);
164                 if (lcp < min(n1, n2))
165                         LZ_ASSERT(T[i1 + lcp] != T[i2 + lcp]);
166         }
167 #endif /* ENABLE_LZ_DEBUG */
168 }
169
170 /* Initialize the SA link array in linear time.
171  *
172  * This is similar to computing the LPF (Longest Previous Factor) array, which
173  * is addressed in several papers.  In particular the algorithms below are based
174  * on Crochemore et al. 2009: "LPF computation revisited".  However, this
175  * match-finder does not actually compute or use the LPF array per se.  Rather,
176  * this function sets up some information necessary to compute the LPF array,
177  * but later lz_sarray_get_matches() actually uses this information to search
178  * the suffix array directly and can keep searching beyond the first (longest)
179  * match whose length would be placed in the LPF array.  This difference from
180  * the theoretical work is necessary because in many real compression formats
181  * matches take variable numbers of bits to encode, so a decent parser needs to
182  * consider more than just the longest match with unspecified offset.
183  *
184  * Note: We cap the lcpprev and lcpnext values to the maximum match length so
185  * that the match-finder need not worry about it later, in the inner loop.
186  *
187  * Note: the LCP array is one of the inputs to this function, but it is used as
188  * temporary space and therefore will be invalidated.
189  */
190 static void
191 init_salink(struct salink link[restrict],
192             lz_sarray_pos_t LCP[restrict],
193             const lz_sarray_pos_t SA[restrict],
194             const u8 T[restrict],
195             lz_sarray_pos_t n,
196             lz_sarray_len_t min_match_len,
197             lz_sarray_len_t max_match_len)
198 {
199         /* Calculate salink.dist_to_next and salink.lcpnext.
200          *
201          * Pass 1 calculates, for each suffix rank, the corresponding
202          * "next_initial" value which is the smallest larger rank that
203          * corresponds to a suffix starting earlier in the string.  It also
204          * calculates "lcpnext_initial", which is the number of bytes shared
205          * with that suffix, although to eliminate checks in
206          * lz_sarray_get_matches(), "lcpnext_initial" is set to 0 if it's less
207          * than the minimum match length or set to the maximum match length if
208          * it's greater than the maximum match length.
209          *
210          * Pass 2 translates each absolute "next_initial", a 4-byte value, into
211          * a relative "dist_to_next", a 1-byte value.  This is done to save
212          * memory.  In the case that the true relative distance cannot be
213          * encoded in 1 byte, it is capped to 255.  This is valid as long as
214          * lz_sarray_get_matches() validates each position before using it.
215          * Note that "lcpnext" need not be updated in this case because it will
216          * not be used until the actual next rank has been found anyway.
217          */
218         link[n - 1].next_initial = LZ_SARRAY_POS_MAX;
219         link[n - 1].lcpnext_initial = 0;
220         for (lz_sarray_pos_t r = n - 2; r != LZ_SARRAY_POS_MAX; r--) {
221                 lz_sarray_pos_t t = r + 1;
222                 lz_sarray_pos_t l = LCP[t];
223                 while (t != LZ_SARRAY_POS_MAX && SA[t] > SA[r]) {
224                         l = min(l, link[t].lcpnext_initial);
225                         t = link[t].next_initial;
226                 }
227                 link[r].next_initial = t;
228
229                 if (l < min_match_len)
230                         l = 0;
231                 else if (l > max_match_len)
232                         l = max_match_len;
233                 link[r].lcpnext_initial = l;
234         }
235         for (lz_sarray_pos_t r = 0; r < n; r++) {
236                 lz_sarray_pos_t next;
237                 lz_sarray_len_t l;
238                 lz_sarray_delta_t dist_to_next;
239
240                 next = link[r].next_initial;
241                 l = link[r].lcpnext_initial;
242
243                 if (next == LZ_SARRAY_POS_MAX)
244                         dist_to_next = 0;
245                 else if (next - r <= LZ_SARRAY_DELTA_MAX)
246                         dist_to_next = next - r;
247                 else
248                         dist_to_next = LZ_SARRAY_DELTA_MAX;
249
250                 link[r].lcpnext = l;
251                 link[r].dist_to_next = dist_to_next;
252         }
253
254         /* Calculate salink.dist_to_prev and salink.lcpprev.
255          *
256          * This is analgous to dist_to_next and lcpnext as described above, but
257          * in the other direction.  That is, here we're interested in, for each
258          * rank, the largest smaller rank that corresponds to a suffix starting
259          * earlier in the string.
260          *
261          * To save memory we don't have a "prev_initial" field, but rather store
262          * those values in the LCP array.
263          */
264         LCP[0] = LZ_SARRAY_POS_MAX;
265         link[0].lcpprev = 0;
266         for (lz_sarray_pos_t r = 1; r < n; r++) {
267                 lz_sarray_pos_t t = r - 1;
268                 lz_sarray_pos_t l = LCP[r];
269                 while (t != LZ_SARRAY_POS_MAX && SA[t] > SA[r]) {
270                         l = min(l, link[t].lcpprev);
271                         t = LCP[t];
272                 }
273                 LCP[r] = t;
274
275                 if (l < min_match_len)
276                         l = 0;
277                 else if (l > max_match_len)
278                         l = max_match_len;
279
280                 link[r].lcpprev = l;
281         }
282         for (lz_sarray_pos_t r = 0; r < n; r++) {
283
284                 lz_sarray_pos_t prev = LCP[r];
285
286                 if (prev == LZ_SARRAY_POS_MAX)
287                         link[r].dist_to_prev = 0;
288                 else if (r - prev <= LZ_SARRAY_DELTA_MAX)
289                         link[r].dist_to_prev = r - prev;
290                 else
291                         link[r].dist_to_prev = LZ_SARRAY_DELTA_MAX;
292         }
293 }
294
295 /* If ENABLE_LZ_DEBUG is defined, verify the values computed by init_salink().
296  *
297  * WARNING: this is for debug use only as it does not necessarily run in linear
298  * time!!!  */
299 static void
300 verify_salink(const struct salink link[],
301               const lz_sarray_pos_t SA[],
302               const u8 T[],
303               lz_sarray_pos_t n,
304               lz_sarray_len_t min_match_len,
305               lz_sarray_len_t max_match_len)
306 {
307 #ifdef ENABLE_LZ_DEBUG
308         for (lz_sarray_pos_t r = 0; r < n; r++) {
309                 for (lz_sarray_pos_t prev = r; ; ) {
310                         if (prev == 0) {
311                                 LZ_ASSERT(link[r].dist_to_prev == 0);
312                                 LZ_ASSERT(link[r].lcpprev == 0);
313                                 break;
314                         }
315
316                         prev--;
317
318                         if (SA[prev] < SA[r]) {
319                                 LZ_ASSERT(link[r].dist_to_prev == min(r - prev, LZ_SARRAY_DELTA_MAX));
320
321                                 lz_sarray_pos_t lcpprev;
322                                 for (lcpprev = 0;
323                                      lcpprev < min(n - SA[prev], n - SA[r]) &&
324                                              T[SA[prev] + lcpprev] == T[SA[r] + lcpprev];
325                                      lcpprev++)
326                                         ;
327                                 if (lcpprev < min_match_len)
328                                         lcpprev = 0;
329                                 else if (lcpprev > max_match_len)
330                                         lcpprev = max_match_len;
331
332                                 LZ_ASSERT(lcpprev == link[r].lcpprev);
333                                 break;
334                         }
335                 }
336
337                 for (lz_sarray_pos_t next = r; ; ) {
338                         if (next == n - 1) {
339                                 LZ_ASSERT(link[r].dist_to_next == 0);
340                                 LZ_ASSERT(link[r].lcpnext == 0);
341                                 break;
342                         }
343
344                         next++;
345
346                         if (SA[next] < SA[r]) {
347                                 LZ_ASSERT(link[r].dist_to_next == min(next - r, LZ_SARRAY_DELTA_MAX));
348
349                                 lz_sarray_pos_t lcpnext;
350                                 for (lcpnext = 0;
351                                      lcpnext < min(n - SA[next], n - SA[r]) &&
352                                              T[SA[next] + lcpnext] == T[SA[r] + lcpnext];
353                                      lcpnext++)
354                                         ;
355                                 if (lcpnext < min_match_len)
356                                         lcpnext = 0;
357                                 else if (lcpnext > max_match_len)
358                                         lcpnext = max_match_len;
359
360                                 LZ_ASSERT(lcpnext == link[r].lcpnext);
361                                 break;
362                         }
363                 }
364         }
365 #endif
366 }
367
368 /*
369  * Initialize the suffix array match-finder.
370  *
371  * @mf
372  *      The suffix array match-finder structure to initialize.  This structure
373  *      is expected to be zeroed before this function is called.  In the case
374  *      that this function fails, lz_sarray_destroy() should be called to free
375  *      any memory that may have been allocated.
376  *
377  * @max_window_size
378  *      The maximum window size to support.  This must be greater than 0.
379  *
380  *      The amount of needed memory will depend on this value; see
381  *      lz_sarray_get_needed_memory() for details.
382  *
383  * @min_match_len
384  *      The minimum length of each match to be found.  Must be greater than 0.
385  *
386  * @max_match_len
387  *      The maximum length of each match to be found.  Must be greater than or
388  *      equal to @min_match_len.
389  *
390  * @max_matches_to_consider
391  *      The maximum number of matches to consider at each position.  This should
392  *      be greater than @max_matches_to_return because @max_matches_to_consider
393  *      counts all the returned matches as well as matches of equal length to
394  *      returned matches that were not returned.  This parameter bounds the
395  *      amount of work the match-finder does at any one position.  This could be
396  *      anywhere from 1 to 100+ depending on the compression ratio and
397  *      performance desired.
398  *
399  * @max_matches_to_return
400  *      Maximum number of matches to return at each position.  Because of the
401  *      suffix array search algorithm, the order in which matches are returned
402  *      will be from longest to shortest, so cut-offs due to this parameter will
403  *      only result in shorter matches being discarded.  This parameter could be
404  *      anywhere from 1 to (@max_match_len - @min_match_len + 1) depending on
405  *      the compression performance desired.  However, making it even moderately
406  *      large (say, greater than 3) may not be very helpful due to the property
407  *      that the matches are returned from longest to shortest.  But the main
408  *      thing to keep in mind is that if the compressor decides to output a
409  *      shorter-than-possible match, ideally it would be best to choose the best
410  *      match of the desired length rather than truncate a longer match to that
411  *      length.
412  *
413  * After initialization, the suffix-array match-finder can be used for any
414  * number of input strings (windows) of length less than or equal to
415  * @max_window_size by successive calls to lz_sarray_load_window().
416  *
417  * Returns %true on success, or %false if sufficient memory could not be
418  * allocated.  See the note for @max_window_size above regarding the needed
419  * memory size.
420  */
421 bool
422 lz_sarray_init(struct lz_sarray *mf,
423                lz_sarray_pos_t max_window_size,
424                lz_sarray_len_t min_match_len,
425                lz_sarray_len_t max_match_len,
426                u32 max_matches_to_consider,
427                u32 max_matches_to_return)
428 {
429         LZ_ASSERT(min_match_len > 0);
430         LZ_ASSERT(max_window_size > 0);
431         LZ_ASSERT(max_match_len >= min_match_len);
432
433         mf->max_window_size = max_window_size;
434         mf->min_match_len = min_match_len;
435         mf->max_match_len = max_match_len;
436         mf->max_matches_to_consider = max_matches_to_consider;
437         mf->max_matches_to_return = max_matches_to_return;
438
439         /* SA and ISA will share the same storage block.  */
440         if ((u64)2 * max_window_size * sizeof(mf->SA[0]) !=
441                  2 * max_window_size * sizeof(mf->SA[0]))
442                 return false;
443         mf->SA = MALLOC(max_window_size * sizeof(mf->SA[0]) +
444                         max(DIVSUFSORT_TMP1_SIZE,
445                             max_window_size * sizeof(mf->SA[0])));
446         if (mf->SA == NULL)
447                 return false;
448
449         if ((u64)max_window_size * sizeof(mf->salink[0]) !=
450                  max_window_size * sizeof(mf->salink[0]))
451                 return false;
452         mf->salink = MALLOC(max(DIVSUFSORT_TMP2_SIZE,
453                                 max_window_size * sizeof(mf->salink[0])));
454         if (mf->salink == NULL)
455                 return false;
456
457         return true;
458 }
459
460 /*
461  * Return the number of bytes of memory that lz_sarray_init() would allocate for
462  * the specified maximum window size.
463  *
464  * This should be (14 * @max_window_size) unless the type definitions have been
465  * changed.
466  */
467 u64
468 lz_sarray_get_needed_memory(lz_sarray_pos_t max_window_size)
469 {
470         u64 size = 0;
471
472         /* SA and ISA: 8 bytes per position  */
473         size += (u64)max_window_size * sizeof(((struct lz_sarray*)0)->SA[0]) +
474                 max(DIVSUFSORT_TMP1_SIZE,
475                     (u64)max_window_size * sizeof(((struct lz_sarray*)0)->SA[0]));
476
477         /* salink: 6 bytes per position  */
478         size += max(DIVSUFSORT_TMP2_SIZE,
479                     (u64)max_window_size * sizeof(((struct lz_sarray*)0)->salink[0]));
480
481         return size;
482 }
483
484 /*
485  * Prepare the suffix array match-finder to scan the specified window for
486  * matches.
487  *
488  * @mf  Suffix array match-finder previously initialized with lz_sarray_init().
489  *
490  * @T   Window, or "block", in which to find matches.
491  *
492  * @n   Size of window in bytes.  This must be positive and less than or equal
493  *      to the @max_window_size passed to lz_sarray_init().
494  *
495  * This function runs in linear time (relative to @n).
496  */
497 void
498 lz_sarray_load_window(struct lz_sarray *mf, const u8 T[], lz_sarray_pos_t n)
499 {
500         lz_sarray_pos_t *ISA, *LCP;
501
502         LZ_ASSERT(n > 0 && n <= mf->max_window_size);
503
504         /* Compute SA (Suffix Array).
505          *
506          * divsufsort() needs temporary space --- one array with 256 spaces and
507          * one array with 65536 spaces.  The implementation of divsufsort() has
508          * been modified from the original to use the provided temporary space
509          * instead of allocating its own.
510          *
511          * We also check at build-time that divsufsort() uses the same integer
512          * size expected by this code.  Unfortunately, divsufsort breaks if
513          * 'sa_idx_t' is defined to be a 16-bit integer; however, that would
514          * limit blocks to only 65536 bytes anyway.  */
515         BUILD_BUG_ON(sizeof(lz_sarray_pos_t) != sizeof(saidx_t));
516
517         divsufsort(T, mf->SA, n, (saidx_t*)&mf->SA[n], (saidx_t*)mf->salink);
518
519         BUILD_BUG_ON(sizeof(bool) > sizeof(mf->salink[0]));
520         verify_suffix_array(mf->SA, T, (bool*)mf->salink, n);
521
522         /* Compute ISA (Inverse Suffix Array) in a preliminary position.
523          *
524          * This is just a trick to save memory.  Since LCP is unneeded after
525          * this function, it can be computed in any available space.  The
526          * storage for the ISA is the best choice because the ISA can be built
527          * quickly in salink for now, then re-built in its real location at the
528          * end.  This is probably worth it because computing the ISA from the SA
529          * is very fast, and since this match-finder is memory-hungry we'd like
530          * to save as much memory as possible.  */
531         BUILD_BUG_ON(sizeof(mf->salink[0]) < sizeof(mf->ISA[0]));
532         ISA = (lz_sarray_pos_t*)mf->salink;
533         compute_inverse_suffix_array(ISA, mf->SA, n);
534
535         /* Compute LCP (Longest Common Prefix) array.  */
536         LCP = mf->SA + n;
537         compute_lcp_array(LCP, mf->SA, ISA, T, n);
538         verify_lcp_array(LCP, mf->SA, T, n);
539
540         /* Initialize suffix array links.  */
541         init_salink(mf->salink, LCP, mf->SA, T, n,
542                     mf->min_match_len, mf->max_match_len);
543         verify_salink(mf->salink, mf->SA, T, n,
544                       mf->min_match_len, mf->max_match_len);
545
546         /* Compute ISA (Inverse Suffix Array) in its final position.  */
547         ISA = mf->SA + n;
548         compute_inverse_suffix_array(ISA, mf->SA, n);
549
550         /* Save new variables and return.  */
551         mf->ISA = ISA;
552         mf->cur_pos = 0;
553         mf->window_size = n;
554 }
555
556 /* Free memory allocated for the suffix array match-finder.  */
557 void
558 lz_sarray_destroy(struct lz_sarray *mf)
559 {
560         FREE(mf->SA);
561         FREE(mf->salink);
562 }