27#ifndef INCLUDED_SDSL_DIVSUFSORT
28#define INCLUDED_SDSL_DIVSUFSORT
42#if !defined(UINT8_MAX)
43#define UINT8_MAX (255)
46#define ALPHABET_SIZE (256)
47#define BUCKET_A_SIZE (ALPHABET_SIZE)
48#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
49#define SS_INSERTIONSORT_THRESHOLD (8)
50#define SS_BLOCKSIZE (1024)
51#define SS_MISORT_STACKSIZE (16)
52#define TR_INSERTIONSORT_THRESHOLD (8)
54template <
typename sa
idx_t>
60 static constexpr uint64_t TR_STACKSIZE = 64;
61 static constexpr uint64_t SS_SMERGE_STACKSIZE = 32;
67 static constexpr uint64_t TR_STACKSIZE = 96;
68 static constexpr uint64_t SS_SMERGE_STACKSIZE = 64;
71#define BUCKET_A(_c0) bucket_A[(_c0)]
72#if ALPHABET_SIZE == 256
73#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
74#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
76#define BUCKET_B(_c0, _c1) (bucket_B[(_c1)*ALPHABET_SIZE + (_c0)])
77#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0)*ALPHABET_SIZE + (_c1)])
80#define STACK_PUSH(_a, _b, _c, _d) \
82 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize++].d = (_d); \
84#define STACK_PUSH5(_a, _b, _c, _d, _e) \
86 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize].d = (_d), \
87 stack[ssize++].e = (_e); \
89#define STACK_POP(_a, _b, _c, _d) \
92 if (ssize == 0) { return; } \
93 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d; \
95#define STACK_POP5(_a, _b, _c, _d, _e) \
98 if (ssize == 0) { return; } \
99 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d, \
100 (_e) = stack[ssize].e; \
103static const int32_t lg_table[256] = {
104 -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,
105 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,
106 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,
107 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,
108 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,
109 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,
110 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
113#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
118 return (n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
119 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]);
120#elif SS_BLOCKSIZE < 256
123 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
130 return (n >> 32) ? ((n >> 48) ? ((n >> 56) ? 56 + lg_table[(n >> 56) & 0xff] : 48 + lg_table[(n >> 48) & 0xff])
131 : ((n >> 40) ? 40 + lg_table[(n >> 40) & 0xff] : 32 + lg_table[(n >> 32) & 0xff]))
132 : ((n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff]
133 : 16 + lg_table[(n >> 16) & 0xff])
134 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff]
135 : 0 + lg_table[(n >> 0) & 0xff]));
136#elif SS_BLOCKSIZE < 256
139 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
147static const int32_t sqq_table[256] = {
148 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, 64, 65, 67, 69, 71, 73,
149 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104,
150 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 128, 128,
151 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 146, 147, 148, 149,
152 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167,
153 167, 168, 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, 181, 181, 182, 183,
154 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
155 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211,
156 212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 224, 224,
157 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236,
158 237, 237, 238, 238, 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, 247, 248,
159 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255
162template <
typename sa
idx_t>
168 e = (x & 0xffff0000) ? ((x & 0xff000000) ? 24 + lg_table[(x >> 24) & 0xff] : 16 + lg_table[(x >> 16) & 0xff])
169 : ((x & 0x0000ff00) ? 8 + lg_table[(x >> 8) & 0xff] : 0 + lg_table[(x >> 0) & 0xff]);
173 y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
174 if (e >= 24) { y = (y + 1 + x / y) >> 1; }
175 y = (y + 1 + x / y) >> 1;
179 y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
183 return sqq_table[x] >> 4;
186 return (x < (y * y)) ? y - 1 : y;
192template <
typename sa
idx_t>
193inline int32_t
ss_compare(
const uint8_t * T,
const saidx_t * p1,
const saidx_t * p2, saidx_t depth)
195 const uint8_t *U1, *U2, *U1n, *U2n;
197 for (U1 = T + depth + *p1, U2 = T + depth + *p2, U1n = T + *(p1 + 1) + 2, U2n = T + *(p2 + 1) + 2;
198 (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
202 return U1 < U1n ? (U2 < U2n ? *U1 - *U2 : 1) : (U2 < U2n ? -1 : 0);
205#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
208template <
typename sa
idx_t>
209inline void ss_insertionsort(
const uint8_t * T,
const saidx_t * PA, saidx_t * first, saidx_t * last, saidx_t depth)
215 for (i = last - 2; first <= i; --i)
217 for (t = *i, j = i + 1; 0 < (r =
ss_compare(T, PA + t, PA + *j, depth));)
221 }
while ((++j < last) && (*j < 0));
222 if (last <= j) {
break; }
224 if (r == 0) { *j = ~*j; }
231#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
233template <
typename sa
idx_t>
234inline void ss_fixdown(
const uint8_t * Td,
const saidx_t * PA, saidx_t * SA, saidx_t i, saidx_t
size)
240 for (v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) <
size; SA[i] = SA[k], i = k)
242 d = Td[PA[SA[k = j++]]];
243 if (d < (e = Td[PA[SA[j]]]))
248 if (d <= c) {
break; }
254template <
typename sa
idx_t>
255inline void ss_heapsort(
const uint8_t * Td,
const saidx_t * PA, saidx_t * SA, saidx_t
size)
264 if (Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) {
std::swap(SA[m], SA[m / 2]); }
267 for (i = m / 2 - 1; 0 <= i; --i) {
ss_fixdown(Td, PA, SA, i, m); }
273 for (i = m - 1; 0 < i; --i)
275 t = SA[0], SA[0] = SA[i];
282template <
typename sa
idx_t>
283inline saidx_t *
ss_median3(
const uint8_t * Td,
const saidx_t * PA, saidx_t * v1, saidx_t * v2, saidx_t * v3)
285 if (Td[PA[*v1]] > Td[PA[*v2]]) {
std::swap(v1, v2); }
286 if (Td[PA[*v2]] > Td[PA[*v3]])
288 if (Td[PA[*v1]] > Td[PA[*v3]]) {
return v1; }
298template <
typename sa
idx_t>
300ss_median5(
const uint8_t * Td,
const saidx_t * PA, saidx_t * v1, saidx_t * v2, saidx_t * v3, saidx_t * v4, saidx_t * v5)
302 if (Td[PA[*v2]] > Td[PA[*v3]]) {
std::swap(v2, v3); }
303 if (Td[PA[*v4]] > Td[PA[*v5]]) {
std::swap(v4, v5); }
304 if (Td[PA[*v2]] > Td[PA[*v4]])
309 if (Td[PA[*v1]] > Td[PA[*v3]]) {
std::swap(v1, v3); }
310 if (Td[PA[*v1]] > Td[PA[*v4]])
315 if (Td[PA[*v3]] > Td[PA[*v4]]) {
return v4; }
320template <
typename sa
idx_t>
321inline saidx_t *
ss_pivot(
const uint8_t * Td,
const saidx_t * PA, saidx_t * first, saidx_t * last)
327 middle = first + t / 2;
331 if (t <= 32) {
return ss_median3(Td, PA, first, middle, last - 1); }
335 return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
339 first =
ss_median3(Td, PA, first, first + t, first + (t << 1));
340 middle =
ss_median3(Td, PA, middle - t, middle, middle + t);
341 last =
ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
342 return ss_median3(Td, PA, first, middle, last);
346template <
typename sa
idx_t>
347inline saidx_t *
ss_partition(
const saidx_t * PA, saidx_t * first, saidx_t * last, saidx_t depth)
351 for (a = first - 1, b = last;;)
353 for (; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
354 for (; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) {}
355 if (b <= a) {
break; }
360 if (first < a) { *first = ~*first; }
365template <
typename sa
idx_t>
366inline void ss_mintrosort(
const uint8_t * T,
const saidx_t * PA, saidx_t * first, saidx_t * last, saidx_t depth)
374 saidx_t *a, *b, *c, *d, *e, *f;
380 for (ssize = 0, limit =
ss_ilg((saidx_t)(last - first));;)
385#if 1 < SS_INSERTIONSORT_THRESHOLD
393 if (limit-- == 0) {
ss_heapsort(Td, PA, first, (saidx_t)(last - first)); }
396 for (a = first + 1, v = Td[PA[*first]]; a < last; ++a)
398 if ((x = Td[PA[*a]]) != v)
400 if (1 < (a - first)) {
break; }
405 if (Td[PA[*first] - 1] < v) { first =
ss_partition(PA, first, a, depth); }
406 if ((a - first) <= (last - a))
411 last = a, depth += 1, limit =
ss_ilg((saidx_t)(a - first));
415 first = a, limit = -1;
423 first = a, limit = -1;
427 last = a, depth += 1, limit =
ss_ilg((saidx_t)(a - first));
439 for (b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) {}
440 if (((a = b) < last) && (x < v))
442 for (; (++b < last) && ((x = Td[PA[*b]]) <= v);)
451 for (c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) {}
452 if ((b < (d = c)) && (x > v))
454 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
466 for (; (++b < c) && ((x = Td[PA[*b]]) <= v);)
474 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
488 if ((s = a - first) > (t = b - a)) { s = t; }
489 for (e = first, f = b - s; 0 < s; --s, ++e, ++f) {
std::swap(*e, *f); }
490 if ((s = d - c) > (t = last - d - 1)) { s = t; }
491 for (e = b, f = last - s; 0 < s; --s, ++e, ++f) {
std::swap(*e, *f); }
493 a = first + (b - a), c = last - (d - c);
494 b = (v <= Td[PA[*a] - 1]) ? a :
ss_partition(PA, a, c, depth);
496 if ((a - first) <= (last - c))
498 if ((last - c) <= (c - b))
504 else if ((a - first) <= (c - b))
514 first = b, last = c, depth += 1, limit =
ss_ilg((saidx_t)(c - b));
519 if ((a - first) <= (c - b))
525 else if ((last - c) <= (c - b))
535 first = b, last = c, depth += 1, limit =
ss_ilg((saidx_t)(c - b));
542 if (Td[PA[*first] - 1] < v)
545 limit =
ss_ilg((saidx_t)(last - first));
556template <
typename sa
idx_t>
560 for (; 0 < n; --n, ++a, ++b) { t = *a, *a = *b, *b = t; }
563template <
typename sa
idx_t>
564inline void ss_rotate(saidx_t * first, saidx_t * middle, saidx_t * last)
568 l = middle - first, r = last - middle;
569 for (; (0 < l) && (0 < r);)
578 a = last - 1, b = middle - 1;
581 *a-- = *b, *b-- = *a;
586 if ((r -= l + 1) <= l) {
break; }
587 a -= 1, b = middle - 1;
594 a = first, b = middle;
597 *a++ = *b, *b++ = *a;
602 if ((l -= r + 1) <= r) {
break; }
611template <
typename sa
idx_t>
630 p = PA + ~*(last - 1);
635 p = PA + *(last - 1);
637 for (a = first, len = middle - first, half = len >> 1, r = -1; 0 < len; len = half, half >>= 1)
640 q =
ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
644 half -= (len & 1) ^ 1;
653 if (r == 0) { *a = ~*a; }
657 if (first == middle) {
break; }
662 while (*--last < 0) {}
664 if (middle == last) {
break; }
669template <
typename sa
idx_t>
678 saidx_t *a, *b, *c, *bufend;
682 bufend = buf + (middle - first) - 1;
685 for (t = *(a = first), b = buf, c = middle;;)
703 *a++ = *c, *c++ = *a;
706 while (b < bufend) { *a++ = *b, *b++ = *a; }
726 *a++ = *c, *c++ = *a;
729 while (b < bufend) { *a++ = *b, *b++ = *a; }
739template <
typename sa
idx_t>
748 const saidx_t *p1, *p2;
749 saidx_t *a, *b, *c, *bufend;
754 bufend = buf + (last - middle) - 1;
767 if (*(middle - 1) < 0)
769 p2 = PA + ~*(middle - 1);
774 p2 = PA + *(middle - 1);
776 for (t = *(a = last - 1), b = bufend, c = middle - 1;;)
784 *a-- = *b, *b-- = *a;
810 *a-- = *c, *c-- = *a;
814 *a-- = *c, *c-- = *a;
817 while (buf < b) { *a-- = *b, *b-- = *a; }
836 *a-- = *b, *b-- = *a;
850 *a-- = *c, *c-- = *a;
854 *a-- = *c, *c-- = *a;
857 while (buf < b) { *a-- = *b, *b-- = *a; }
884template <
typename sa
idx_t>
894#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
895#define MERGE_CHECK(a, b, c) \
897 if (((c)&1) || (((c)&2) && (ss_compare(T, PA + GETIDX(*((a)-1)), PA + *(a), depth) == 0))) { *(a) = ~*(a); } \
898 if (((c)&4) && ((ss_compare(T, PA + GETIDX(*((b)-1)), PA + *(b), depth) == 0))) { *(b) = ~*(b); } \
905 saidx_t *l, *r, *lm, *rm;
906 saidx_t m, len, half;
910 for (check = 0, ssize = 0;;)
912 if ((last - middle) <= bufsize)
914 if ((first < middle) && (middle < last)) {
ss_mergebackward(T, PA, first, middle, last, buf, depth); }
920 if ((middle - first) <= bufsize)
922 if (first < middle) {
ss_mergeforward(T, PA, first, middle, last, buf, depth); }
928 for (m = 0, len = std::min(middle - first, last - middle), half = len >> 1; 0 < len; len = half, half >>= 1)
933 half -= (len & 1) ^ 1;
939 lm = middle - m, rm = middle + m;
941 l = r = middle, next = 0;
956 for (; *r < 0; ++r) {}
961 if ((l - first) <= (last - r))
963 STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
964 middle = lm, last = l, check = (check & 3) | (next & 4);
968 if ((next & 2) && (r == middle)) { next ^= 6; }
969 STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
970 first = r, middle = rm, check = (next & 3) | (check & 4);
975 if (
ss_compare(T, PA +
GETIDX(*(middle - 1)), PA + *middle, depth) == 0) { *middle = ~*middle; }
987template <
typename sa
idx_t>
1000 saidx_t *b, *middle, *curbuf;
1001 saidx_t j, k, curbufsize, limit;
1005 if (lastsuffix != 0) { ++first; }
1007#if SS_BLOCKSIZE == 0
1010 if ((bufsize <
SS_BLOCKSIZE) && (bufsize < (last - first)) && (bufsize < (limit =
ss_isqrt(last - first))))
1013 buf = middle = last - limit, bufsize = limit;
1017 middle = last, limit = 0;
1021#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1023#elif 1 < SS_BLOCKSIZE
1028 if (curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
1029 for (b = a, k =
SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1)
1031 ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
1034#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1036#elif 1 < SS_BLOCKSIZE
1043 ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
1049#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1051#elif 1 < SS_BLOCKSIZE
1058 if (lastsuffix != 0)
1062 PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
1063 for (a = first, i = *(first - 1); (a < last) && ((*a < 0) || (0 <
ss_compare(T, &(PAi[0]), PA + *a, depth)));
1074 return (n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
1075 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]);
1080 return (n >> 32) ? ((n >> 48) ? ((n >> 56) ? 56 + lg_table[(n >> 56) & 0xff] : 48 + lg_table[(n >> 48) & 0xff])
1081 : ((n >> 40) ? 40 + lg_table[(n >> 40) & 0xff] : 32 + lg_table[(n >> 32) & 0xff]))
1082 : ((n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff]
1083 : 16 + lg_table[(n >> 16) & 0xff])
1084 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff]
1085 : 0 + lg_table[(n >> 0) & 0xff]));
1089template <
typename sa
idx_t>
1095 for (a = first + 1; a < last; ++a)
1097 for (t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);)
1101 }
while ((first <= --b) && (*b < 0));
1102 if (b < first) {
break; }
1104 if (r == 0) { *b = ~*b; }
1109template <
typename sa
idx_t>
1116 for (v = SA[i], c = ISAd[v]; (j = 2 * i + 1) <
size; SA[i] = SA[k], i = k)
1118 d = ISAd[SA[k = j++]];
1119 if (d < (e = ISAd[SA[j]]))
1124 if (d <= c) {
break; }
1130template <
typename sa
idx_t>
1137 if ((
size % 2) == 0)
1140 if (ISAd[SA[m / 2]] < ISAd[SA[m]]) {
std::swap(SA[m], SA[m / 2]); }
1143 for (i = m / 2 - 1; 0 <= i; --i) {
tr_fixdown(ISAd, SA, i, m); }
1144 if ((
size % 2) == 0)
1149 for (i = m - 1; 0 < i; --i)
1151 t = SA[0], SA[0] = SA[i];
1158template <
typename sa
idx_t>
1159inline saidx_t *
tr_median3(
const saidx_t * ISAd, saidx_t * v1, saidx_t * v2, saidx_t * v3)
1161 if (ISAd[*v1] > ISAd[*v2]) {
std::swap(v1, v2); }
1162 if (ISAd[*v2] > ISAd[*v3])
1164 if (ISAd[*v1] > ISAd[*v3]) {
return v1; }
1174template <
typename sa
idx_t>
1175inline saidx_t *
tr_median5(
const saidx_t * ISAd, saidx_t * v1, saidx_t * v2, saidx_t * v3, saidx_t * v4, saidx_t * v5)
1177 if (ISAd[*v2] > ISAd[*v3]) {
std::swap(v2, v3); }
1178 if (ISAd[*v4] > ISAd[*v5]) {
std::swap(v4, v5); }
1179 if (ISAd[*v2] > ISAd[*v4])
1184 if (ISAd[*v1] > ISAd[*v3]) {
std::swap(v1, v3); }
1185 if (ISAd[*v1] > ISAd[*v4])
1190 if (ISAd[*v3] > ISAd[*v4]) {
return v4; }
1195template <
typename sa
idx_t>
1196inline saidx_t *
tr_pivot(
const saidx_t * ISAd, saidx_t * first, saidx_t * last)
1202 middle = first + t / 2;
1206 if (t <= 32) {
return tr_median3(ISAd, first, middle, last - 1); }
1210 return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1214 first =
tr_median3(ISAd, first, first + t, first + (t << 1));
1215 middle =
tr_median3(ISAd, middle - t, middle, middle + t);
1216 last =
tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
1217 return tr_median3(ISAd, first, middle, last);
1220template <
typename sa
idx_t>
1229template <
typename sa
idx_t>
1233template <
typename sa
idx_t>
1240template <
typename sa
idx_t>
1243 if (size <= budget->remain)
1258template <
typename sa
idx_t>
1267 saidx_t *a, *b, *c, *d, *e, *f;
1271 for (b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) {}
1272 if (((a = b) < last) && (x < v))
1274 for (; (++b < last) && ((x = ISAd[*b]) <= v);)
1283 for (c = last; (b < --c) && ((x = ISAd[*c]) == v);) {}
1284 if ((b < (d = c)) && (x > v))
1286 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1298 for (; (++b < c) && ((x = ISAd[*b]) <= v);)
1306 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1319 if ((s = a - first) > (t = b - a)) { s = t; }
1320 for (e = first, f = b - s; 0 < s; --s, ++e, ++f) {
std::swap(*e, *f); }
1321 if ((s = d - c) > (t = last - d - 1)) { s = t; }
1322 for (e = b, f = last - s; 0 < s; --s, ++e, ++f) {
std::swap(*e, *f); }
1323 first += (b - a), last -= (d - c);
1325 *pa = first, *pb = last;
1328template <
typename sa
idx_t>
1330tr_copy(saidx_t * ISA,
const saidx_t * SA, saidx_t * first, saidx_t * a, saidx_t * b, saidx_t * last, saidx_t depth)
1338 for (c = first, d = a - 1; c <= d; ++c)
1340 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1346 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1348 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1356template <
typename sa
idx_t>
1367 saidx_t rank, lastrank, newrank = -1;
1371 for (c = first, d = a - 1; c <= d; ++c)
1373 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1376 rank = ISA[s + depth];
1377 if (lastrank != rank)
1387 for (e = d; first <= e; --e)
1390 if (lastrank != rank)
1395 if (newrank != rank) { ISA[*e] = newrank; }
1399 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1401 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1404 rank = ISA[s + depth];
1405 if (lastrank != rank)
1415template <
typename sa
idx_t>
1417 const saidx_t * ISAd,
1431 saidx_t incr = ISAd - ISA;
1432 int32_t limit, next;
1433 int32_t ssize, trlink = -1;
1435 for (ssize = 0, limit =
tr_ilg((saidx_t)(last - first));;)
1443 tr_partition(ISAd - incr, first, first, last, &a, &b, (saidx_t)(last - SA - 1));
1448 for (c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1452 for (c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
1459 STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1462 if ((a - first) <= (last - b))
1464 if (1 < (a - first))
1467 last = a, limit =
tr_ilg((saidx_t)(a - first));
1469 else if (1 < (last - b))
1471 first = b, limit =
tr_ilg((saidx_t)(last - b));
1475 STACK_POP5(ISAd, first, last, limit, trlink);
1483 first = b, limit =
tr_ilg((saidx_t)(last - b));
1485 else if (1 < (a - first))
1487 last = a, limit =
tr_ilg((saidx_t)(a - first));
1491 STACK_POP5(ISAd, first, last, limit, trlink);
1495 else if (limit == -2)
1498 a = stack[--ssize].b, b = stack[ssize].c;
1499 if (stack[ssize].d == 0) {
tr_copy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA)); }
1502 if (0 <= trlink) { stack[trlink].d = -1; }
1503 tr_partialcopy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA));
1505 STACK_POP5(ISAd, first, last, limit, trlink);
1515 }
while ((++a < last) && (0 <= *a));
1524 next = (ISA[*a] != ISAd[*a]) ?
tr_ilg((saidx_t)(a - first + 1)) : -1;
1527 for (b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; }
1533 if ((a - first) <= (last - a))
1536 ISAd += incr, last = a, limit = next;
1543 first = a, limit = -3;
1547 ISAd += incr, last = a, limit = next;
1553 if (0 <= trlink) { stack[trlink].d = -1; }
1554 if (1 < (last - a)) { first = a, limit = -3; }
1557 STACK_POP5(ISAd, first, last, limit, trlink);
1563 STACK_POP5(ISAd, first, last, limit, trlink);
1578 tr_heapsort(ISAd, first, (saidx_t)(last - first));
1579 for (a = last - 1; first < a; a = b)
1581 for (x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
1594 if ((last - first) != (b - a))
1596 next = (ISA[*a] != v) ?
tr_ilg((saidx_t)(b - a)) : -1;
1599 for (c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1602 for (c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
1608 if ((a - first) <= (last - b))
1610 if ((last - b) <= (b - a))
1612 if (1 < (a - first))
1618 else if (1 < (last - b))
1625 ISAd += incr, first = a, last = b, limit = next;
1628 else if ((a - first) <= (b - a))
1630 if (1 < (a - first))
1639 ISAd += incr, first = a, last = b, limit = next;
1646 ISAd += incr, first = a, last = b, limit = next;
1651 if ((a - first) <= (b - a))
1659 else if (1 < (a - first))
1666 ISAd += incr, first = a, last = b, limit = next;
1669 else if ((last - b) <= (b - a))
1680 ISAd += incr, first = a, last = b, limit = next;
1687 ISAd += incr, first = a, last = b, limit = next;
1693 if ((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
1694 if ((a - first) <= (last - b))
1696 if (1 < (a - first))
1701 else if (1 < (last - b))
1707 STACK_POP5(ISAd, first, last, limit, trlink);
1717 else if (1 < (a - first))
1723 STACK_POP5(ISAd, first, last, limit, trlink);
1732 limit =
tr_ilg((saidx_t)(last - first)), ISAd += incr;
1736 if (0 <= trlink) { stack[trlink].d = -1; }
1737 STACK_POP5(ISAd, first, last, limit, trlink);
1746template <
typename sa
idx_t>
1747inline void trsort(saidx_t * ISA, saidx_t * SA, saidx_t n, saidx_t depth)
1750 saidx_t *first, *last;
1752 saidx_t t, skip, unsorted;
1756 for (ISAd = ISA + depth; - n < *SA; ISAd += ISAd - ISA)
1762 if ((t = *first) < 0)
1771 *(first + skip) = skip;
1774 last = SA + ISA[t] + 1;
1775 if (1 < (last - first))
1779 if (budget.
count != 0) { unsorted += budget.
count; }
1782 skip = first - last;
1785 else if ((last - first) == 1)
1791 }
while (first < (SA + n));
1792 if (skip != 0) { *(first + skip) = skip; }
1793 if (unsorted == 0) {
break; }
1798template <
typename sa
idx_t>
1799inline saidx_t
sort_typeBstar(
const uint8_t * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n)
1801 saidx_t *PAb, *ISAb, *buf;
1806 saidx_t i, j, k, t, m, bufsize;
1820 for (i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;)
1825 }
while ((0 <= --i) && ((c0 = T[i]) >= c1));
1832 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { ++
BUCKET_B(c0, c1); }
1861 for (i = m - 2; 0 <= i; --i)
1863 t = PAb[i], c0 = T[t], c1 = T[t + 1];
1866 t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
1871 tmp = omp_get_max_threads();
1872 buf = SA + m, bufsize = (n - (2 * m)) / tmp;
1874#pragma omp parallel default(shared) private(curbuf, k, l, d0, d1, tmp)
1876 tmp = omp_get_thread_num();
1877 curbuf = buf + tmp * bufsize;
1881#pragma omp critical(sssort_lock)
1891 if (--d0 < 0) {
break; }
1893 }
while (((l - k) <= 1) && (0 < (l = k)));
1894 c0 = d0, c1 = d1, j = k;
1897 if (l == 0) {
break; }
1898 sssort(T, PAb, SA + k, SA + l, curbuf, bufsize, (saidx_t)2, n, *(SA + k) == (m - 1));
1902 buf = SA + m, bufsize = n - (2 * m);
1908 if (1 < (j - i)) {
sssort(T, PAb, SA + i, SA + j, buf, bufsize, (saidx_t)2, n, *(SA + i) == (m - 1)); }
1914 for (i = m - 1; 0 <= i; --i)
1921 }
while ((0 <= --i) && (0 <= SA[i]));
1923 if (i <= 0) {
break; }
1927 ISAb[SA[i] = ~SA[i]] = j;
1928 }
while (SA[--i] < 0);
1933 trsort(ISAb, SA, m, (saidx_t)1);
1936 for (i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;)
1938 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) {}
1942 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) {}
1943 SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
1958 for (i = t, j =
BUCKET_BSTAR(c0, c1); j <= k; --i, --k) { SA[i] = SA[k]; }
1969template <
typename sa
idx_t>
1970inline void construct_SA(
const uint8_t * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n, saidx_t m)
1983 for (i = SA +
BUCKET_BSTAR(c1, c1 + 1), j = SA +
BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
1988 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1989 assert(T[s - 1] <= T[s]);
1992 if ((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1995 if (0 <= c2) {
BUCKET_B(c2, c1) = k - SA; }
2003 assert(((s == 0) && (T[s] == c1)) || (s < 0));
2013 *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
2015 for (i = SA, j = SA + n; i < j; ++i)
2019 assert(T[s - 1] >= T[s]);
2021 if ((s == 0) || (T[s - 1] < c0)) { s = ~s; }
2040template <
typename sa
idx_t>
2048 saidx_t *i, *j, *k, *orig;
2059 for (i = SA +
BUCKET_BSTAR(c1, c1 + 1), j = SA +
BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
2064 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
2065 assert(T[s - 1] <= T[s]);
2067 *j = ~((saidx_t)c0);
2068 if ((0 < s) && (T[s - 1] > c0)) { s = ~s; }
2071 if (0 <= c2) {
BUCKET_B(c2, c1) = k - SA; }
2094 *k++ = (T[n - 2] < c2) ? ~((saidx_t)T[n - 2]) : (n - 1);
2096 for (i = SA, j = SA + n, orig = SA; i < j; ++i)
2100 assert(T[s - 1] >= T[s]);
2103 if ((0 < s) && (T[s - 1] < c0)) { s = ~((saidx_t)T[s - 1]); }
2127template <
typename sa
idx_t>
2130 saidx_t *bucket_A, *bucket_B;
2135 if ((T == NULL) || (SA == NULL) || (n < 0)) {
return -1; }
2148 SA[m ^ 1] = 0, SA[m] = 1;
2152 bucket_A = (saidx_t *)malloc(
BUCKET_A_SIZE *
sizeof(saidx_t));
2153 bucket_B = (saidx_t *)malloc(
BUCKET_B_SIZE *
sizeof(saidx_t));
2156 if ((bucket_A != NULL) && (bucket_B != NULL))
2213template <
typename sa
idx_t>
2214inline int _compare(
const uint8_t * T, saidx_t Tsize,
const uint8_t * P, saidx_t Psize, saidx_t suf, saidx_t * match)
2218 for (i = suf + *match, j = *match, r = 0; (i < Tsize) && (j < Psize) && ((r = T[i] - P[j]) == 0); ++i, ++j) {}
2220 return (r == 0) ? -(j != Psize) : r;
#define STACK_POP5(_a, _b, _c, _d, _e)
#define STACK_PUSH5(_a, _b, _c, _d, _e)
#define TR_INSERTIONSORT_THRESHOLD
#define BUCKET_BSTAR(_c0, _c1)
#define SS_INSERTIONSORT_THRESHOLD
#define SS_MISORT_STACKSIZE
#define STACK_POP(_a, _b, _c, _d)
#define BUCKET_B(_c0, _c1)
#define STACK_PUSH(_a, _b, _c, _d)
#define MERGE_CHECK(a, b, c)
Namespace for the succinct data structure library.
void trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth)
saidx_t * ss_median3(const uint8_t *Td, const saidx_t *PA, saidx_t *v1, saidx_t *v2, saidx_t *v3)
saidx_t * tr_median3(const saidx_t *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3)
void ss_mergeforward(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t depth)
void tr_introsort(saidx_t *ISA, const saidx_t *ISAd, saidx_t *SA, saidx_t *first, saidx_t *last, trbudget_t< saidx_t > *budget)
saidx_t * ss_median5(const uint8_t *Td, const saidx_t *PA, saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5)
void ss_mergebackward(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t depth)
void ss_rotate(saidx_t *first, saidx_t *middle, saidx_t *last)
int32_t trbudget_check(trbudget_t< saidx_t > *budget, saidx_t size)
void ss_insertionsort(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *last, saidx_t depth)
saidx_t sort_typeBstar(const uint8_t *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n)
int32_t divsufsort(const uint8_t *T, saidx_t *SA, saidx_t n)
void ss_fixdown(const uint8_t *Td, const saidx_t *PA, saidx_t *SA, saidx_t i, saidx_t size)
void swap(int_vector_reference< t_int_vector > x, int_vector_reference< t_int_vector > y) noexcept
void ss_inplacemerge(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t depth)
void trbudget_init(trbudget_t< saidx_t > *budget, saidx_t chance, saidx_t incval)
void tr_fixdown(const saidx_t *ISAd, saidx_t *SA, saidx_t i, saidx_t size)
int _compare(const uint8_t *T, saidx_t Tsize, const uint8_t *P, saidx_t Psize, saidx_t suf, saidx_t *match)
saidx_t * tr_pivot(const saidx_t *ISAd, saidx_t *first, saidx_t *last)
void tr_copy(saidx_t *ISA, const saidx_t *SA, saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last, saidx_t depth)
saidx_t * tr_median5(const saidx_t *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5)
int32_t ss_ilg(int32_t n)
void construct_SA(const uint8_t *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n, saidx_t m)
int32_t ss_compare(const uint8_t *T, const saidx_t *p1, const saidx_t *p2, saidx_t depth)
saidx_t ss_isqrt(saidx_t x)
void ss_mintrosort(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *last, saidx_t depth)
void sssort(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *last, saidx_t *buf, saidx_t bufsize, saidx_t depth, saidx_t n, int32_t lastsuffix)
void ss_heapsort(const uint8_t *Td, const saidx_t *PA, saidx_t *SA, saidx_t size)
void tr_insertionsort(const saidx_t *ISAd, saidx_t *first, saidx_t *last)
saidx_t * ss_partition(const saidx_t *PA, saidx_t *first, saidx_t *last, saidx_t depth)
void ss_swapmerge(const uint8_t *T, const saidx_t *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t bufsize, saidx_t depth)
void tr_heapsort(const saidx_t *ISAd, saidx_t *SA, saidx_t size)
saidx_t * ss_pivot(const uint8_t *Td, const saidx_t *PA, saidx_t *first, saidx_t *last)
int32_t divsufsort64(const uint8_t *T, int64_t *SA, int64_t n)
int_vector ::size_type size(const range_type &r)
Size of a range.
void tr_partition(const saidx_t *ISAd, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t **pa, saidx_t **pb, saidx_t v)
void tr_partialcopy(saidx_t *ISA, const saidx_t *SA, saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last, saidx_t depth)
void ss_blockswap(saidx_t *a, saidx_t *b, saidx_t n)
int32_t tr_ilg(int32_t n)
saidx_t construct_BWT(const uint8_t *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n, saidx_t m)