From 4f6a8f018505317a116447d1fff83451f2d79808 Mon Sep 17 00:00:00 2001 From: ddidderr Date: Sat, 11 Jul 2026 14:21:09 +0200 Subject: [PATCH] feat(rust): port divsufsort Move the dictionary builder's suffix-array construction from lib/dictBuilder/divsufsort.c to rust/src/divsufsort.rs, the first dictBuilder module to migrate. It rides on the dict-builder cargo feature dimension introduced by the previous commit. divsufsort() is a self-contained algorithm (two-stage sort of type-B* substrings via sssort, rank refinement via trsort, then induced sorting of the full array), so its context-free signature allows a direct symbol takeover: the Rust #[no_mangle] export provides the existing `divsufsort` symbol and the C file becomes a declaration-only shim that just keeps the header's prototypes in the build. Only divsufsort() moved; divbwt() has no callers anywhere in zstd, so it is now declaration-only, keeping the Rust export surface minimal. The unused openMP parameter is retained for signature compatibility (zstd never defines LIBBSC_OPENMP). The port is a mechanical translation of the exact configuration zstd compiles: ALPHABET_SIZE=256, SS_INSERTIONSORT_THRESHOLD=8, SS_BLOCKSIZE=1024, SS_MISORT_STACKSIZE=16, SS_SMERGE_STACKSIZE=32, TR_STACKSIZE=64. Every C `int*` cursor into the SA buffer becomes an `isize` index into a single `&mut [i32]` slice, preserving the pointer arithmetic (including transient one-before-the-range cursors and the bitwise-complement rank marking) while staying bounds-checked; all value arithmetic keeps C int semantics. The C -1/-2 error results are preserved, with Vec::try_reserve_exact standing in for the bucket-array malloc failure path. Behavior is bit-identical by construction and by measurement (see test plan); runtime on an 11 MB training buffer is within ~5% of the C build end-to-end. Users see no behavioral change: dictionaries trained through ZDICT_trainFromBuffer_legacy() are byte-identical to the C build. The only external difference is that the never-called `divbwt` symbol is no longer defined in the library. Test plan: - cd rust && cargo fmt --check && cargo clippy --all-targets -- -D warnings && cargo test --all-targets && cargo build --release (125 tests pass; new unit tests cover empty/one/two-byte inputs, all-equal bytes, an exact hand-computed "abracadabra" SA, and fixed-seed LCG buffers at 256-, 4-, and 2-symbol alphabets verified against a naive reference sort plus permutation/sorted invariants) - Feature matrix: cargo build --release --no-default-features --features compression,decompression (and decompression-only, compression-only, compression,dict-builder); `divsufsort` is exported only when dict-builder is enabled - make -C tests fuzzer && ./tests/fuzzer -i1 --no-big-tests (includes ZDICT training tests): pass - make -C tests test-rust-lib-smoke: pass - make -C tests test-invalidDictionaries: pass - make -C programs zstd zstd-dictBuilder zstd-small zstd-compress zstd-decompress: build; compress/decompress round-trip verified - Byte-identity vs pristine C build (commit 959e4852): a harness calling ZDICT_trainFromBuffer_legacy() (the only zstd path reaching divsufsort) and divsufsort() directly, linked against both libzstd.a builds, produces byte-identical dictionaries (80,288 B and full 112,640 B capacity) and byte-identical suffix arrays on a 1 MB source set and an 11 MB binary/repetitive set; a differential driver over 148 random and structured buffers (sizes 3..6000, alphabets 1..256, Fibonacci word, sawtooth, 6 KB near-constant) shows zero mismatches. The CLI --train path could not be exercised because the Rust CLI frontend rejects --train in both the pristine and ported builds (a pre-existing migration gap unrelated to this change). --- lib/dictBuilder/divsufsort.c | 1891 +---------------------- rust/README.md | 13 +- rust/src/divsufsort.rs | 2756 ++++++++++++++++++++++++++++++++++ rust/src/lib.rs | 2 + 4 files changed, 2772 insertions(+), 1890 deletions(-) create mode 100644 rust/src/divsufsort.rs diff --git a/lib/dictBuilder/divsufsort.c b/lib/dictBuilder/divsufsort.c index a2870fb3b..133c08bf4 100644 --- a/lib/dictBuilder/divsufsort.c +++ b/lib/dictBuilder/divsufsort.c @@ -24,1890 +24,9 @@ * OTHER DEALINGS IN THE SOFTWARE. */ -/*- Compiler specifics -*/ -#ifdef __clang__ -#pragma clang diagnostic ignored "-Wshorten-64-to-32" -#endif - -#if defined(_MSC_VER) -# pragma warning(disable : 4244) -# pragma warning(disable : 4127) /* C4127 : Condition expression is constant */ -#endif - - -/*- Dependencies -*/ -#include -#include -#include - +/* divsufsort() is implemented in rust/src/divsufsort.rs, which provides the + * symbol directly. This translation unit keeps the header's prototypes in + * the build so the dictionary builder continues to compile against the + * original interface. divbwt() has no callers in zstd and is declaration- + * only; it moves to Rust if a user ever appears. */ #include "divsufsort.h" - -/*- Constants -*/ -#if defined(INLINE) -# undef INLINE -#endif -#if !defined(INLINE) -# define INLINE __inline -#endif -#if defined(ALPHABET_SIZE) && (ALPHABET_SIZE < 1) -# undef ALPHABET_SIZE -#endif -#if !defined(ALPHABET_SIZE) -# define ALPHABET_SIZE (256) -#endif -#define BUCKET_A_SIZE (ALPHABET_SIZE) -#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE) -#if defined(SS_INSERTIONSORT_THRESHOLD) -# if SS_INSERTIONSORT_THRESHOLD < 1 -# undef SS_INSERTIONSORT_THRESHOLD -# define SS_INSERTIONSORT_THRESHOLD (1) -# endif -#else -# define SS_INSERTIONSORT_THRESHOLD (8) -#endif -#if defined(SS_BLOCKSIZE) -# if SS_BLOCKSIZE < 0 -# undef SS_BLOCKSIZE -# define SS_BLOCKSIZE (0) -# elif 32768 <= SS_BLOCKSIZE -# undef SS_BLOCKSIZE -# define SS_BLOCKSIZE (32767) -# endif -#else -# define SS_BLOCKSIZE (1024) -#endif -/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */ -#if SS_BLOCKSIZE == 0 -# define SS_MISORT_STACKSIZE (96) -#elif SS_BLOCKSIZE <= 4096 -# define SS_MISORT_STACKSIZE (16) -#else -# define SS_MISORT_STACKSIZE (24) -#endif -#define SS_SMERGE_STACKSIZE (32) -#define TR_INSERTIONSORT_THRESHOLD (8) -#define TR_STACKSIZE (64) - - -/*- Macros -*/ -#ifndef SWAP -# define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0) -#endif /* SWAP */ -#ifndef MIN -# define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b)) -#endif /* MIN */ -#ifndef MAX -# define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b)) -#endif /* MAX */ -#define STACK_PUSH(_a, _b, _c, _d)\ - do {\ - assert(ssize < STACK_SIZE);\ - stack[ssize].a = (_a), stack[ssize].b = (_b),\ - stack[ssize].c = (_c), stack[ssize++].d = (_d);\ - } while(0) -#define STACK_PUSH5(_a, _b, _c, _d, _e)\ - do {\ - assert(ssize < STACK_SIZE);\ - stack[ssize].a = (_a), stack[ssize].b = (_b),\ - stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\ - } while(0) -#define STACK_POP(_a, _b, _c, _d)\ - do {\ - assert(0 <= ssize);\ - if(ssize == 0) { return; }\ - (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ - (_c) = stack[ssize].c, (_d) = stack[ssize].d;\ - } while(0) -#define STACK_POP5(_a, _b, _c, _d, _e)\ - do {\ - assert(0 <= ssize);\ - if(ssize == 0) { return; }\ - (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ - (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\ - } while(0) -#define BUCKET_A(_c0) bucket_A[(_c0)] -#if ALPHABET_SIZE == 256 -#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)]) -#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)]) -#else -#define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)]) -#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)]) -#endif - - -/*- Private Functions -*/ - -static const int lg_table[256]= { - -1,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, - 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, - 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, - 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7 -}; - -#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) - -static INLINE -int -ss_ilg(int n) { -#if SS_BLOCKSIZE == 0 - return (n & 0xffff0000) ? - ((n & 0xff000000) ? - 24 + lg_table[(n >> 24) & 0xff] : - 16 + lg_table[(n >> 16) & 0xff]) : - ((n & 0x0000ff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]); -#elif SS_BLOCKSIZE < 256 - return lg_table[n]; -#else - return (n & 0xff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]; -#endif -} - -#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ - -#if SS_BLOCKSIZE != 0 - -static const int sqq_table[256] = { - 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, - 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, - 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109, -110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, -128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, -143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, -156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, -169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, -181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, -192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, -202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211, -212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, -221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, -230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, -239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, -247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255 -}; - -static INLINE -int -ss_isqrt(int x) { - int y, e; - - if(x >= (SS_BLOCKSIZE * SS_BLOCKSIZE)) { return SS_BLOCKSIZE; } - e = (x & 0xffff0000) ? - ((x & 0xff000000) ? - 24 + lg_table[(x >> 24) & 0xff] : - 16 + lg_table[(x >> 16) & 0xff]) : - ((x & 0x0000ff00) ? - 8 + lg_table[(x >> 8) & 0xff] : - 0 + lg_table[(x >> 0) & 0xff]); - - if(e >= 16) { - y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7); - if(e >= 24) { y = (y + 1 + x / y) >> 1; } - y = (y + 1 + x / y) >> 1; - } else if(e >= 8) { - y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1; - } else { - return sqq_table[x] >> 4; - } - - return (x < (y * y)) ? y - 1 : y; -} - -#endif /* SS_BLOCKSIZE != 0 */ - - -/*---------------------------------------------------------------------------*/ - -/* Compares two suffixes. */ -static INLINE -int -ss_compare(const unsigned char *T, - const int *p1, const int *p2, - int depth) { - const unsigned char *U1, *U2, *U1n, *U2n; - - for(U1 = T + depth + *p1, - U2 = T + depth + *p2, - U1n = T + *(p1 + 1) + 2, - U2n = T + *(p2 + 1) + 2; - (U1 < U1n) && (U2 < U2n) && (*U1 == *U2); - ++U1, ++U2) { - } - - return U1 < U1n ? - (U2 < U2n ? *U1 - *U2 : 1) : - (U2 < U2n ? -1 : 0); -} - - -/*---------------------------------------------------------------------------*/ - -#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) - -/* Insertionsort for small size groups */ -static -void -ss_insertionsort(const unsigned char *T, const int *PA, - int *first, int *last, int depth) { - int *i, *j; - int t; - int r; - - for(i = last - 2; first <= i; --i) { - for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) { - do { *(j - 1) = *j; } while((++j < last) && (*j < 0)); - if(last <= j) { break; } - } - if(r == 0) { *j = ~*j; } - *(j - 1) = t; - } -} - -#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */ - - -/*---------------------------------------------------------------------------*/ - -#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) - -static INLINE -void -ss_fixdown(const unsigned char *Td, const int *PA, - int *SA, int i, int size) { - int j, k; - int v; - int c, d, e; - - for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { - d = Td[PA[SA[k = j++]]]; - if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; } - if(d <= c) { break; } - } - SA[i] = v; -} - -/* Simple top-down heapsort. */ -static -void -ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) { - int i, m; - int t; - - m = size; - if((size % 2) == 0) { - m--; - if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); } - } - - for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); } - if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); } - for(i = m - 1; 0 < i; --i) { - t = SA[0], SA[0] = SA[i]; - ss_fixdown(Td, PA, SA, 0, i); - SA[i] = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Returns the median of three elements. */ -static INLINE -int * -ss_median3(const unsigned char *Td, const int *PA, - int *v1, int *v2, int *v3) { - int *t; - if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); } - if(Td[PA[*v2]] > Td[PA[*v3]]) { - if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; } - else { return v3; } - } - return v2; -} - -/* Returns the median of five elements. */ -static INLINE -int * -ss_median5(const unsigned char *Td, const int *PA, - int *v1, int *v2, int *v3, int *v4, int *v5) { - int *t; - if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); } - if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); } - if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); } - if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); } - if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); } - if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; } - return v3; -} - -/* Returns the pivot element. */ -static INLINE -int * -ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) { - int *middle; - int t; - - t = last - first; - middle = first + t / 2; - - if(t <= 512) { - if(t <= 32) { - return ss_median3(Td, PA, first, middle, last - 1); - } else { - t >>= 2; - return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1); - } - } - t >>= 3; - first = ss_median3(Td, PA, first, first + t, first + (t << 1)); - middle = ss_median3(Td, PA, middle - t, middle, middle + t); - last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1); - return ss_median3(Td, PA, first, middle, last); -} - - -/*---------------------------------------------------------------------------*/ - -/* Binary partition for substrings. */ -static INLINE -int * -ss_partition(const int *PA, - int *first, int *last, int depth) { - int *a, *b; - int t; - for(a = first - 1, b = last;;) { - for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; } - for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { } - if(b <= a) { break; } - t = ~*b; - *b = *a; - *a = t; - } - if(first < a) { *first = ~*first; } - return a; -} - -/* Multikey introsort for medium size groups. */ -static -void -ss_mintrosort(const unsigned char *T, const int *PA, - int *first, int *last, - int depth) { -#define STACK_SIZE SS_MISORT_STACKSIZE - struct { int *a, *b, c; int d; } stack[STACK_SIZE]; - const unsigned char *Td; - int *a, *b, *c, *d, *e, *f; - int s, t; - int ssize; - int limit; - int v, x = 0; - - for(ssize = 0, limit = ss_ilg(last - first);;) { - - if((last - first) <= SS_INSERTIONSORT_THRESHOLD) { -#if 1 < SS_INSERTIONSORT_THRESHOLD - if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); } -#endif - STACK_POP(first, last, depth, limit); - continue; - } - - Td = T + depth; - if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); } - if(limit < 0) { - for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) { - if((x = Td[PA[*a]]) != v) { - if(1 < (a - first)) { break; } - v = x; - first = a; - } - } - if(Td[PA[*first] - 1] < v) { - first = ss_partition(PA, first, a, depth); - } - if((a - first) <= (last - a)) { - if(1 < (a - first)) { - STACK_PUSH(a, last, depth, -1); - last = a, depth += 1, limit = ss_ilg(a - first); - } else { - first = a, limit = -1; - } - } else { - if(1 < (last - a)) { - STACK_PUSH(first, a, depth + 1, ss_ilg(a - first)); - first = a, limit = -1; - } else { - last = a, depth += 1, limit = ss_ilg(a - first); - } - } - continue; - } - - /* choose pivot */ - a = ss_pivot(Td, PA, first, last); - v = Td[PA[*a]]; - SWAP(*first, *a); - - /* partition */ - for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { } - if(((a = b) < last) && (x < v)) { - for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - } - for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { } - if((b < (d = c)) && (x > v)) { - for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - for(; b < c;) { - SWAP(*b, *c); - for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - - if(a <= d) { - c = b - 1; - - if((s = a - first) > (t = b - a)) { s = t; } - for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - if((s = d - c) > (t = last - d - 1)) { s = t; } - for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - - a = first + (b - a), c = last - (d - c); - b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth); - - if((a - first) <= (last - c)) { - if((last - c) <= (c - b)) { - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - STACK_PUSH(c, last, depth, limit); - last = a; - } else if((a - first) <= (c - b)) { - STACK_PUSH(c, last, depth, limit); - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - last = a; - } else { - STACK_PUSH(c, last, depth, limit); - STACK_PUSH(first, a, depth, limit); - first = b, last = c, depth += 1, limit = ss_ilg(c - b); - } - } else { - if((a - first) <= (c - b)) { - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - STACK_PUSH(first, a, depth, limit); - first = c; - } else if((last - c) <= (c - b)) { - STACK_PUSH(first, a, depth, limit); - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - first = c; - } else { - STACK_PUSH(first, a, depth, limit); - STACK_PUSH(c, last, depth, limit); - first = b, last = c, depth += 1, limit = ss_ilg(c - b); - } - } - } else { - limit += 1; - if(Td[PA[*first] - 1] < v) { - first = ss_partition(PA, first, last, depth); - limit = ss_ilg(last - first); - } - depth += 1; - } - } -#undef STACK_SIZE -} - -#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ - - -/*---------------------------------------------------------------------------*/ - -#if SS_BLOCKSIZE != 0 - -static INLINE -void -ss_blockswap(int *a, int *b, int n) { - int t; - for(; 0 < n; --n, ++a, ++b) { - t = *a, *a = *b, *b = t; - } -} - -static INLINE -void -ss_rotate(int *first, int *middle, int *last) { - int *a, *b, t; - int l, r; - l = middle - first, r = last - middle; - for(; (0 < l) && (0 < r);) { - if(l == r) { ss_blockswap(first, middle, l); break; } - if(l < r) { - a = last - 1, b = middle - 1; - t = *a; - do { - *a-- = *b, *b-- = *a; - if(b < first) { - *a = t; - last = a; - if((r -= l + 1) <= l) { break; } - a -= 1, b = middle - 1; - t = *a; - } - } while(1); - } else { - a = first, b = middle; - t = *a; - do { - *a++ = *b, *b++ = *a; - if(last <= b) { - *a = t; - first = a + 1; - if((l -= r + 1) <= r) { break; } - a += 1, b = middle; - t = *a; - } - } while(1); - } - } -} - - -/*---------------------------------------------------------------------------*/ - -static -void -ss_inplacemerge(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int depth) { - const int *p; - int *a, *b; - int len, half; - int q, r; - int x; - - for(;;) { - if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); } - else { x = 0; p = PA + *(last - 1); } - for(a = first, len = middle - first, half = len >> 1, r = -1; - 0 < len; - len = half, half >>= 1) { - b = a + half; - q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth); - if(q < 0) { - a = b + 1; - half -= (len & 1) ^ 1; - } else { - r = q; - } - } - if(a < middle) { - if(r == 0) { *a = ~*a; } - ss_rotate(a, middle, last); - last -= middle - a; - middle = a; - if(first == middle) { break; } - } - --last; - if(x != 0) { while(*--last < 0) { } } - if(middle == last) { break; } - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Merge-forward with internal buffer. */ -static -void -ss_mergeforward(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int depth) { - int *a, *b, *c, *bufend; - int t; - int r; - - bufend = buf + (middle - first) - 1; - ss_blockswap(buf, first, middle - first); - - for(t = *(a = first), b = buf, c = middle;;) { - r = ss_compare(T, PA + *b, PA + *c, depth); - if(r < 0) { - do { - *a++ = *b; - if(bufend <= b) { *bufend = t; return; } - *b++ = *a; - } while(*b < 0); - } else if(r > 0) { - do { - *a++ = *c, *c++ = *a; - if(last <= c) { - while(b < bufend) { *a++ = *b, *b++ = *a; } - *a = *b, *b = t; - return; - } - } while(*c < 0); - } else { - *c = ~*c; - do { - *a++ = *b; - if(bufend <= b) { *bufend = t; return; } - *b++ = *a; - } while(*b < 0); - - do { - *a++ = *c, *c++ = *a; - if(last <= c) { - while(b < bufend) { *a++ = *b, *b++ = *a; } - *a = *b, *b = t; - return; - } - } while(*c < 0); - } - } -} - -/* Merge-backward with internal buffer. */ -static -void -ss_mergebackward(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int depth) { - const int *p1, *p2; - int *a, *b, *c, *bufend; - int t; - int r; - int x; - - bufend = buf + (last - middle) - 1; - ss_blockswap(buf, middle, last - middle); - - x = 0; - if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; } - else { p1 = PA + *bufend; } - if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; } - else { p2 = PA + *(middle - 1); } - for(t = *(a = last - 1), b = bufend, c = middle - 1;;) { - r = ss_compare(T, p1, p2, depth); - if(0 < r) { - if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } - *a-- = *b; - if(b <= buf) { *buf = t; break; } - *b-- = *a; - if(*b < 0) { p1 = PA + ~*b; x |= 1; } - else { p1 = PA + *b; } - } else if(r < 0) { - if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } - *a-- = *c, *c-- = *a; - if(c < first) { - while(buf < b) { *a-- = *b, *b-- = *a; } - *a = *b, *b = t; - break; - } - if(*c < 0) { p2 = PA + ~*c; x |= 2; } - else { p2 = PA + *c; } - } else { - if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } - *a-- = ~*b; - if(b <= buf) { *buf = t; break; } - *b-- = *a; - if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } - *a-- = *c, *c-- = *a; - if(c < first) { - while(buf < b) { *a-- = *b, *b-- = *a; } - *a = *b, *b = t; - break; - } - if(*b < 0) { p1 = PA + ~*b; x |= 1; } - else { p1 = PA + *b; } - if(*c < 0) { p2 = PA + ~*c; x |= 2; } - else { p2 = PA + *c; } - } - } -} - -/* D&C based merge. */ -static -void -ss_swapmerge(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int bufsize, int depth) { -#define STACK_SIZE SS_SMERGE_STACKSIZE -#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a))) -#define MERGE_CHECK(a, b, c)\ - do {\ - if(((c) & 1) ||\ - (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\ - *(a) = ~*(a);\ - }\ - if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\ - *(b) = ~*(b);\ - }\ - } while(0) - struct { int *a, *b, *c; int d; } stack[STACK_SIZE]; - int *l, *r, *lm, *rm; - int m, len, half; - int ssize; - int check, next; - - for(check = 0, ssize = 0;;) { - if((last - middle) <= bufsize) { - if((first < middle) && (middle < last)) { - ss_mergebackward(T, PA, first, middle, last, buf, depth); - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - continue; - } - - if((middle - first) <= bufsize) { - if(first < middle) { - ss_mergeforward(T, PA, first, middle, last, buf, depth); - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - continue; - } - - for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1; - 0 < len; - len = half, half >>= 1) { - if(ss_compare(T, PA + GETIDX(*(middle + m + half)), - PA + GETIDX(*(middle - m - half - 1)), depth) < 0) { - m += half + 1; - half -= (len & 1) ^ 1; - } - } - - if(0 < m) { - lm = middle - m, rm = middle + m; - ss_blockswap(lm, middle, m); - l = r = middle, next = 0; - if(rm < last) { - if(*rm < 0) { - *rm = ~*rm; - if(first < lm) { for(; *--l < 0;) { } next |= 4; } - next |= 1; - } else if(first < lm) { - for(; *r < 0; ++r) { } - next |= 2; - } - } - - if((l - first) <= (last - r)) { - STACK_PUSH(r, rm, last, (next & 3) | (check & 4)); - middle = lm, last = l, check = (check & 3) | (next & 4); - } else { - if((next & 2) && (r == middle)) { next ^= 6; } - STACK_PUSH(first, lm, l, (check & 3) | (next & 4)); - first = r, middle = rm, check = (next & 3) | (check & 4); - } - } else { - if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) { - *middle = ~*middle; - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - } - } -#undef STACK_SIZE -} - -#endif /* SS_BLOCKSIZE != 0 */ - - -/*---------------------------------------------------------------------------*/ - -/* Substring sort */ -static -void -sssort(const unsigned char *T, const int *PA, - int *first, int *last, - int *buf, int bufsize, - int depth, int n, int lastsuffix) { - int *a; -#if SS_BLOCKSIZE != 0 - int *b, *middle, *curbuf; - int j, k, curbufsize, limit; -#endif - int i; - - if(lastsuffix != 0) { ++first; } - -#if SS_BLOCKSIZE == 0 - ss_mintrosort(T, PA, first, last, depth); -#else - if((bufsize < SS_BLOCKSIZE) && - (bufsize < (last - first)) && - (bufsize < (limit = ss_isqrt(last - first)))) { - if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; } - buf = middle = last - limit, bufsize = limit; - } else { - middle = last, limit = 0; - } - for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) { -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth); -#endif - curbufsize = last - (a + SS_BLOCKSIZE); - curbuf = a + SS_BLOCKSIZE; - if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; } - for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) { - ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth); - } - } -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, a, middle, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, a, middle, depth); -#endif - for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) { - if(i & 1) { - ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth); - a -= k; - } - } - if(limit != 0) { -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, middle, last, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, middle, last, depth); -#endif - ss_inplacemerge(T, PA, first, middle, last, depth); - } -#endif - - if(lastsuffix != 0) { - /* Insert last type B* suffix. */ - int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2; - for(a = first, i = *(first - 1); - (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth))); - ++a) { - *(a - 1) = *a; - } - *(a - 1) = i; - } -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -int -tr_ilg(int n) { - return (n & 0xffff0000) ? - ((n & 0xff000000) ? - 24 + lg_table[(n >> 24) & 0xff] : - 16 + lg_table[(n >> 16) & 0xff]) : - ((n & 0x0000ff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]); -} - - -/*---------------------------------------------------------------------------*/ - -/* Simple insertionsort for small size groups. */ -static -void -tr_insertionsort(const int *ISAd, int *first, int *last) { - int *a, *b; - int t, r; - - for(a = first + 1; a < last; ++a) { - for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) { - do { *(b + 1) = *b; } while((first <= --b) && (*b < 0)); - if(b < first) { break; } - } - if(r == 0) { *b = ~*b; } - *(b + 1) = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -void -tr_fixdown(const int *ISAd, int *SA, int i, int size) { - int j, k; - int v; - int c, d, e; - - for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { - d = ISAd[SA[k = j++]]; - if(d < (e = ISAd[SA[j]])) { k = j; d = e; } - if(d <= c) { break; } - } - SA[i] = v; -} - -/* Simple top-down heapsort. */ -static -void -tr_heapsort(const int *ISAd, int *SA, int size) { - int i, m; - int t; - - m = size; - if((size % 2) == 0) { - m--; - if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); } - } - - for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); } - if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); } - for(i = m - 1; 0 < i; --i) { - t = SA[0], SA[0] = SA[i]; - tr_fixdown(ISAd, SA, 0, i); - SA[i] = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Returns the median of three elements. */ -static INLINE -int * -tr_median3(const int *ISAd, int *v1, int *v2, int *v3) { - int *t; - if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); } - if(ISAd[*v2] > ISAd[*v3]) { - if(ISAd[*v1] > ISAd[*v3]) { return v1; } - else { return v3; } - } - return v2; -} - -/* Returns the median of five elements. */ -static INLINE -int * -tr_median5(const int *ISAd, - int *v1, int *v2, int *v3, int *v4, int *v5) { - int *t; - if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); } - if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); } - if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4); SWAP(v3, v5); } - if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3); } - if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4); SWAP(v3, v5); } - if(ISAd[*v3] > ISAd[*v4]) { return v4; } - return v3; -} - -/* Returns the pivot element. */ -static INLINE -int * -tr_pivot(const int *ISAd, int *first, int *last) { - int *middle; - int t; - - t = last - first; - middle = first + t / 2; - - if(t <= 512) { - if(t <= 32) { - return tr_median3(ISAd, first, middle, last - 1); - } else { - t >>= 2; - return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1); - } - } - t >>= 3; - first = tr_median3(ISAd, first, first + t, first + (t << 1)); - middle = tr_median3(ISAd, middle - t, middle, middle + t); - last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1); - return tr_median3(ISAd, first, middle, last); -} - - -/*---------------------------------------------------------------------------*/ - -typedef struct _trbudget_t trbudget_t; -struct _trbudget_t { - int chance; - int remain; - int incval; - int count; -}; - -static INLINE -void -trbudget_init(trbudget_t *budget, int chance, int incval) { - budget->chance = chance; - budget->remain = budget->incval = incval; -} - -static INLINE -int -trbudget_check(trbudget_t *budget, int size) { - if(size <= budget->remain) { budget->remain -= size; return 1; } - if(budget->chance == 0) { budget->count += size; return 0; } - budget->remain += budget->incval - size; - budget->chance -= 1; - return 1; -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -void -tr_partition(const int *ISAd, - int *first, int *middle, int *last, - int **pa, int **pb, int v) { - int *a, *b, *c, *d, *e, *f; - int t, s; - int x = 0; - - for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { } - if(((a = b) < last) && (x < v)) { - for(; (++b < last) && ((x = ISAd[*b]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - } - for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { } - if((b < (d = c)) && (x > v)) { - for(; (b < --c) && ((x = ISAd[*c]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - for(; b < c;) { - SWAP(*b, *c); - for(; (++b < c) && ((x = ISAd[*b]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - for(; (b < --c) && ((x = ISAd[*c]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - - if(a <= d) { - c = b - 1; - if((s = a - first) > (t = b - a)) { s = t; } - for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - if((s = d - c) > (t = last - d - 1)) { s = t; } - for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - first += (b - a), last -= (d - c); - } - *pa = first, *pb = last; -} - -static -void -tr_copy(int *ISA, const int *SA, - int *first, int *a, int *b, int *last, - int depth) { - /* sort suffixes of middle partition - by using sorted order of suffixes of left and right partition. */ - int *c, *d, *e; - int s, v; - - v = b - SA - 1; - for(c = first, d = a - 1; c <= d; ++c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *++d = s; - ISA[s] = d - SA; - } - } - for(c = last - 1, e = d + 1, d = b; e < d; --c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *--d = s; - ISA[s] = d - SA; - } - } -} - -static -void -tr_partialcopy(int *ISA, const int *SA, - int *first, int *a, int *b, int *last, - int depth) { - int *c, *d, *e; - int s, v; - int rank, lastrank, newrank = -1; - - v = b - SA - 1; - lastrank = -1; - for(c = first, d = a - 1; c <= d; ++c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *++d = s; - rank = ISA[s + depth]; - if(lastrank != rank) { lastrank = rank; newrank = d - SA; } - ISA[s] = newrank; - } - } - - lastrank = -1; - for(e = d; first <= e; --e) { - rank = ISA[*e]; - if(lastrank != rank) { lastrank = rank; newrank = e - SA; } - if(newrank != rank) { ISA[*e] = newrank; } - } - - lastrank = -1; - for(c = last - 1, e = d + 1, d = b; e < d; --c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *--d = s; - rank = ISA[s + depth]; - if(lastrank != rank) { lastrank = rank; newrank = d - SA; } - ISA[s] = newrank; - } - } -} - -static -void -tr_introsort(int *ISA, const int *ISAd, - int *SA, int *first, int *last, - trbudget_t *budget) { -#define STACK_SIZE TR_STACKSIZE - struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE]; - int *a, *b, *c; - int t; - int v, x = 0; - int incr = ISAd - ISA; - int limit, next; - int ssize, trlink = -1; - - for(ssize = 0, limit = tr_ilg(last - first);;) { - - if(limit < 0) { - if(limit == -1) { - /* tandem repeat partition */ - tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1); - - /* update ranks */ - if(a < last) { - for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } - } - if(b < last) { - for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } - } - - /* push */ - if(1 < (b - a)) { - STACK_PUSH5(NULL, a, b, 0, 0); - STACK_PUSH5(ISAd - incr, first, last, -2, trlink); - trlink = ssize - 2; - } - if((a - first) <= (last - b)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink); - last = a, limit = tr_ilg(a - first); - } else if(1 < (last - b)) { - first = b, limit = tr_ilg(last - b); - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } else { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink); - first = b, limit = tr_ilg(last - b); - } else if(1 < (a - first)) { - last = a, limit = tr_ilg(a - first); - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } else if(limit == -2) { - /* tandem repeat copy */ - a = stack[--ssize].b, b = stack[ssize].c; - if(stack[ssize].d == 0) { - tr_copy(ISA, SA, first, a, b, last, ISAd - ISA); - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA); - } - STACK_POP5(ISAd, first, last, limit, trlink); - } else { - /* sorted partition */ - if(0 <= *first) { - a = first; - do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a)); - first = a; - } - if(first < last) { - a = first; do { *a = ~*a; } while(*++a < 0); - next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1; - if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } } - - /* push */ - if(trbudget_check(budget, a - first)) { - if((a - first) <= (last - a)) { - STACK_PUSH5(ISAd, a, last, -3, trlink); - ISAd += incr, last = a, limit = next; - } else { - if(1 < (last - a)) { - STACK_PUSH5(ISAd + incr, first, a, next, trlink); - first = a, limit = -3; - } else { - ISAd += incr, last = a, limit = next; - } - } - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - if(1 < (last - a)) { - first = a, limit = -3; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - continue; - } - - if((last - first) <= TR_INSERTIONSORT_THRESHOLD) { - tr_insertionsort(ISAd, first, last); - limit = -3; - continue; - } - - if(limit-- == 0) { - tr_heapsort(ISAd, first, last - first); - for(a = last - 1; first < a; a = b) { - for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; } - } - limit = -3; - continue; - } - - /* choose pivot */ - a = tr_pivot(ISAd, first, last); - SWAP(*first, *a); - v = ISAd[*first]; - - /* partition */ - tr_partition(ISAd, first, first + 1, last, &a, &b, v); - if((last - first) != (b - a)) { - next = (ISA[*a] != v) ? tr_ilg(b - a) : -1; - - /* update ranks */ - for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } - if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } } - - /* push */ - if((1 < (b - a)) && (trbudget_check(budget, b - a))) { - if((a - first) <= (last - b)) { - if((last - b) <= (b - a)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - STACK_PUSH5(ISAd, b, last, limit, trlink); - last = a; - } else if(1 < (last - b)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - first = b; - } else { - ISAd += incr, first = a, last = b, limit = next; - } - } else if((a - first) <= (b - a)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, limit, trlink); - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - last = a; - } else { - STACK_PUSH5(ISAd, b, last, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - STACK_PUSH5(ISAd, b, last, limit, trlink); - STACK_PUSH5(ISAd, first, a, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - if((a - first) <= (b - a)) { - if(1 < (last - b)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - STACK_PUSH5(ISAd, first, a, limit, trlink); - first = b; - } else if(1 < (a - first)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - last = a; - } else { - ISAd += incr, first = a, last = b, limit = next; - } - } else if((last - b) <= (b - a)) { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, limit, trlink); - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - first = b; - } else { - STACK_PUSH5(ISAd, first, a, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - STACK_PUSH5(ISAd, first, a, limit, trlink); - STACK_PUSH5(ISAd, b, last, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } - } else { - if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; } - if((a - first) <= (last - b)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, limit, trlink); - last = a; - } else if(1 < (last - b)) { - first = b; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } else { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, limit, trlink); - first = b; - } else if(1 < (a - first)) { - last = a; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } - } else { - if(trbudget_check(budget, last - first)) { - limit = tr_ilg(last - first), ISAd += incr; - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } -#undef STACK_SIZE -} - - - -/*---------------------------------------------------------------------------*/ - -/* Tandem repeat sort */ -static -void -trsort(int *ISA, int *SA, int n, int depth) { - int *ISAd; - int *first, *last; - trbudget_t budget; - int t, skip, unsorted; - - trbudget_init(&budget, tr_ilg(n) * 2 / 3, n); -/* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */ - for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) { - first = SA; - skip = 0; - unsorted = 0; - do { - if((t = *first) < 0) { first -= t; skip += t; } - else { - if(skip != 0) { *(first + skip) = skip; skip = 0; } - last = SA + ISA[t] + 1; - if(1 < (last - first)) { - budget.count = 0; - tr_introsort(ISA, ISAd, SA, first, last, &budget); - if(budget.count != 0) { unsorted += budget.count; } - else { skip = first - last; } - } else if((last - first) == 1) { - skip = -1; - } - first = last; - } - } while(first < (SA + n)); - if(skip != 0) { *(first + skip) = skip; } - if(unsorted == 0) { break; } - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Sorts suffixes of type B*. */ -static -int -sort_typeBstar(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int openMP) { - int *PAb, *ISAb, *buf; -#ifdef LIBBSC_OPENMP - int *curbuf; - int l; -#endif - int i, j, k, t, m, bufsize; - int c0, c1; -#ifdef LIBBSC_OPENMP - int d0, d1; -#endif - (void)openMP; - - /* Initialize bucket arrays. */ - for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; } - for(i = 0; i < BUCKET_B_SIZE; ++i) { bucket_B[i] = 0; } - - /* Count the number of occurrences of the first one or two characters of each - type A, B and B* suffix. Moreover, store the beginning position of all - type B* suffixes into the array SA. */ - for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) { - /* type A suffix. */ - do { ++BUCKET_A(c1 = c0); } while((0 <= --i) && ((c0 = T[i]) >= c1)); - if(0 <= i) { - /* type B* suffix. */ - ++BUCKET_BSTAR(c0, c1); - SA[--m] = i; - /* type B suffix. */ - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { - ++BUCKET_B(c0, c1); - } - } - } - m = n - m; -/* -note: - A type B* suffix is lexicographically smaller than a type B suffix that - begins with the same first two characters. -*/ - - /* Calculate the index of start/end point of each bucket. */ - for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0) { - t = i + BUCKET_A(c0); - BUCKET_A(c0) = i + j; /* start point */ - i = t + BUCKET_B(c0, c0); - for(c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1) { - j += BUCKET_BSTAR(c0, c1); - BUCKET_BSTAR(c0, c1) = j; /* end point */ - i += BUCKET_B(c0, c1); - } - } - - if(0 < m) { - /* Sort the type B* suffixes by their first two characters. */ - PAb = SA + n - m; ISAb = SA + m; - for(i = m - 2; 0 <= i; --i) { - t = PAb[i], c0 = T[t], c1 = T[t + 1]; - SA[--BUCKET_BSTAR(c0, c1)] = i; - } - t = PAb[m - 1], c0 = T[t], c1 = T[t + 1]; - SA[--BUCKET_BSTAR(c0, c1)] = m - 1; - - /* Sort the type B* substrings using sssort. */ -#ifdef LIBBSC_OPENMP - if (openMP) - { - buf = SA + m; - c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m; -#pragma omp parallel default(shared) private(bufsize, curbuf, k, l, d0, d1) - { - bufsize = (n - (2 * m)) / omp_get_num_threads(); - curbuf = buf + omp_get_thread_num() * bufsize; - k = 0; - for(;;) { - #pragma omp critical(sssort_lock) - { - if(0 < (l = j)) { - d0 = c0, d1 = c1; - do { - k = BUCKET_BSTAR(d0, d1); - if(--d1 <= d0) { - d1 = ALPHABET_SIZE - 1; - if(--d0 < 0) { break; } - } - } while(((l - k) <= 1) && (0 < (l = k))); - c0 = d0, c1 = d1, j = k; - } - } - if(l == 0) { break; } - sssort(T, PAb, SA + k, SA + l, - curbuf, bufsize, 2, n, *(SA + k) == (m - 1)); - } - } - } - else - { - buf = SA + m, bufsize = n - (2 * m); - for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) { - for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) { - i = BUCKET_BSTAR(c0, c1); - if(1 < (j - i)) { - sssort(T, PAb, SA + i, SA + j, - buf, bufsize, 2, n, *(SA + i) == (m - 1)); - } - } - } - } -#else - buf = SA + m, bufsize = n - (2 * m); - for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) { - for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) { - i = BUCKET_BSTAR(c0, c1); - if(1 < (j - i)) { - sssort(T, PAb, SA + i, SA + j, - buf, bufsize, 2, n, *(SA + i) == (m - 1)); - } - } - } -#endif - - /* Compute ranks of type B* substrings. */ - for(i = m - 1; 0 <= i; --i) { - if(0 <= SA[i]) { - j = i; - do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i])); - SA[i + 1] = i - j; - if(i <= 0) { break; } - } - j = i; - do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0); - ISAb[SA[i]] = j; - } - - /* Construct the inverse suffix array of type B* suffixes using trsort. */ - trsort(ISAb, SA, m, 1); - - /* Set the sorted order of type B* suffixes. */ - for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) { - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { } - if(0 <= i) { - t = i; - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { } - SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t; - } - } - - /* Calculate the index of start/end point of each bucket. */ - BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */ - for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) { - i = BUCKET_A(c0 + 1) - 1; - for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) { - t = i - BUCKET_B(c0, c1); - BUCKET_B(c0, c1) = i; /* end point */ - - /* Move all type B* suffixes to the correct position. */ - for(i = t, j = BUCKET_BSTAR(c0, c1); - j <= k; - --i, --k) { SA[i] = SA[k]; } - } - BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */ - BUCKET_B(c0, c0) = i; /* end point */ - } - } - - return m; -} - -/* Constructs the suffix array by using the sorted order of type B* suffixes. */ -static -void -construct_SA(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int m) { - int *i, *j, *k; - int s; - int c0, c1, c2; - - if(0 < m) { - /* Construct the sorted order of type B suffixes by using - the sorted order of type B* suffixes. */ - for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) { - /* Scan the suffix array from right to left. */ - for(i = SA + BUCKET_BSTAR(c1, c1 + 1), - j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; - i <= j; - --j) { - if(0 < (s = *j)) { - assert(T[s] == c1); - assert(((s + 1) < n) && (T[s] <= T[s + 1])); - assert(T[s - 1] <= T[s]); - *j = ~s; - c0 = T[--s]; - if((0 < s) && (T[s - 1] > c0)) { s = ~s; } - if(c0 != c2) { - if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; } - k = SA + BUCKET_B(c2 = c0, c1); - } - assert(k < j); assert(k != NULL); - *k-- = s; - } else { - assert(((s == 0) && (T[s] == c1)) || (s < 0)); - *j = ~s; - } - } - } - } - - /* Construct the suffix array by using - the sorted order of type B suffixes. */ - k = SA + BUCKET_A(c2 = T[n - 1]); - *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1); - /* Scan the suffix array from left to right. */ - for(i = SA, j = SA + n; i < j; ++i) { - if(0 < (s = *i)) { - assert(T[s - 1] >= T[s]); - c0 = T[--s]; - if((s == 0) || (T[s - 1] < c0)) { s = ~s; } - if(c0 != c2) { - BUCKET_A(c2) = k - SA; - k = SA + BUCKET_A(c2 = c0); - } - assert(i < k); - *k++ = s; - } else { - assert(s < 0); - *i = ~s; - } - } -} - -/* Constructs the burrows-wheeler transformed string directly - by using the sorted order of type B* suffixes. */ -static -int -construct_BWT(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int m) { - int *i, *j, *k, *orig; - int s; - int c0, c1, c2; - - if(0 < m) { - /* Construct the sorted order of type B suffixes by using - the sorted order of type B* suffixes. */ - for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) { - /* Scan the suffix array from right to left. */ - for(i = SA + BUCKET_BSTAR(c1, c1 + 1), - j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; - i <= j; - --j) { - if(0 < (s = *j)) { - assert(T[s] == c1); - assert(((s + 1) < n) && (T[s] <= T[s + 1])); - assert(T[s - 1] <= T[s]); - c0 = T[--s]; - *j = ~((int)c0); - if((0 < s) && (T[s - 1] > c0)) { s = ~s; } - if(c0 != c2) { - if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; } - k = SA + BUCKET_B(c2 = c0, c1); - } - assert(k < j); assert(k != NULL); - *k-- = s; - } else if(s != 0) { - *j = ~s; -#ifndef NDEBUG - } else { - assert(T[s] == c1); -#endif - } - } - } - } - - /* Construct the BWTed string by using - the sorted order of type B suffixes. */ - k = SA + BUCKET_A(c2 = T[n - 1]); - *k++ = (T[n - 2] < c2) ? ~((int)T[n - 2]) : (n - 1); - /* Scan the suffix array from left to right. */ - for(i = SA, j = SA + n, orig = SA; i < j; ++i) { - if(0 < (s = *i)) { - assert(T[s - 1] >= T[s]); - c0 = T[--s]; - *i = c0; - if((0 < s) && (T[s - 1] < c0)) { s = ~((int)T[s - 1]); } - if(c0 != c2) { - BUCKET_A(c2) = k - SA; - k = SA + BUCKET_A(c2 = c0); - } - assert(i < k); - *k++ = s; - } else if(s != 0) { - *i = ~s; - } else { - orig = i; - } - } - - return orig - SA; -} - -/* Constructs the burrows-wheeler transformed string directly - by using the sorted order of type B* suffixes. */ -static -int -construct_BWT_indexes(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int m, - unsigned char * num_indexes, int * indexes) { - int *i, *j, *k, *orig; - int s; - int c0, c1, c2; - - int mod = n / 8; - { - mod |= mod >> 1; mod |= mod >> 2; - mod |= mod >> 4; mod |= mod >> 8; - mod |= mod >> 16; mod >>= 1; - - *num_indexes = (unsigned char)((n - 1) / (mod + 1)); - } - - if(0 < m) { - /* Construct the sorted order of type B suffixes by using - the sorted order of type B* suffixes. */ - for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) { - /* Scan the suffix array from right to left. */ - for(i = SA + BUCKET_BSTAR(c1, c1 + 1), - j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; - i <= j; - --j) { - if(0 < (s = *j)) { - assert(T[s] == c1); - assert(((s + 1) < n) && (T[s] <= T[s + 1])); - assert(T[s - 1] <= T[s]); - - if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = j - SA; - - c0 = T[--s]; - *j = ~((int)c0); - if((0 < s) && (T[s - 1] > c0)) { s = ~s; } - if(c0 != c2) { - if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; } - k = SA + BUCKET_B(c2 = c0, c1); - } - assert(k < j); assert(k != NULL); - *k-- = s; - } else if(s != 0) { - *j = ~s; -#ifndef NDEBUG - } else { - assert(T[s] == c1); -#endif - } - } - } - } - - /* Construct the BWTed string by using - the sorted order of type B suffixes. */ - k = SA + BUCKET_A(c2 = T[n - 1]); - if (T[n - 2] < c2) { - if (((n - 1) & mod) == 0) indexes[(n - 1) / (mod + 1) - 1] = k - SA; - *k++ = ~((int)T[n - 2]); - } - else { - *k++ = n - 1; - } - - /* Scan the suffix array from left to right. */ - for(i = SA, j = SA + n, orig = SA; i < j; ++i) { - if(0 < (s = *i)) { - assert(T[s - 1] >= T[s]); - - if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = i - SA; - - c0 = T[--s]; - *i = c0; - if(c0 != c2) { - BUCKET_A(c2) = k - SA; - k = SA + BUCKET_A(c2 = c0); - } - assert(i < k); - if((0 < s) && (T[s - 1] < c0)) { - if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = k - SA; - *k++ = ~((int)T[s - 1]); - } else - *k++ = s; - } else if(s != 0) { - *i = ~s; - } else { - orig = i; - } - } - - return orig - SA; -} - - -/*---------------------------------------------------------------------------*/ - -/*- Function -*/ - -int -divsufsort(const unsigned char *T, int *SA, int n, int openMP) { - int *bucket_A, *bucket_B; - int m; - int err = 0; - - /* Check arguments. */ - if((T == NULL) || (SA == NULL) || (n < 0)) { return -1; } - else if(n == 0) { return 0; } - else if(n == 1) { SA[0] = 0; return 0; } - else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; } - - bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int)); - bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int)); - - /* Suffixsort. */ - if((bucket_A != NULL) && (bucket_B != NULL)) { - m = sort_typeBstar(T, SA, bucket_A, bucket_B, n, openMP); - construct_SA(T, SA, bucket_A, bucket_B, n, m); - } else { - err = -2; - } - - free(bucket_B); - free(bucket_A); - - return err; -} - -int -divbwt(const unsigned char *T, unsigned char *U, int *A, int n, unsigned char * num_indexes, int * indexes, int openMP) { - int *B; - int *bucket_A, *bucket_B; - int m, pidx, i; - - /* Check arguments. */ - if((T == NULL) || (U == NULL) || (n < 0)) { return -1; } - else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; } - - if((B = A) == NULL) { B = (int *)malloc((size_t)(n + 1) * sizeof(int)); } - bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int)); - bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int)); - - /* Burrows-Wheeler Transform. */ - if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) { - m = sort_typeBstar(T, B, bucket_A, bucket_B, n, openMP); - - if (num_indexes == NULL || indexes == NULL) { - pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m); - } else { - pidx = construct_BWT_indexes(T, B, bucket_A, bucket_B, n, m, num_indexes, indexes); - } - - /* Copy to output string. */ - U[0] = T[n - 1]; - for(i = 0; i < pidx; ++i) { U[i + 1] = (unsigned char)B[i]; } - for(i += 1; i < n; ++i) { U[i] = (unsigned char)B[i]; } - pidx += 1; - } else { - pidx = -2; - } - - free(bucket_B); - free(bucket_A); - if(A == NULL) { free(B); } - - return pidx; -} diff --git a/rust/README.md b/rust/README.md index cfab35a4e..d3ff3a5b2 100644 --- a/rust/README.md +++ b/rust/README.md @@ -42,6 +42,10 @@ zstd ABI: the dynamic-programming optimal parser itself remains in C for now. - `zstd_ldm` implements long-distance-match parameter selection, table maintenance, sequence generation, and sequence consumption. +- Dictionary building + - `divsufsort` constructs the suffix array that drives the legacy `ZDICT` + trainer (`ZDICT_trainFromBuffer_legacy`). The sample analysis and + dictionary assembly in `zdict.c`, `cover.c`, and `fastcover.c` remain C. - Runtime support - `threading` provides platform pthread wrappers required by zstd headers. - `pool` implements the bounded worker pool used by multithreaded compression. @@ -60,10 +64,11 @@ zstd ABI: `fileio` backend still owns file opening, safe replacement, sparse writes, metadata, and streaming I/O. -The optimal block matcher, high-level frame compression, dictionary-building, -legacy decoding callbacks, and the CLI file-I/O backend are still C. They must -move before the rewrite is complete. Keeping that boundary explicit prevents a -passing hybrid build from being mistaken for the final all-Rust result. +The optimal block matcher, high-level frame compression, dictionary-building +except suffix-array construction, legacy decoding callbacks, and the CLI +file-I/O backend are still C. They must move before the rewrite is complete. +Keeping that boundary explicit prevents a passing hybrid build from being +mistaken for the final all-Rust result. ## Compatibility boundary diff --git a/rust/src/divsufsort.rs b/rust/src/divsufsort.rs new file mode 100644 index 000000000..a9b54a3e1 --- /dev/null +++ b/rust/src/divsufsort.rs @@ -0,0 +1,2756 @@ +#![allow(clippy::missing_safety_doc)] +#![allow(clippy::too_many_arguments)] + +//! Suffix-array construction for the dictionary builder. +//! +//! Port of `lib/dictBuilder/divsufsort.c` (libdivsufsort-lite, Copyright (c) +//! 2003-2008 Yuta Mori, MIT license) in the exact configuration zstd compiles +//! it with: `ALPHABET_SIZE = 256`, `SS_INSERTIONSORT_THRESHOLD = 8`, +//! `SS_BLOCKSIZE = 1024`, and no OpenMP. Only `divsufsort()` is exported; +//! `divbwt()` has no callers anywhere in zstd and was not ported. +//! +//! The C implementation walks raw `int*` cursors through the caller's SA +//! buffer, including transient one-before-the-range positions. Every such +//! cursor is translated to an `isize` index into one `&mut [i32]` slice +//! covering the whole buffer, so all arithmetic — including the +//! bitwise-complement rank marking and the C `int` value semantics — matches +//! the original exactly while staying bounds-checked. + +use std::os::raw::c_int; +use std::slice; + +const BUCKET_A_SIZE: usize = 256; /* ALPHABET_SIZE */ +const BUCKET_B_SIZE: usize = 256 * 256; /* ALPHABET_SIZE * ALPHABET_SIZE */ +const ALPHABET_SIZE: i32 = 256; +const SS_INSERTIONSORT_THRESHOLD: isize = 8; +const SS_BLOCKSIZE: isize = 1024; +/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */ +const SS_MISORT_STACKSIZE: usize = 16; +const SS_SMERGE_STACKSIZE: usize = 32; +const TR_INSERTIONSORT_THRESHOLD: isize = 8; +const TR_STACKSIZE: usize = 64; + +#[rustfmt::skip] +static LG_TABLE: [i32; 256] = [ + -1,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, + 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, + 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, + 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, +]; + +#[rustfmt::skip] +static SQQ_TABLE: [i32; 256] = [ + 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, + 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, + 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109, +110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, +128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, +143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, +156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, +169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, +181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, +192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, +202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211, +212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, +221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, +230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, +239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, +247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255, +]; + +/* `ss_ilg` in its `256 <= SS_BLOCKSIZE` configuration. */ +#[inline] +fn ss_ilg(n: isize) -> i32 { + let n = n as i32; + if n & 0xff00 != 0 { + 8 + LG_TABLE[((n >> 8) & 0xff) as usize] + } else { + LG_TABLE[(n & 0xff) as usize] + } +} + +#[inline] +fn ss_isqrt(x: isize) -> isize { + if x >= SS_BLOCKSIZE * SS_BLOCKSIZE { + return SS_BLOCKSIZE; + } + let x = x as i32; + let e = if (x as u32) & 0xffff_0000 != 0 { + if (x as u32) & 0xff00_0000 != 0 { + 24 + LG_TABLE[((x >> 24) & 0xff) as usize] + } else { + 16 + LG_TABLE[((x >> 16) & 0xff) as usize] + } + } else if x & 0xff00 != 0 { + 8 + LG_TABLE[((x >> 8) & 0xff) as usize] + } else { + LG_TABLE[(x & 0xff) as usize] + }; + + let mut y; + if e >= 16 { + y = SQQ_TABLE[(x >> ((e - 6) - (e & 1))) as usize] << ((e >> 1) - 7); + if e >= 24 { + y = (y + 1 + x / y) >> 1; + } + y = (y + 1 + x / y) >> 1; + } else if e >= 8 { + y = (SQQ_TABLE[(x >> ((e - 6) - (e & 1))) as usize] >> (7 - (e >> 1))) + 1; + } else { + return (SQQ_TABLE[x as usize] >> 4) as isize; + } + + (if x < y * y { y - 1 } else { y }) as isize +} + +/* --------------------------------------------------------------------- */ + +/// Compares two suffixes. `(p10, p11)` and `(p20, p21)` are the `p[0]`/`p[1]` +/// pairs the C routine reads through its `const int*` arguments; passing the +/// values directly also serves `sssort()`'s local two-element `PAi` array. +#[inline] +fn ss_compare(t: &[u8], p10: i32, p11: i32, p20: i32, p21: i32, depth: i32) -> i32 { + let mut u1 = (depth + p10) as isize; + let mut u2 = (depth + p20) as isize; + let u1n = (p11 + 2) as isize; + let u2n = (p21 + 2) as isize; + + while u1 < u1n && u2 < u2n && t[u1 as usize] == t[u2 as usize] { + u1 += 1; + u2 += 1; + } + + if u1 < u1n { + if u2 < u2n { + t[u1 as usize] as i32 - t[u2 as usize] as i32 + } else { + 1 + } + } else if u2 < u2n { + -1 + } else { + 0 + } +} + +/// `ss_compare(T, p1, p2, depth)` for pointers `p1`/`p2` into the SA buffer. +#[inline] +fn ss_compare_pa(t: &[u8], sa: &[i32], p1: isize, p2: isize, depth: i32) -> i32 { + ss_compare( + t, + sa[p1 as usize], + sa[(p1 + 1) as usize], + sa[p2 as usize], + sa[(p2 + 1) as usize], + depth, + ) +} + +/* --------------------------------------------------------------------- */ + +/* Insertionsort for small size groups */ +fn ss_insertionsort(t: &[u8], sa: &mut [i32], pa: isize, first: isize, last: isize, depth: i32) { + let mut i = last - 2; + while first <= i { + let t0 = sa[i as usize]; + let mut j = i + 1; + let mut r; + loop { + r = ss_compare_pa(t, sa, pa + t0 as isize, pa + sa[j as usize] as isize, depth); + if r <= 0 { + break; + } + loop { + sa[(j - 1) as usize] = sa[j as usize]; + j += 1; + if !(j < last && sa[j as usize] < 0) { + break; + } + } + if last <= j { + break; + } + } + if r == 0 { + sa[j as usize] = !sa[j as usize]; + } + sa[(j - 1) as usize] = t0; + i -= 1; + } +} + +/* --------------------------------------------------------------------- */ + +/// `Td[PA[SA[p]]]` — the depth-`td` sorting key of the suffix stored at `p`. +#[inline(always)] +fn ss_key(t: &[u8], sa: &[i32], td: isize, pa: isize, p: isize) -> i32 { + t[(td + sa[(pa + sa[p as usize] as isize) as usize] as isize) as usize] as i32 +} + +/// `Td[v]` for an already-loaded SA element `v` (`Td[PA[v]]` in C). +#[inline(always)] +fn ss_key_of(t: &[u8], sa: &[i32], td: isize, pa: isize, v: i32) -> i32 { + t[(td + sa[(pa + v as isize) as usize] as isize) as usize] as i32 +} + +/// `Td[PA[SA[p]] - 1]` — the character preceding the depth-`td` key. +#[inline(always)] +fn ss_key_pred(t: &[u8], sa: &[i32], td: isize, pa: isize, p: isize) -> i32 { + t[(td + sa[(pa + sa[p as usize] as isize) as usize] as isize - 1) as usize] as i32 +} + +fn ss_fixdown(t: &[u8], td: isize, sa: &mut [i32], pa: isize, base: isize, i: isize, size: isize) { + let mut i = i; + let v = sa[(base + i) as usize]; + let c = ss_key_of(t, sa, td, pa, v); + loop { + let mut j = 2 * i + 1; + if j >= size { + break; + } + let mut k = j; + j += 1; + let mut d = ss_key(t, sa, td, pa, base + k); + let e = ss_key(t, sa, td, pa, base + j); + if d < e { + k = j; + d = e; + } + if d <= c { + break; + } + sa[(base + i) as usize] = sa[(base + k) as usize]; + i = k; + } + sa[(base + i) as usize] = v; +} + +/* Simple top-down heapsort. */ +fn ss_heapsort(t: &[u8], td: isize, sa: &mut [i32], pa: isize, base: isize, size: isize) { + let mut m = size; + if size % 2 == 0 { + m -= 1; + if ss_key(t, sa, td, pa, base + m / 2) < ss_key(t, sa, td, pa, base + m) { + sa.swap((base + m) as usize, (base + m / 2) as usize); + } + } + + let mut i = m / 2 - 1; + while 0 <= i { + ss_fixdown(t, td, sa, pa, base, i, m); + i -= 1; + } + if size % 2 == 0 { + sa.swap(base as usize, (base + m) as usize); + ss_fixdown(t, td, sa, pa, base, 0, m); + } + let mut i = m - 1; + while 0 < i { + let t0 = sa[base as usize]; + sa[base as usize] = sa[(base + i) as usize]; + ss_fixdown(t, td, sa, pa, base, 0, i); + sa[(base + i) as usize] = t0; + i -= 1; + } +} + +/* --------------------------------------------------------------------- */ + +/* Returns the median of three elements. */ +#[inline] +fn ss_median3( + t: &[u8], + sa: &[i32], + td: isize, + pa: isize, + v1: isize, + v2: isize, + v3: isize, +) -> isize { + let mut v1 = v1; + let mut v2 = v2; + if ss_key(t, sa, td, pa, v1) > ss_key(t, sa, td, pa, v2) { + std::mem::swap(&mut v1, &mut v2); + } + if ss_key(t, sa, td, pa, v2) > ss_key(t, sa, td, pa, v3) { + if ss_key(t, sa, td, pa, v1) > ss_key(t, sa, td, pa, v3) { + return v1; + } + return v3; + } + v2 +} + +/* Returns the median of five elements. */ +#[inline] +fn ss_median5( + t: &[u8], + sa: &[i32], + td: isize, + pa: isize, + v1: isize, + v2: isize, + v3: isize, + v4: isize, + v5: isize, +) -> isize { + let mut v1 = v1; + let mut v2 = v2; + let mut v3 = v3; + let mut v4 = v4; + let mut v5 = v5; + if ss_key(t, sa, td, pa, v2) > ss_key(t, sa, td, pa, v3) { + std::mem::swap(&mut v2, &mut v3); + } + if ss_key(t, sa, td, pa, v4) > ss_key(t, sa, td, pa, v5) { + std::mem::swap(&mut v4, &mut v5); + } + if ss_key(t, sa, td, pa, v2) > ss_key(t, sa, td, pa, v4) { + std::mem::swap(&mut v2, &mut v4); + std::mem::swap(&mut v3, &mut v5); + } + if ss_key(t, sa, td, pa, v1) > ss_key(t, sa, td, pa, v3) { + std::mem::swap(&mut v1, &mut v3); + } + if ss_key(t, sa, td, pa, v1) > ss_key(t, sa, td, pa, v4) { + std::mem::swap(&mut v1, &mut v4); + std::mem::swap(&mut v3, &mut v5); + } + if ss_key(t, sa, td, pa, v3) > ss_key(t, sa, td, pa, v4) { + return v4; + } + v3 +} + +/* Returns the pivot element. */ +#[inline] +fn ss_pivot(t: &[u8], sa: &[i32], td: isize, pa: isize, first: isize, last: isize) -> isize { + let mut t0 = last - first; + let middle = first + t0 / 2; + + if t0 <= 512 { + if t0 <= 32 { + return ss_median3(t, sa, td, pa, first, middle, last - 1); + } + t0 >>= 2; + return ss_median5( + t, + sa, + td, + pa, + first, + first + t0, + middle, + last - 1 - t0, + last - 1, + ); + } + t0 >>= 3; + let first = ss_median3(t, sa, td, pa, first, first + t0, first + (t0 << 1)); + let middle = ss_median3(t, sa, td, pa, middle - t0, middle, middle + t0); + let last = ss_median3(t, sa, td, pa, last - 1 - (t0 << 1), last - 1 - t0, last - 1); + ss_median3(t, sa, td, pa, first, middle, last) +} + +/* --------------------------------------------------------------------- */ + +/* Binary partition for substrings. */ +/* The `>= x + 1` comparison deliberately mirrors the C expression shape. */ +#[allow(clippy::int_plus_one)] +fn ss_partition(sa: &mut [i32], pa: isize, first: isize, last: isize, depth: i32) -> isize { + let mut a = first - 1; + let mut b = last; + loop { + loop { + a += 1; + if !(a < b) { + break; + } + if !(sa[(pa + sa[a as usize] as isize) as usize] + depth + >= sa[(pa + sa[a as usize] as isize + 1) as usize] + 1) + { + break; + } + sa[a as usize] = !sa[a as usize]; + } + loop { + b -= 1; + if !(a < b) { + break; + } + if !(sa[(pa + sa[b as usize] as isize) as usize] + depth + < sa[(pa + sa[b as usize] as isize + 1) as usize] + 1) + { + break; + } + } + if b <= a { + break; + } + let t0 = !sa[b as usize]; + sa[b as usize] = sa[a as usize]; + sa[a as usize] = t0; + } + if first < a { + sa[first as usize] = !sa[first as usize]; + } + a +} + +/* Multikey introsort for medium size groups. */ +fn ss_mintrosort(t: &[u8], sa: &mut [i32], pa: isize, first: isize, last: isize, depth: i32) { + let mut stack = [(0isize, 0isize, 0i32, 0i32); SS_MISORT_STACKSIZE]; + let mut ssize = 0usize; + let mut first = first; + let mut last = last; + let mut depth = depth; + let mut limit = ss_ilg(last - first); + let mut x: i32 = 0; + + loop { + if last - first <= SS_INSERTIONSORT_THRESHOLD { + if 1 < last - first { + ss_insertionsort(t, sa, pa, first, last, depth); + } + /* STACK_POP */ + if ssize == 0 { + return; + } + ssize -= 1; + (first, last, depth, limit) = stack[ssize]; + continue; + } + + let td = depth as isize; + if limit == 0 { + ss_heapsort(t, td, sa, pa, first, last - first); + } + limit -= 1; + if limit < 0 { + let mut a = first + 1; + let mut v = ss_key(t, sa, td, pa, first); + while a < last { + x = ss_key(t, sa, td, pa, a); + if x != v { + if 1 < a - first { + break; + } + v = x; + first = a; + } + a += 1; + } + if ss_key_pred(t, sa, td, pa, first) < v { + first = ss_partition(sa, pa, first, a, depth); + } + if a - first <= last - a { + if 1 < a - first { + stack[ssize] = (a, last, depth, -1); + ssize += 1; + last = a; + depth += 1; + limit = ss_ilg(a - first); + } else { + first = a; + limit = -1; + } + } else if 1 < last - a { + stack[ssize] = (first, a, depth + 1, ss_ilg(a - first)); + ssize += 1; + first = a; + limit = -1; + } else { + last = a; + depth += 1; + limit = ss_ilg(a - first); + } + continue; + } + + /* choose pivot */ + let mut a = ss_pivot(t, sa, td, pa, first, last); + let v = ss_key(t, sa, td, pa, a); + sa.swap(first as usize, a as usize); + + /* partition */ + let mut b = first; + loop { + b += 1; + if !(b < last) { + break; + } + x = ss_key(t, sa, td, pa, b); + if x != v { + break; + } + } + a = b; + if a < last && x < v { + loop { + b += 1; + if !(b < last) { + break; + } + x = ss_key(t, sa, td, pa, b); + if !(x <= v) { + break; + } + if x == v { + sa.swap(b as usize, a as usize); + a += 1; + } + } + } + let mut c = last; + loop { + c -= 1; + if !(b < c) { + break; + } + x = ss_key(t, sa, td, pa, c); + if x != v { + break; + } + } + let mut d = c; + if b < d && x > v { + loop { + c -= 1; + if !(b < c) { + break; + } + x = ss_key(t, sa, td, pa, c); + if !(x >= v) { + break; + } + if x == v { + sa.swap(c as usize, d as usize); + d -= 1; + } + } + } + while b < c { + sa.swap(b as usize, c as usize); + loop { + b += 1; + if !(b < c) { + break; + } + x = ss_key(t, sa, td, pa, b); + if !(x <= v) { + break; + } + if x == v { + sa.swap(b as usize, a as usize); + a += 1; + } + } + loop { + c -= 1; + if !(b < c) { + break; + } + x = ss_key(t, sa, td, pa, c); + if !(x >= v) { + break; + } + if x == v { + sa.swap(c as usize, d as usize); + d -= 1; + } + } + } + + if a <= d { + c = b - 1; + + let mut s = a - first; + let t0 = b - a; + if s > t0 { + s = t0; + } + let mut e = first; + let mut f = b - s; + while 0 < s { + sa.swap(e as usize, f as usize); + s -= 1; + e += 1; + f += 1; + } + let mut s = d - c; + let t0 = last - d - 1; + if s > t0 { + s = t0; + } + let mut e = b; + let mut f = last - s; + while 0 < s { + sa.swap(e as usize, f as usize); + s -= 1; + e += 1; + f += 1; + } + + a = first + (b - a); + c = last - (d - c); + b = if v <= ss_key_pred(t, sa, td, pa, a) { + a + } else { + ss_partition(sa, pa, a, c, depth) + }; + + if a - first <= last - c { + if last - c <= c - b { + stack[ssize] = (b, c, depth + 1, ss_ilg(c - b)); + ssize += 1; + stack[ssize] = (c, last, depth, limit); + ssize += 1; + last = a; + } else if a - first <= c - b { + stack[ssize] = (c, last, depth, limit); + ssize += 1; + stack[ssize] = (b, c, depth + 1, ss_ilg(c - b)); + ssize += 1; + last = a; + } else { + stack[ssize] = (c, last, depth, limit); + ssize += 1; + stack[ssize] = (first, a, depth, limit); + ssize += 1; + first = b; + last = c; + depth += 1; + limit = ss_ilg(c - b); + } + } else if a - first <= c - b { + stack[ssize] = (b, c, depth + 1, ss_ilg(c - b)); + ssize += 1; + stack[ssize] = (first, a, depth, limit); + ssize += 1; + first = c; + } else if last - c <= c - b { + stack[ssize] = (first, a, depth, limit); + ssize += 1; + stack[ssize] = (b, c, depth + 1, ss_ilg(c - b)); + ssize += 1; + first = c; + } else { + stack[ssize] = (first, a, depth, limit); + ssize += 1; + stack[ssize] = (c, last, depth, limit); + ssize += 1; + first = b; + last = c; + depth += 1; + limit = ss_ilg(c - b); + } + } else { + limit += 1; + if ss_key_pred(t, sa, td, pa, first) < v { + first = ss_partition(sa, pa, first, last, depth); + limit = ss_ilg(last - first); + } + depth += 1; + } + } +} + +/* --------------------------------------------------------------------- */ + +#[inline] +fn ss_blockswap(sa: &mut [i32], a: isize, b: isize, n: isize) { + let mut a = a; + let mut b = b; + let mut n = n; + while 0 < n { + sa.swap(a as usize, b as usize); + n -= 1; + a += 1; + b += 1; + } +} + +#[inline] +fn ss_rotate(sa: &mut [i32], first: isize, middle: isize, last: isize) { + let mut first = first; + let mut last = last; + let mut l = middle - first; + let mut r = last - middle; + while 0 < l && 0 < r { + if l == r { + ss_blockswap(sa, first, middle, l); + break; + } + if l < r { + let mut a = last - 1; + let mut b = middle - 1; + let mut t0 = sa[a as usize]; + loop { + sa[a as usize] = sa[b as usize]; + a -= 1; + sa[b as usize] = sa[a as usize]; + b -= 1; + if b < first { + sa[a as usize] = t0; + last = a; + r -= l + 1; + if r <= l { + break; + } + a -= 1; + b = middle - 1; + t0 = sa[a as usize]; + } + } + } else { + let mut a = first; + let mut b = middle; + let mut t0 = sa[a as usize]; + loop { + sa[a as usize] = sa[b as usize]; + a += 1; + sa[b as usize] = sa[a as usize]; + b += 1; + if last <= b { + sa[a as usize] = t0; + first = a + 1; + l -= r + 1; + if l <= r { + break; + } + a += 1; + b = middle; + t0 = sa[a as usize]; + } + } + } + } +} + +/* --------------------------------------------------------------------- */ + +fn ss_inplacemerge( + t: &[u8], + sa: &mut [i32], + pa: isize, + first: isize, + middle: isize, + last: isize, + depth: i32, +) { + let mut middle = middle; + let mut last = last; + loop { + let x: i32; + let p: isize; + if sa[(last - 1) as usize] < 0 { + x = 1; + p = pa + (!sa[(last - 1) as usize]) as isize; + } else { + x = 0; + p = pa + sa[(last - 1) as usize] as isize; + } + let mut a = first; + let mut len = middle - first; + let mut half = len >> 1; + let mut r: i32 = -1; + while 0 < len { + let b = a + half; + let bv = sa[b as usize]; + let q = ss_compare_pa( + t, + sa, + pa + (if 0 <= bv { bv } else { !bv }) as isize, + p, + depth, + ); + if q < 0 { + a = b + 1; + half -= (len & 1) ^ 1; + } else { + r = q; + } + len = half; + half >>= 1; + } + if a < middle { + if r == 0 { + sa[a as usize] = !sa[a as usize]; + } + ss_rotate(sa, a, middle, last); + last -= middle - a; + middle = a; + if first == middle { + break; + } + } + last -= 1; + if x != 0 { + loop { + last -= 1; + if !(sa[last as usize] < 0) { + break; + } + } + } + if middle == last { + break; + } + } +} + +/* --------------------------------------------------------------------- */ + +/* Merge-forward with internal buffer. */ +fn ss_mergeforward( + t: &[u8], + sa: &mut [i32], + pa: isize, + first: isize, + middle: isize, + last: isize, + buf: isize, + depth: i32, +) { + let bufend = buf + (middle - first) - 1; + ss_blockswap(sa, buf, first, middle - first); + + let mut a = first; + let t0 = sa[a as usize]; + let mut b = buf; + let mut c = middle; + loop { + let r = ss_compare_pa( + t, + sa, + pa + sa[b as usize] as isize, + pa + sa[c as usize] as isize, + depth, + ); + if r < 0 { + loop { + sa[a as usize] = sa[b as usize]; + a += 1; + if bufend <= b { + sa[bufend as usize] = t0; + return; + } + sa[b as usize] = sa[a as usize]; + b += 1; + if !(sa[b as usize] < 0) { + break; + } + } + } else if r > 0 { + loop { + sa[a as usize] = sa[c as usize]; + a += 1; + sa[c as usize] = sa[a as usize]; + c += 1; + if last <= c { + while b < bufend { + sa[a as usize] = sa[b as usize]; + a += 1; + sa[b as usize] = sa[a as usize]; + b += 1; + } + sa[a as usize] = sa[b as usize]; + sa[b as usize] = t0; + return; + } + if !(sa[c as usize] < 0) { + break; + } + } + } else { + sa[c as usize] = !sa[c as usize]; + loop { + sa[a as usize] = sa[b as usize]; + a += 1; + if bufend <= b { + sa[bufend as usize] = t0; + return; + } + sa[b as usize] = sa[a as usize]; + b += 1; + if !(sa[b as usize] < 0) { + break; + } + } + loop { + sa[a as usize] = sa[c as usize]; + a += 1; + sa[c as usize] = sa[a as usize]; + c += 1; + if last <= c { + while b < bufend { + sa[a as usize] = sa[b as usize]; + a += 1; + sa[b as usize] = sa[a as usize]; + b += 1; + } + sa[a as usize] = sa[b as usize]; + sa[b as usize] = t0; + return; + } + if !(sa[c as usize] < 0) { + break; + } + } + } + } +} + +/* Merge-backward with internal buffer. */ +fn ss_mergebackward( + t: &[u8], + sa: &mut [i32], + pa: isize, + first: isize, + middle: isize, + last: isize, + buf: isize, + depth: i32, +) { + let bufend = buf + (last - middle) - 1; + ss_blockswap(sa, buf, middle, last - middle); + + let mut x = 0i32; + let mut p1: isize; + let mut p2: isize; + if sa[bufend as usize] < 0 { + p1 = pa + (!sa[bufend as usize]) as isize; + x |= 1; + } else { + p1 = pa + sa[bufend as usize] as isize; + } + if sa[(middle - 1) as usize] < 0 { + p2 = pa + (!sa[(middle - 1) as usize]) as isize; + x |= 2; + } else { + p2 = pa + sa[(middle - 1) as usize] as isize; + } + let mut a = last - 1; + let t0 = sa[a as usize]; + let mut b = bufend; + let mut c = middle - 1; + loop { + let r = ss_compare_pa(t, sa, p1, p2, depth); + if 0 < r { + if x & 1 != 0 { + loop { + sa[a as usize] = sa[b as usize]; + a -= 1; + sa[b as usize] = sa[a as usize]; + b -= 1; + if !(sa[b as usize] < 0) { + break; + } + } + x ^= 1; + } + sa[a as usize] = sa[b as usize]; + a -= 1; + if b <= buf { + sa[buf as usize] = t0; + break; + } + sa[b as usize] = sa[a as usize]; + b -= 1; + if sa[b as usize] < 0 { + p1 = pa + (!sa[b as usize]) as isize; + x |= 1; + } else { + p1 = pa + sa[b as usize] as isize; + } + } else if r < 0 { + if x & 2 != 0 { + loop { + sa[a as usize] = sa[c as usize]; + a -= 1; + sa[c as usize] = sa[a as usize]; + c -= 1; + if !(sa[c as usize] < 0) { + break; + } + } + x ^= 2; + } + sa[a as usize] = sa[c as usize]; + a -= 1; + sa[c as usize] = sa[a as usize]; + c -= 1; + if c < first { + while buf < b { + sa[a as usize] = sa[b as usize]; + a -= 1; + sa[b as usize] = sa[a as usize]; + b -= 1; + } + sa[a as usize] = sa[b as usize]; + sa[b as usize] = t0; + break; + } + if sa[c as usize] < 0 { + p2 = pa + (!sa[c as usize]) as isize; + x |= 2; + } else { + p2 = pa + sa[c as usize] as isize; + } + } else { + if x & 1 != 0 { + loop { + sa[a as usize] = sa[b as usize]; + a -= 1; + sa[b as usize] = sa[a as usize]; + b -= 1; + if !(sa[b as usize] < 0) { + break; + } + } + x ^= 1; + } + sa[a as usize] = !sa[b as usize]; + a -= 1; + if b <= buf { + sa[buf as usize] = t0; + break; + } + sa[b as usize] = sa[a as usize]; + b -= 1; + if x & 2 != 0 { + loop { + sa[a as usize] = sa[c as usize]; + a -= 1; + sa[c as usize] = sa[a as usize]; + c -= 1; + if !(sa[c as usize] < 0) { + break; + } + } + x ^= 2; + } + sa[a as usize] = sa[c as usize]; + a -= 1; + sa[c as usize] = sa[a as usize]; + c -= 1; + if c < first { + while buf < b { + sa[a as usize] = sa[b as usize]; + a -= 1; + sa[b as usize] = sa[a as usize]; + b -= 1; + } + sa[a as usize] = sa[b as usize]; + sa[b as usize] = t0; + break; + } + if sa[b as usize] < 0 { + p1 = pa + (!sa[b as usize]) as isize; + x |= 1; + } else { + p1 = pa + sa[b as usize] as isize; + } + if sa[c as usize] < 0 { + p2 = pa + (!sa[c as usize]) as isize; + x |= 2; + } else { + p2 = pa + sa[c as usize] as isize; + } + } + } +} + +/// `GETIDX` — undoes the "already merged" complement marking. +#[inline(always)] +fn getidx(a: i32) -> i32 { + if 0 <= a { + a + } else { + !a + } +} + +/// `MERGE_CHECK` — restores or sets the complement marks after a merge. +#[inline] +fn ss_merge_check(t: &[u8], sa: &mut [i32], pa: isize, a: isize, b: isize, c: i32, depth: i32) { + if (c & 1) != 0 + || ((c & 2) != 0 + && ss_compare_pa( + t, + sa, + pa + getidx(sa[(a - 1) as usize]) as isize, + pa + sa[a as usize] as isize, + depth, + ) == 0) + { + sa[a as usize] = !sa[a as usize]; + } + if (c & 4) != 0 + && ss_compare_pa( + t, + sa, + pa + getidx(sa[(b - 1) as usize]) as isize, + pa + sa[b as usize] as isize, + depth, + ) == 0 + { + sa[b as usize] = !sa[b as usize]; + } +} + +/* D&C based merge. */ +fn ss_swapmerge( + t: &[u8], + sa: &mut [i32], + pa: isize, + first: isize, + middle: isize, + last: isize, + buf: isize, + bufsize: isize, + depth: i32, +) { + let mut stack = [(0isize, 0isize, 0isize, 0i32); SS_SMERGE_STACKSIZE]; + let mut ssize = 0usize; + let mut first = first; + let mut middle = middle; + let mut last = last; + let mut check = 0i32; + + loop { + if last - middle <= bufsize { + if first < middle && middle < last { + ss_mergebackward(t, sa, pa, first, middle, last, buf, depth); + } + ss_merge_check(t, sa, pa, first, last, check, depth); + if ssize == 0 { + return; + } + ssize -= 1; + (first, middle, last, check) = stack[ssize]; + continue; + } + + if middle - first <= bufsize { + if first < middle { + ss_mergeforward(t, sa, pa, first, middle, last, buf, depth); + } + ss_merge_check(t, sa, pa, first, last, check, depth); + if ssize == 0 { + return; + } + ssize -= 1; + (first, middle, last, check) = stack[ssize]; + continue; + } + + let mut m: isize = 0; + let mut len = std::cmp::min(middle - first, last - middle); + let mut half = len >> 1; + while 0 < len { + if ss_compare_pa( + t, + sa, + pa + getidx(sa[(middle + m + half) as usize]) as isize, + pa + getidx(sa[(middle - m - half - 1) as usize]) as isize, + depth, + ) < 0 + { + m += half + 1; + half -= (len & 1) ^ 1; + } + len = half; + half >>= 1; + } + + if 0 < m { + let lm = middle - m; + let rm = middle + m; + ss_blockswap(sa, lm, middle, m); + let mut l = middle; + let mut r = middle; + let mut next = 0i32; + if rm < last { + if sa[rm as usize] < 0 { + sa[rm as usize] = !sa[rm as usize]; + if first < lm { + loop { + l -= 1; + if !(sa[l as usize] < 0) { + break; + } + } + next |= 4; + } + next |= 1; + } else if first < lm { + while sa[r as usize] < 0 { + r += 1; + } + next |= 2; + } + } + + if l - first <= last - r { + stack[ssize] = (r, rm, last, (next & 3) | (check & 4)); + ssize += 1; + middle = lm; + last = l; + check = (check & 3) | (next & 4); + } else { + if (next & 2) != 0 && r == middle { + next ^= 6; + } + stack[ssize] = (first, lm, l, (check & 3) | (next & 4)); + ssize += 1; + first = r; + middle = rm; + check = (next & 3) | (check & 4); + } + } else { + if ss_compare_pa( + t, + sa, + pa + getidx(sa[(middle - 1) as usize]) as isize, + pa + sa[middle as usize] as isize, + depth, + ) == 0 + { + sa[middle as usize] = !sa[middle as usize]; + } + ss_merge_check(t, sa, pa, first, last, check, depth); + if ssize == 0 { + return; + } + ssize -= 1; + (first, middle, last, check) = stack[ssize]; + } + } +} + +/* --------------------------------------------------------------------- */ + +/* Substring sort */ +fn sssort( + t: &[u8], + sa: &mut [i32], + pa: isize, + first: isize, + last: isize, + buf: isize, + bufsize: isize, + depth: i32, + n: isize, + lastsuffix: bool, +) { + let mut first = first; + let mut buf = buf; + let mut bufsize = bufsize; + + if lastsuffix { + first += 1; + } + + let mut limit: isize = 0; + let mut middle = last; + if bufsize < SS_BLOCKSIZE && bufsize < last - first { + limit = ss_isqrt(last - first); + if bufsize < limit { + if SS_BLOCKSIZE < limit { + limit = SS_BLOCKSIZE; + } + middle = last - limit; + buf = middle; + bufsize = limit; + } else { + limit = 0; + } + } + let mut a = first; + let mut i: isize = 0; + while SS_BLOCKSIZE < middle - a { + ss_mintrosort(t, sa, pa, a, a + SS_BLOCKSIZE, depth); + let mut curbufsize = last - (a + SS_BLOCKSIZE); + let mut curbuf = a + SS_BLOCKSIZE; + if curbufsize <= bufsize { + curbufsize = bufsize; + curbuf = buf; + } + let mut b = a; + let mut k = SS_BLOCKSIZE; + let mut j = i; + while j & 1 != 0 { + ss_swapmerge(t, sa, pa, b - k, b, b + k, curbuf, curbufsize, depth); + b -= k; + k <<= 1; + j >>= 1; + } + a += SS_BLOCKSIZE; + i += 1; + } + ss_mintrosort(t, sa, pa, a, middle, depth); + let mut k = SS_BLOCKSIZE; + while i != 0 { + if i & 1 != 0 { + ss_swapmerge(t, sa, pa, a - k, a, middle, buf, bufsize, depth); + a -= k; + } + k <<= 1; + i >>= 1; + } + if limit != 0 { + ss_mintrosort(t, sa, pa, middle, last, depth); + ss_inplacemerge(t, sa, pa, first, middle, last, depth); + } + + if lastsuffix { + /* Insert last type B* suffix. */ + let pai0 = sa[(pa + sa[(first - 1) as usize] as isize) as usize]; + let pai1 = (n - 2) as i32; + let i0 = sa[(first - 1) as usize]; + let mut a = first; + while a < last { + let av = sa[a as usize]; + if !(av < 0 + || 0 < ss_compare( + t, + pai0, + pai1, + sa[(pa + av as isize) as usize], + sa[(pa + av as isize + 1) as usize], + depth, + )) + { + break; + } + sa[(a - 1) as usize] = av; + a += 1; + } + sa[(a - 1) as usize] = i0; + } +} + +/* --------------------------------------------------------------------- */ + +#[inline] +fn tr_ilg(n: isize) -> i32 { + let n = n as i32; + if (n as u32) & 0xffff_0000 != 0 { + if (n as u32) & 0xff00_0000 != 0 { + 24 + LG_TABLE[((n >> 24) & 0xff) as usize] + } else { + 16 + LG_TABLE[((n >> 16) & 0xff) as usize] + } + } else if n & 0xff00 != 0 { + 8 + LG_TABLE[((n >> 8) & 0xff) as usize] + } else { + LG_TABLE[(n & 0xff) as usize] + } +} + +/* --------------------------------------------------------------------- */ + +/// `ISAd[SA[p]]` — the depth-offset rank of the suffix stored at `p`. +#[inline(always)] +fn tr_key(sa: &[i32], isad: isize, p: isize) -> i32 { + sa[(isad + sa[p as usize] as isize) as usize] +} + +/* Simple insertionsort for small size groups. */ +fn tr_insertionsort(sa: &mut [i32], isad: isize, first: isize, last: isize) { + let mut a = first + 1; + while a < last { + let t0 = sa[a as usize]; + let mut b = a - 1; + let mut r; + loop { + r = sa[(isad + t0 as isize) as usize] - tr_key(sa, isad, b); + if !(0 > r) { + break; + } + loop { + sa[(b + 1) as usize] = sa[b as usize]; + b -= 1; + if !(first <= b && sa[b as usize] < 0) { + break; + } + } + if b < first { + break; + } + } + if r == 0 { + sa[b as usize] = !sa[b as usize]; + } + sa[(b + 1) as usize] = t0; + a += 1; + } +} + +/* --------------------------------------------------------------------- */ + +fn tr_fixdown(sa: &mut [i32], isad: isize, base: isize, i: isize, size: isize) { + let mut i = i; + let v = sa[(base + i) as usize]; + let c = sa[(isad + v as isize) as usize]; + loop { + let mut j = 2 * i + 1; + if j >= size { + break; + } + let mut k = j; + j += 1; + let mut d = tr_key(sa, isad, base + k); + let e = tr_key(sa, isad, base + j); + if d < e { + k = j; + d = e; + } + if d <= c { + break; + } + sa[(base + i) as usize] = sa[(base + k) as usize]; + i = k; + } + sa[(base + i) as usize] = v; +} + +/* Simple top-down heapsort. */ +fn tr_heapsort(sa: &mut [i32], isad: isize, base: isize, size: isize) { + let mut m = size; + if size % 2 == 0 { + m -= 1; + if tr_key(sa, isad, base + m / 2) < tr_key(sa, isad, base + m) { + sa.swap((base + m) as usize, (base + m / 2) as usize); + } + } + + let mut i = m / 2 - 1; + while 0 <= i { + tr_fixdown(sa, isad, base, i, m); + i -= 1; + } + if size % 2 == 0 { + sa.swap(base as usize, (base + m) as usize); + tr_fixdown(sa, isad, base, 0, m); + } + let mut i = m - 1; + while 0 < i { + let t0 = sa[base as usize]; + sa[base as usize] = sa[(base + i) as usize]; + tr_fixdown(sa, isad, base, 0, i); + sa[(base + i) as usize] = t0; + i -= 1; + } +} + +/* --------------------------------------------------------------------- */ + +/* Returns the median of three elements. */ +#[inline] +fn tr_median3(sa: &[i32], isad: isize, v1: isize, v2: isize, v3: isize) -> isize { + let mut v1 = v1; + let mut v2 = v2; + if tr_key(sa, isad, v1) > tr_key(sa, isad, v2) { + std::mem::swap(&mut v1, &mut v2); + } + if tr_key(sa, isad, v2) > tr_key(sa, isad, v3) { + if tr_key(sa, isad, v1) > tr_key(sa, isad, v3) { + return v1; + } + return v3; + } + v2 +} + +/* Returns the median of five elements. */ +#[inline] +fn tr_median5( + sa: &[i32], + isad: isize, + v1: isize, + v2: isize, + v3: isize, + v4: isize, + v5: isize, +) -> isize { + let mut v1 = v1; + let mut v2 = v2; + let mut v3 = v3; + let mut v4 = v4; + let mut v5 = v5; + if tr_key(sa, isad, v2) > tr_key(sa, isad, v3) { + std::mem::swap(&mut v2, &mut v3); + } + if tr_key(sa, isad, v4) > tr_key(sa, isad, v5) { + std::mem::swap(&mut v4, &mut v5); + } + if tr_key(sa, isad, v2) > tr_key(sa, isad, v4) { + std::mem::swap(&mut v2, &mut v4); + std::mem::swap(&mut v3, &mut v5); + } + if tr_key(sa, isad, v1) > tr_key(sa, isad, v3) { + std::mem::swap(&mut v1, &mut v3); + } + if tr_key(sa, isad, v1) > tr_key(sa, isad, v4) { + std::mem::swap(&mut v1, &mut v4); + std::mem::swap(&mut v3, &mut v5); + } + if tr_key(sa, isad, v3) > tr_key(sa, isad, v4) { + return v4; + } + v3 +} + +/* Returns the pivot element. */ +#[inline] +fn tr_pivot(sa: &[i32], isad: isize, first: isize, last: isize) -> isize { + let mut t0 = last - first; + let middle = first + t0 / 2; + + if t0 <= 512 { + if t0 <= 32 { + return tr_median3(sa, isad, first, middle, last - 1); + } + t0 >>= 2; + return tr_median5(sa, isad, first, first + t0, middle, last - 1 - t0, last - 1); + } + t0 >>= 3; + let first = tr_median3(sa, isad, first, first + t0, first + (t0 << 1)); + let middle = tr_median3(sa, isad, middle - t0, middle, middle + t0); + let last = tr_median3(sa, isad, last - 1 - (t0 << 1), last - 1 - t0, last - 1); + tr_median3(sa, isad, first, middle, last) +} + +/* --------------------------------------------------------------------- */ + +struct TrBudget { + chance: i32, + remain: i32, + incval: i32, + count: i32, +} + +impl TrBudget { + fn new(chance: i32, incval: i32) -> Self { + TrBudget { + chance, + remain: incval, + incval, + count: 0, + } + } + + fn check(&mut self, size: isize) -> bool { + let size = size as i32; + if size <= self.remain { + self.remain -= size; + return true; + } + if self.chance == 0 { + self.count += size; + return false; + } + self.remain += self.incval - size; + self.chance -= 1; + true + } +} + +/* --------------------------------------------------------------------- */ + +fn tr_partition( + sa: &mut [i32], + isad: isize, + first: isize, + middle: isize, + last: isize, + v: i32, +) -> (isize, isize) { + let mut first = first; + let mut last = last; + let mut x: i32 = 0; + + let mut b = middle - 1; + loop { + b += 1; + if !(b < last) { + break; + } + x = tr_key(sa, isad, b); + if x != v { + break; + } + } + let mut a = b; + if a < last && x < v { + loop { + b += 1; + if !(b < last) { + break; + } + x = tr_key(sa, isad, b); + if !(x <= v) { + break; + } + if x == v { + sa.swap(b as usize, a as usize); + a += 1; + } + } + } + let mut c = last; + loop { + c -= 1; + if !(b < c) { + break; + } + x = tr_key(sa, isad, c); + if x != v { + break; + } + } + let mut d = c; + if b < d && x > v { + loop { + c -= 1; + if !(b < c) { + break; + } + x = tr_key(sa, isad, c); + if !(x >= v) { + break; + } + if x == v { + sa.swap(c as usize, d as usize); + d -= 1; + } + } + } + while b < c { + sa.swap(b as usize, c as usize); + loop { + b += 1; + if !(b < c) { + break; + } + x = tr_key(sa, isad, b); + if !(x <= v) { + break; + } + if x == v { + sa.swap(b as usize, a as usize); + a += 1; + } + } + loop { + c -= 1; + if !(b < c) { + break; + } + x = tr_key(sa, isad, c); + if !(x >= v) { + break; + } + if x == v { + sa.swap(c as usize, d as usize); + d -= 1; + } + } + } + + if a <= d { + c = b - 1; + let mut s = a - first; + let t0 = b - a; + if s > t0 { + s = t0; + } + let mut e = first; + let mut f = b - s; + while 0 < s { + sa.swap(e as usize, f as usize); + s -= 1; + e += 1; + f += 1; + } + let mut s = d - c; + let t0 = last - d - 1; + if s > t0 { + s = t0; + } + let mut e = b; + let mut f = last - s; + while 0 < s { + sa.swap(e as usize, f as usize); + s -= 1; + e += 1; + f += 1; + } + first += b - a; + last -= d - c; + } + (first, last) +} + +/* sort suffixes of middle partition by using sorted order of suffixes of + * left and right partition. */ +fn tr_copy( + sa: &mut [i32], + isa: isize, + first: isize, + a: isize, + b: isize, + last: isize, + depth: isize, +) { + /* All cursor arithmetic is relative to the slice start, which is the C + * routine's `SA` pointer, so `x - SA` becomes plain `x`. */ + let v = (b - 1) as i32; + + let mut c = first; + let mut d = a - 1; + while c <= d { + let s = sa[c as usize] - depth as i32; + if 0 <= s && sa[(isa + s as isize) as usize] == v { + d += 1; + sa[d as usize] = s; + sa[(isa + s as isize) as usize] = d as i32; + } + c += 1; + } + let mut c = last - 1; + let e = d + 1; + let mut d = b; + while e < d { + let s = sa[c as usize] - depth as i32; + if 0 <= s && sa[(isa + s as isize) as usize] == v { + d -= 1; + sa[d as usize] = s; + sa[(isa + s as isize) as usize] = d as i32; + } + c -= 1; + } +} + +fn tr_partialcopy( + sa: &mut [i32], + isa: isize, + first: isize, + a: isize, + b: isize, + last: isize, + depth: isize, +) { + let v = (b - 1) as i32; + let mut newrank: i32 = -1; + + let mut lastrank: i32 = -1; + let mut c = first; + let mut d = a - 1; + while c <= d { + let s = sa[c as usize] - depth as i32; + if 0 <= s && sa[(isa + s as isize) as usize] == v { + d += 1; + sa[d as usize] = s; + let rank = sa[(isa + s as isize + depth) as usize]; + if lastrank != rank { + lastrank = rank; + newrank = d as i32; + } + sa[(isa + s as isize) as usize] = newrank; + } + c += 1; + } + + let mut lastrank: i32 = -1; + let mut e = d; + while first <= e { + let rank = sa[(isa + sa[e as usize] as isize) as usize]; + if lastrank != rank { + lastrank = rank; + newrank = e as i32; + } + if newrank != rank { + sa[(isa + sa[e as usize] as isize) as usize] = newrank; + } + e -= 1; + } + + let mut lastrank: i32 = -1; + let mut c = last - 1; + let e = d + 1; + let mut d = b; + while e < d { + let s = sa[c as usize] - depth as i32; + if 0 <= s && sa[(isa + s as isize) as usize] == v { + d -= 1; + sa[d as usize] = s; + let rank = sa[(isa + s as isize + depth) as usize]; + if lastrank != rank { + lastrank = rank; + newrank = d as i32; + } + sa[(isa + s as isize) as usize] = newrank; + } + c -= 1; + } +} + +fn tr_introsort( + sa: &mut [i32], + isa: isize, + isad: isize, + first: isize, + last: isize, + budget: &mut TrBudget, +) { + /* Stack frames are (ISAd, first, last, limit, trlink); the tandem-repeat + * copy frame stores its `(a, b)` pair in the pointer fields with a zero + * placeholder where C pushes a NULL ISAd. */ + let mut stack = [(0isize, 0isize, 0isize, 0i32, 0i32); TR_STACKSIZE]; + let mut ssize = 0usize; + let mut trlink: i32 = -1; + let mut isad = isad; + let mut first = first; + let mut last = last; + let incr = isad - isa; + let mut limit = tr_ilg(last - first); + + loop { + if limit < 0 { + if limit == -1 { + /* tandem repeat partition */ + let (a, b) = tr_partition(sa, isad - incr, first, first, last, (last - 1) as i32); + + /* update ranks */ + if a < last { + let v = (a - 1) as i32; + let mut c = first; + while c < a { + sa[(isa + sa[c as usize] as isize) as usize] = v; + c += 1; + } + } + if b < last { + let v = (b - 1) as i32; + let mut c = a; + while c < b { + sa[(isa + sa[c as usize] as isize) as usize] = v; + c += 1; + } + } + + /* push */ + if 1 < b - a { + stack[ssize] = (0, a, b, 0, 0); + ssize += 1; + stack[ssize] = (isad - incr, first, last, -2, trlink); + ssize += 1; + trlink = ssize as i32 - 2; + } + if a - first <= last - b { + if 1 < a - first { + stack[ssize] = (isad, b, last, tr_ilg(last - b), trlink); + ssize += 1; + last = a; + limit = tr_ilg(a - first); + } else if 1 < last - b { + first = b; + limit = tr_ilg(last - b); + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } else if 1 < last - b { + stack[ssize] = (isad, first, a, tr_ilg(a - first), trlink); + ssize += 1; + first = b; + limit = tr_ilg(last - b); + } else if 1 < a - first { + last = a; + limit = tr_ilg(a - first); + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } else if limit == -2 { + /* tandem repeat copy */ + ssize -= 1; + let a = stack[ssize].1; + let b = stack[ssize].2; + if stack[ssize].3 == 0 { + tr_copy(sa, isa, first, a, b, last, isad - isa); + } else { + if 0 <= trlink { + stack[trlink as usize].3 = -1; + } + tr_partialcopy(sa, isa, first, a, b, last, isad - isa); + } + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } else { + /* sorted partition */ + if 0 <= sa[first as usize] { + let mut a = first; + loop { + sa[(isa + sa[a as usize] as isize) as usize] = a as i32; + a += 1; + if !(a < last && 0 <= sa[a as usize]) { + break; + } + } + first = a; + } + if first < last { + let mut a = first; + loop { + sa[a as usize] = !sa[a as usize]; + a += 1; + if !(sa[a as usize] < 0) { + break; + } + } + let next = + if sa[(isa + sa[a as usize] as isize) as usize] != tr_key(sa, isad, a) { + tr_ilg(a - first + 1) + } else { + -1 + }; + a += 1; + if a < last { + let v = (a - 1) as i32; + let mut b = first; + while b < a { + sa[(isa + sa[b as usize] as isize) as usize] = v; + b += 1; + } + } + + /* push */ + if budget.check(a - first) { + if a - first <= last - a { + stack[ssize] = (isad, a, last, -3, trlink); + ssize += 1; + isad += incr; + last = a; + limit = next; + } else if 1 < last - a { + stack[ssize] = (isad + incr, first, a, next, trlink); + ssize += 1; + first = a; + limit = -3; + } else { + isad += incr; + last = a; + limit = next; + } + } else { + if 0 <= trlink { + stack[trlink as usize].3 = -1; + } + if 1 < last - a { + first = a; + limit = -3; + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } + continue; + } + + if last - first <= TR_INSERTIONSORT_THRESHOLD { + tr_insertionsort(sa, isad, first, last); + limit = -3; + continue; + } + + /* C decrements `limit` here (`limit-- == 0`); the decrement is + * observable only on the not-taken path because the taken path + * overwrites `limit` with -3. */ + if limit == 0 { + tr_heapsort(sa, isad, first, last - first); + let mut a = last - 1; + while first < a { + let x = tr_key(sa, isad, a); + let mut b = a - 1; + while first <= b && tr_key(sa, isad, b) == x { + sa[b as usize] = !sa[b as usize]; + b -= 1; + } + a = b; + } + limit = -3; + continue; + } + limit -= 1; + + /* choose pivot */ + let a = tr_pivot(sa, isad, first, last); + sa.swap(first as usize, a as usize); + let v = tr_key(sa, isad, first); + + /* partition */ + let (a, b) = tr_partition(sa, isad, first, first + 1, last, v); + if last - first != b - a { + let next = if sa[(isa + sa[a as usize] as isize) as usize] != v { + tr_ilg(b - a) + } else { + -1 + }; + + /* update ranks */ + { + let vv = (a - 1) as i32; + let mut c = first; + while c < a { + sa[(isa + sa[c as usize] as isize) as usize] = vv; + c += 1; + } + } + if b < last { + let vv = (b - 1) as i32; + let mut c = a; + while c < b { + sa[(isa + sa[c as usize] as isize) as usize] = vv; + c += 1; + } + } + + /* push */ + if 1 < b - a && budget.check(b - a) { + if a - first <= last - b { + if last - b <= b - a { + if 1 < a - first { + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + last = a; + } else if 1 < last - b { + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + first = b; + } else { + isad += incr; + first = a; + last = b; + limit = next; + } + } else if a - first <= b - a { + if 1 < a - first { + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + last = a; + } else { + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + isad += incr; + first = a; + last = b; + limit = next; + } + } else { + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + isad += incr; + first = a; + last = b; + limit = next; + } + } else if a - first <= b - a { + if 1 < last - b { + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + first = b; + } else if 1 < a - first { + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + last = a; + } else { + isad += incr; + first = a; + last = b; + limit = next; + } + } else if last - b <= b - a { + if 1 < last - b { + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + stack[ssize] = (isad + incr, a, b, next, trlink); + ssize += 1; + first = b; + } else { + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + isad += incr; + first = a; + last = b; + limit = next; + } + } else { + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + isad += incr; + first = a; + last = b; + limit = next; + } + } else { + if 1 < b - a && 0 <= trlink { + stack[trlink as usize].3 = -1; + } + if a - first <= last - b { + if 1 < a - first { + stack[ssize] = (isad, b, last, limit, trlink); + ssize += 1; + last = a; + } else if 1 < last - b { + first = b; + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } else if 1 < last - b { + stack[ssize] = (isad, first, a, limit, trlink); + ssize += 1; + first = b; + } else if 1 < a - first { + last = a; + } else { + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } + } else if budget.check(last - first) { + limit = tr_ilg(last - first); + isad += incr; + } else { + if 0 <= trlink { + stack[trlink as usize].3 = -1; + } + if ssize == 0 { + return; + } + ssize -= 1; + (isad, first, last, limit, trlink) = stack[ssize]; + } + } +} + +/* --------------------------------------------------------------------- */ + +/* Tandem repeat sort */ +fn trsort(sa: &mut [i32], isa: isize, n: isize, depth: isize) { + let mut budget = TrBudget::new(tr_ilg(n) * 2 / 3, n as i32); + /* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */ + let mut isad = isa + depth; + while -(n as i32) < sa[0] { + let mut first: isize = 0; + let mut skip: isize = 0; + let mut unsorted: i32 = 0; + loop { + let t0 = sa[first as usize]; + if t0 < 0 { + first -= t0 as isize; + skip += t0 as isize; + } else { + if skip != 0 { + sa[(first + skip) as usize] = skip as i32; + skip = 0; + } + let last = sa[(isa + t0 as isize) as usize] as isize + 1; + if 1 < last - first { + budget.count = 0; + tr_introsort(sa, isa, isad, first, last, &mut budget); + if budget.count != 0 { + unsorted += budget.count; + } else { + skip = first - last; + } + } else if last - first == 1 { + skip = -1; + } + first = last; + } + if !(first < n) { + break; + } + } + if skip != 0 { + sa[(first + skip) as usize] = skip as i32; + } + if unsorted == 0 { + break; + } + isad += isad - isa; + } +} + +/* --------------------------------------------------------------------- */ + +/// `BUCKET_B(c0, c1)` for the 256-symbol alphabet. +#[inline(always)] +fn bb(c0: i32, c1: i32) -> usize { + (((c1 as u32) << 8) | c0 as u32) as usize +} + +/// `BUCKET_BSTAR(c0, c1)` for the 256-symbol alphabet. +#[inline(always)] +fn bstar(c0: i32, c1: i32) -> usize { + (((c0 as u32) << 8) | c1 as u32) as usize +} + +/* Sorts suffixes of type B*. */ +fn sort_type_bstar( + t: &[u8], + sa: &mut [i32], + bucket_a: &mut [i32], + bucket_b: &mut [i32], + n: isize, +) -> isize { + /* Initialize bucket arrays. */ + for slot in bucket_a.iter_mut() { + *slot = 0; + } + for slot in bucket_b.iter_mut() { + *slot = 0; + } + + /* Count the number of occurrences of the first one or two characters of + each type A, B and B* suffix. Moreover, store the beginning position of + all type B* suffixes into the array SA. */ + let mut i = n - 1; + let mut m = n; + let mut c0 = t[(n - 1) as usize] as i32; + let mut c1; + while 0 <= i { + /* type A suffix. */ + loop { + c1 = c0; + bucket_a[c1 as usize] += 1; + i -= 1; + if 0 <= i { + c0 = t[i as usize] as i32; + if c0 >= c1 { + continue; + } + } + break; + } + if 0 <= i { + /* type B* suffix. */ + bucket_b[bstar(c0, c1)] += 1; + m -= 1; + sa[m as usize] = i as i32; + /* type B suffix. */ + i -= 1; + c1 = c0; + while 0 <= i { + c0 = t[i as usize] as i32; + if !(c0 <= c1) { + break; + } + bucket_b[bb(c0, c1)] += 1; + i -= 1; + c1 = c0; + } + } + } + let m = n - m; + /* + note: + A type B* suffix is lexicographically smaller than a type B suffix that + begins with the same first two characters. + */ + + /* Calculate the index of start/end point of each bucket. */ + { + let mut i: i32 = 0; + let mut j: i32 = 0; + for c0 in 0..ALPHABET_SIZE { + let t0 = i + bucket_a[c0 as usize]; + bucket_a[c0 as usize] = i + j; /* start point */ + i = t0 + bucket_b[bb(c0, c0)]; + for c1 in (c0 + 1)..ALPHABET_SIZE { + j += bucket_b[bstar(c0, c1)]; + bucket_b[bstar(c0, c1)] = j; /* end point */ + i += bucket_b[bb(c0, c1)]; + } + } + } + + if 0 < m { + /* Sort the type B* suffixes by their first two characters. */ + let pab = n - m; + let isab = m; + let mut i = m - 2; + while 0 <= i { + let t0 = sa[(pab + i) as usize]; + let c0 = t[t0 as usize] as i32; + let c1 = t[(t0 + 1) as usize] as i32; + bucket_b[bstar(c0, c1)] -= 1; + sa[bucket_b[bstar(c0, c1)] as usize] = i as i32; + i -= 1; + } + { + let t0 = sa[(pab + m - 1) as usize]; + let c0 = t[t0 as usize] as i32; + let c1 = t[(t0 + 1) as usize] as i32; + bucket_b[bstar(c0, c1)] -= 1; + sa[bucket_b[bstar(c0, c1)] as usize] = (m - 1) as i32; + } + + /* Sort the type B* substrings using sssort. */ + let buf = m; + let bufsize = n - 2 * m; + let mut c0 = ALPHABET_SIZE - 2; + let mut j = m; + while 0 < j { + let mut c1 = ALPHABET_SIZE - 1; + while c0 < c1 { + let i = bucket_b[bstar(c0, c1)] as isize; + if 1 < j - i { + sssort( + t, + sa, + pab, + i, + j, + buf, + bufsize, + 2, + n, + sa[i as usize] == (m - 1) as i32, + ); + } + j = i; + c1 -= 1; + } + c0 -= 1; + } + + /* Compute ranks of type B* substrings. */ + let mut i = m - 1; + while 0 <= i { + if 0 <= sa[i as usize] { + let j = i; + loop { + sa[(isab + sa[i as usize] as isize) as usize] = i as i32; + i -= 1; + if !(0 <= i && 0 <= sa[i as usize]) { + break; + } + } + sa[(i + 1) as usize] = (i - j) as i32; + if i <= 0 { + break; + } + } + let j = i; + loop { + sa[i as usize] = !sa[i as usize]; + sa[(isab + sa[i as usize] as isize) as usize] = j as i32; + i -= 1; + if !(sa[i as usize] < 0) { + break; + } + } + sa[(isab + sa[i as usize] as isize) as usize] = j as i32; + i -= 1; + } + + /* Construct the inverse suffix array of type B* suffixes using + trsort. */ + trsort(sa, isab, m, 1); + + /* Set the sorted order of type B* suffixes. */ + let mut i = n - 1; + let mut j = m; + let mut c0 = t[(n - 1) as usize] as i32; + while 0 <= i { + i -= 1; + let mut c1 = c0; + while 0 <= i { + c0 = t[i as usize] as i32; + if !(c0 >= c1) { + break; + } + i -= 1; + c1 = c0; + } + if 0 <= i { + let t0 = i; + i -= 1; + c1 = c0; + while 0 <= i { + c0 = t[i as usize] as i32; + if !(c0 <= c1) { + break; + } + i -= 1; + c1 = c0; + } + j -= 1; + sa[sa[(isab + j) as usize] as usize] = if t0 == 0 || 1 < t0 - i { + t0 as i32 + } else { + !(t0 as i32) + }; + } + } + + /* Calculate the index of start/end point of each bucket. */ + bucket_b[bb(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1)] = n as i32; /* end point */ + let mut k = m - 1; + let mut c0 = ALPHABET_SIZE - 2; + while 0 <= c0 { + let mut i = bucket_a[(c0 + 1) as usize] as isize - 1; + let mut c1 = ALPHABET_SIZE - 1; + while c0 < c1 { + let t0 = i - bucket_b[bb(c0, c1)] as isize; + bucket_b[bb(c0, c1)] = i as i32; /* end point */ + + /* Move all type B* suffixes to the correct position. */ + i = t0; + let j = bucket_b[bstar(c0, c1)] as isize; + while j <= k { + sa[i as usize] = sa[k as usize]; + i -= 1; + k -= 1; + } + c1 -= 1; + } + bucket_b[bstar(c0, c0 + 1)] = (i - bucket_b[bb(c0, c0)] as isize + 1) as i32; /* start point */ + bucket_b[bb(c0, c0)] = i as i32; /* end point */ + c0 -= 1; + } + } + + m +} + +/* Constructs the suffix array by using the sorted order of type B* + * suffixes. */ +fn construct_sa( + t: &[u8], + sa: &mut [i32], + bucket_a: &mut [i32], + bucket_b: &mut [i32], + n: isize, + m: isize, +) { + if 0 < m { + /* Construct the sorted order of type B suffixes by using + the sorted order of type B* suffixes. */ + let mut c1 = ALPHABET_SIZE - 2; + while 0 <= c1 { + /* Scan the suffix array from right to left. */ + let i = bucket_b[bstar(c1, c1 + 1)] as isize; + let mut j = bucket_a[(c1 + 1) as usize] as isize - 1; + let mut k: isize = 0; + let mut c2: i32 = -1; + while i <= j { + let mut s = sa[j as usize]; + if 0 < s { + debug_assert_eq!(t[s as usize] as i32, c1); + debug_assert!((s as isize + 1) < n && t[s as usize] <= t[(s + 1) as usize]); + debug_assert!(t[(s - 1) as usize] <= t[s as usize]); + sa[j as usize] = !s; + s -= 1; + let c0 = t[s as usize] as i32; + if 0 < s && (t[(s - 1) as usize] as i32) > c0 { + s = !s; + } + if c0 != c2 { + if 0 <= c2 { + bucket_b[bb(c2, c1)] = k as i32; + } + c2 = c0; + k = bucket_b[bb(c2, c1)] as isize; + } + debug_assert!(k < j); + sa[k as usize] = s; + k -= 1; + } else { + debug_assert!((s == 0 && t[s as usize] as i32 == c1) || s < 0); + sa[j as usize] = !s; + } + j -= 1; + } + c1 -= 1; + } + } + + /* Construct the suffix array by using the sorted order of type B + suffixes. */ + let mut c2 = t[(n - 1) as usize] as i32; + let mut k = bucket_a[c2 as usize] as isize; + sa[k as usize] = if (t[(n - 2) as usize] as i32) < c2 { + !((n - 1) as i32) + } else { + (n - 1) as i32 + }; + k += 1; + /* Scan the suffix array from left to right. */ + let mut i: isize = 0; + let j = n; + while i < j { + let mut s = sa[i as usize]; + if 0 < s { + debug_assert!(t[(s - 1) as usize] >= t[s as usize]); + s -= 1; + let c0 = t[s as usize] as i32; + if s == 0 || (t[(s - 1) as usize] as i32) < c0 { + s = !s; + } + if c0 != c2 { + bucket_a[c2 as usize] = k as i32; + c2 = c0; + k = bucket_a[c2 as usize] as isize; + } + debug_assert!(i < k); + sa[k as usize] = s; + k += 1; + } else { + debug_assert!(s < 0); + sa[i as usize] = !s; + } + i += 1; + } +} + +/* --------------------------------------------------------------------- */ + +/// Rust implementation of the `divsufsort()` entry point used by +/// `ZDICT_trainFromBuffer_legacy()`. +/// +/// Integration removes the C function body, so this direct export provides +/// the existing library symbol without a wrapper. The `open_mp` parameter is +/// accepted for signature compatibility only: zstd never defines +/// `LIBBSC_OPENMP`, so the C implementation ignored it as well. +/// +/// Returns 0 on success, -1 for invalid arguments, and -2 when the bucket +/// work arrays cannot be allocated, exactly like the C routine. +#[no_mangle] +pub unsafe extern "C" fn divsufsort( + t: *const u8, + sa: *mut c_int, + n: c_int, + open_mp: c_int, +) -> c_int { + let _ = open_mp; + + /* Check arguments. */ + if t.is_null() || sa.is_null() || n < 0 { + return -1; + } + if n == 0 { + return 0; + } + + let text = unsafe { slice::from_raw_parts(t, n as usize) }; + let suffix = unsafe { slice::from_raw_parts_mut(sa, n as usize) }; + if n == 1 { + suffix[0] = 0; + return 0; + } + if n == 2 { + let m = usize::from(text[0] < text[1]); + suffix[m ^ 1] = 0; + suffix[m] = 1; + return 0; + } + + let mut bucket_a: Vec = Vec::new(); + let mut bucket_b: Vec = Vec::new(); + if bucket_a.try_reserve_exact(BUCKET_A_SIZE).is_err() + || bucket_b.try_reserve_exact(BUCKET_B_SIZE).is_err() + { + /* Match the C implementation's -2 result when malloc fails. */ + return -2; + } + bucket_a.resize(BUCKET_A_SIZE, 0); + bucket_b.resize(BUCKET_B_SIZE, 0); + + /* Suffixsort. */ + let m = sort_type_bstar(text, suffix, &mut bucket_a, &mut bucket_b, n as isize); + construct_sa(text, suffix, &mut bucket_a, &mut bucket_b, n as isize, m); + 0 +} + +#[cfg(test)] +mod tests { + use super::*; + use std::ptr; + + fn build_sa(text: &[u8]) -> Vec { + let mut sa = vec![0i32; text.len()]; + let result = unsafe { divsufsort(text.as_ptr(), sa.as_mut_ptr(), text.len() as c_int, 0) }; + assert_eq!(result, 0); + sa + } + + /// Trivial O(n^2 log n) reference: sort the suffix start positions by the + /// suffixes themselves. + fn reference_sa(text: &[u8]) -> Vec { + let mut sa: Vec = (0..text.len() as i32).collect(); + sa.sort_by(|&a, &b| text[a as usize..].cmp(&text[b as usize..])); + sa + } + + /// Suffix-array invariants: a permutation of `0..n` whose suffixes are in + /// strictly increasing lexicographic order. + fn assert_valid_sa(text: &[u8], sa: &[i32]) { + assert_eq!(sa.len(), text.len()); + let mut seen = vec![false; text.len()]; + for &p in sa { + let p = usize::try_from(p).expect("suffix index must be non-negative"); + assert!(p < text.len(), "suffix index {p} out of range"); + assert!(!seen[p], "duplicate suffix index {p}"); + seen[p] = true; + } + for pair in sa.windows(2) { + assert!( + text[pair[0] as usize..] < text[pair[1] as usize..], + "suffixes {} and {} are not in sorted order", + pair[0], + pair[1] + ); + } + } + + /// Fixed-seed numerical-recipes LCG, used to generate reproducible + /// pseudo-random sample buffers. + fn lcg_bytes(len: usize, seed: u32, alphabet: u32) -> Vec { + let mut state = seed; + (0..len) + .map(|_| { + state = state.wrapping_mul(1_664_525).wrapping_add(1_013_904_223); + ((state >> 24) % alphabet) as u8 + }) + .collect() + } + + #[test] + fn rejects_invalid_arguments() { + let text = [0u8; 1]; + let mut sa = [0i32; 1]; + assert_eq!( + unsafe { divsufsort(ptr::null(), sa.as_mut_ptr(), 1, 0) }, + -1 + ); + assert_eq!( + unsafe { divsufsort(text.as_ptr(), ptr::null_mut(), 1, 0) }, + -1 + ); + assert_eq!( + unsafe { divsufsort(text.as_ptr(), sa.as_mut_ptr(), -1, 0) }, + -1 + ); + } + + #[test] + fn sorts_trivial_inputs() { + /* empty */ + let text = [0u8; 1]; + let mut sa = [i32::MIN; 1]; + assert_eq!( + unsafe { divsufsort(text.as_ptr(), sa.as_mut_ptr(), 0, 0) }, + 0 + ); + assert_eq!(sa[0], i32::MIN, "n == 0 must not touch the output"); + + /* single byte */ + assert_eq!(build_sa(b"z"), [0]); + + /* two bytes: ascending, descending, and equal */ + assert_eq!(build_sa(b"ab"), [0, 1]); + assert_eq!(build_sa(b"ba"), [1, 0]); + assert_eq!(build_sa(b"aa"), [1, 0]); + } + + #[test] + fn sorts_all_equal_bytes() { + let text = vec![b'q'; 10_000]; + let sa = build_sa(&text); + /* For a constant text the shortest suffix sorts first. */ + let expected: Vec = (0..text.len() as i32).rev().collect(); + assert_eq!(sa, expected); + } + + #[test] + fn sorts_abracadabra_exactly() { + /* Hand-computed: a(10) abra(7) abracadabra(0) acadabra(3) adabra(5) + * bra(8) bracadabra(1) cadabra(4) dabra(6) ra(9) racadabra(2). */ + assert_eq!(build_sa(b"abracadabra"), [10, 7, 0, 3, 5, 8, 1, 4, 6, 9, 2]); + } + + #[test] + fn matches_reference_on_periodic_text() { + /* Tandem repeats exercise trsort's repeat partitioning. */ + let text: Vec = b"ab".iter().copied().cycle().take(4096).collect(); + let sa = build_sa(&text); + assert_valid_sa(&text, &sa); + assert_eq!(sa, reference_sa(&text)); + } + + #[test] + fn matches_reference_on_random_bytes() { + let text = lcg_bytes(8192, 0x0BAD_5EED, 256); + let sa = build_sa(&text); + assert_valid_sa(&text, &sa); + assert_eq!(sa, reference_sa(&text)); + } + + #[test] + fn matches_reference_on_low_alphabet_text() { + /* A four-symbol alphabet produces the large first-two-character + * buckets that reach sssort's block merging and the deeper trsort + * paths. */ + let text = lcg_bytes(16_384, 0xDEAD_BEEF, 4); + let sa = build_sa(&text); + assert_valid_sa(&text, &sa); + assert_eq!(sa, reference_sa(&text)); + } +} diff --git a/rust/src/lib.rs b/rust/src/lib.rs index c415f1531..4a2988345 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -5,6 +5,8 @@ pub mod bitstream; pub mod common; pub mod cpu; pub mod debug; +#[cfg(feature = "dict-builder")] +pub mod divsufsort; pub mod entropy_common; pub mod errors; #[cfg(feature = "compression")]