27#ifndef INCLUDED_SDSL_DIVSUFSORT
28#define INCLUDED_SDSL_DIVSUFSORT
44#if !defined(UINT8_MAX)
45# define UINT8_MAX (255)
48#define ALPHABET_SIZE (256)
49#define BUCKET_A_SIZE (ALPHABET_SIZE)
50#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
51#define SS_INSERTIONSORT_THRESHOLD (8)
52#define SS_BLOCKSIZE (1024)
53#define SS_MISORT_STACKSIZE (16)
54#define TR_INSERTIONSORT_THRESHOLD (8)
56template <
typename sa
idx_t>
62 static constexpr uint64_t TR_STACKSIZE = 64;
63 static constexpr uint64_t SS_SMERGE_STACKSIZE = 32;
69 static constexpr uint64_t TR_STACKSIZE = 96;
70 static constexpr uint64_t SS_SMERGE_STACKSIZE = 64;
73#define BUCKET_A(_c0) bucket_A[(_c0)]
74#if ALPHABET_SIZE == 256
75# define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
76# define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
78# define BUCKET_B(_c0, _c1) (bucket_B[(_c1)*ALPHABET_SIZE + (_c0)])
79# define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0)*ALPHABET_SIZE + (_c1)])
82#define STACK_PUSH(_a, _b, _c, _d) \
85 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize++].d = (_d); \
88#define STACK_PUSH5(_a, _b, _c, _d, _e) \
91 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize].d = (_d), \
92 stack[ssize++].e = (_e); \
95#define STACK_POP(_a, _b, _c, _d) \
103 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d; \
106#define STACK_POP5(_a, _b, _c, _d, _e) \
109 assert(0 <= ssize); \
114 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d, \
115 (_e) = stack[ssize].e; \
119static const int32_t lg_table[256] = {
120 -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,
121 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,
122 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,
123 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,
124 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,
125 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,
126 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};
128#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
132# if SS_BLOCKSIZE == 0
133 return (n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
134 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]);
135# elif SS_BLOCKSIZE < 256
138 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
144# if SS_BLOCKSIZE == 0
145 return (n >> 32) ? ((n >> 48) ? ((n >> 56) ? 56 + lg_table[(n >> 56) & 0xff] : 48 + lg_table[(n >> 48) & 0xff])
146 : ((n >> 40) ? 40 + lg_table[(n >> 40) & 0xff] : 32 + lg_table[(n >> 32) & 0xff]))
148 ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
149 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]));
150# elif SS_BLOCKSIZE < 256
153 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
161static const int32_t sqq_table[256] = {
162 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, 64, 65, 67, 69, 71, 73,
163 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104,
164 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 128, 128,
165 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 146, 147, 148, 149,
166 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167,
167 167, 168, 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, 181, 181, 182, 183,
168 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
169 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211,
170 212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 224, 224,
171 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236,
172 237, 237, 238, 238, 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, 247, 248,
173 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255};
175template <
typename sa
idx_t>
184 e = (x & 0xffff0000) ? ((x & 0xff000000) ? 24 + lg_table[(x >> 24) & 0xff] : 16 + lg_table[(x >> 16) & 0xff])
185 : ((x & 0x0000ff00) ? 8 + lg_table[(x >> 8) & 0xff] : 0 + lg_table[(x >> 0) & 0xff]);
189 y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
192 y = (y + 1 + x / y) >> 1;
194 y = (y + 1 + x / y) >> 1;
198 y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
202 return sqq_table[x] >> 4;
205 return (x < (y * y)) ? y - 1 : y;
211template <
typename sa
idx_t>
212inline int32_t
ss_compare(uint8_t
const * T, saidx_t
const * p1, saidx_t
const * p2, saidx_t depth)
214 uint8_t
const *U1, *U2, *U1n, *U2n;
216 for (U1 = T + depth + *p1, U2 = T + depth + *p2, U1n = T + *(p1 + 1) + 2, U2n = T + *(p2 + 1) + 2;
217 (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
221 return U1 < U1n ? (U2 < U2n ? *U1 - *U2 : 1) : (U2 < U2n ? -1 : 0);
224#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
227template <
typename sa
idx_t>
228inline void ss_insertionsort(uint8_t
const * T, saidx_t
const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
234 for (i = last - 2; first <= i; --i)
236 for (t = *i, j = i + 1; 0 < (r =
ss_compare(T, PA + t, PA + *j, depth));)
242 while ((++j < last) && (*j < 0));
258#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
260template <
typename sa
idx_t>
261inline void ss_fixdown(uint8_t
const * Td, saidx_t
const * PA, saidx_t * SA, saidx_t i, saidx_t
size)
267 for (v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) <
size; SA[i] = SA[k], i = k)
269 d = Td[PA[SA[k = j++]]];
270 if (d < (e = Td[PA[SA[j]]]))
284template <
typename sa
idx_t>
285inline void ss_heapsort(uint8_t
const * Td, saidx_t
const * PA, saidx_t * SA, saidx_t
size)
294 if (Td[PA[SA[m / 2]]] < Td[PA[SA[m]]])
296 std::swap(SA[m], SA[m / 2]);
300 for (i = m / 2 - 1; 0 <= i; --i)
306 std::swap(SA[0], SA[m]);
309 for (i = m - 1; 0 < i; --i)
311 t = SA[0], SA[0] = SA[i];
318template <
typename sa
idx_t>
319inline saidx_t *
ss_median3(uint8_t
const * Td, saidx_t
const * PA, saidx_t * v1, saidx_t * v2, saidx_t * v3)
321 if (Td[PA[*v1]] > Td[PA[*v2]])
325 if (Td[PA[*v2]] > Td[PA[*v3]])
327 if (Td[PA[*v1]] > Td[PA[*v3]])
340template <
typename sa
idx_t>
342ss_median5(uint8_t
const * Td, saidx_t
const * PA, saidx_t * v1, saidx_t * v2, saidx_t * v3, saidx_t * v4, saidx_t * v5)
344 if (Td[PA[*v2]] > Td[PA[*v3]])
348 if (Td[PA[*v4]] > Td[PA[*v5]])
352 if (Td[PA[*v2]] > Td[PA[*v4]])
357 if (Td[PA[*v1]] > Td[PA[*v3]])
361 if (Td[PA[*v1]] > Td[PA[*v4]])
366 if (Td[PA[*v3]] > Td[PA[*v4]])
374template <
typename sa
idx_t>
375inline saidx_t *
ss_pivot(uint8_t
const * Td, saidx_t
const * PA, saidx_t * first, saidx_t * last)
381 middle = first + t / 2;
387 return ss_median3(Td, PA, first, middle, last - 1);
392 return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
396 first =
ss_median3(Td, PA, first, first + t, first + (t << 1));
397 middle =
ss_median3(Td, PA, middle - t, middle, middle + t);
398 last =
ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
399 return ss_median3(Td, PA, first, middle, last);
403template <
typename sa
idx_t>
404inline saidx_t *
ss_partition(saidx_t
const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
408 for (a = first - 1, b = last;;)
410 for (; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));)
414 for (; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));)
432template <
typename sa
idx_t>
433inline void ss_mintrosort(uint8_t
const * T, saidx_t
const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
441 saidx_t *a, *b, *c, *d, *e, *f;
447 for (ssize = 0, limit =
ss_ilg((saidx_t)(last - first));;)
452# if 1 < SS_INSERTIONSORT_THRESHOLD
453 if (1 < (last - first))
465 ss_heapsort(Td, PA, first, (saidx_t)(last - first));
469 for (a = first + 1, v = Td[PA[*first]]; a < last; ++a)
471 if ((x = Td[PA[*a]]) != v)
481 if (Td[PA[*first] - 1] < v)
485 if ((a - first) <= (last - a))
490 last = a, depth += 1, limit =
ss_ilg((saidx_t)(a - first));
494 first = a, limit = -1;
502 first = a, limit = -1;
506 last = a, depth += 1, limit =
ss_ilg((saidx_t)(a - first));
515 std::swap(*first, *a);
518 for (b = first; (++b < last) && ((x = Td[PA[*b]]) == v);)
520 if (((a = b) < last) && (x < v))
522 for (; (++b < last) && ((x = Td[PA[*b]]) <= v);)
531 for (c = last; (b < --c) && ((x = Td[PA[*c]]) == v);)
533 if ((b < (d = c)) && (x > v))
535 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
547 for (; (++b < c) && ((x = Td[PA[*b]]) <= v);)
555 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
569 if ((s = a - first) > (t = b - a))
573 for (e = first, f = b - s; 0 < s; --s, ++e, ++f)
577 if ((s = d - c) > (t = last - d - 1))
581 for (e = b, f = last - s; 0 < s; --s, ++e, ++f)
586 a = first + (b - a), c = last - (d - c);
587 b = (v <= Td[PA[*a] - 1]) ? a :
ss_partition(PA, a, c, depth);
589 if ((a - first) <= (last - c))
591 if ((last - c) <= (c - b))
597 else if ((a - first) <= (c - b))
607 first = b, last = c, depth += 1, limit =
ss_ilg((saidx_t)(c - b));
612 if ((a - first) <= (c - b))
618 else if ((last - c) <= (c - b))
628 first = b, last = c, depth += 1, limit =
ss_ilg((saidx_t)(c - b));
635 if (Td[PA[*first] - 1] < v)
638 limit =
ss_ilg((saidx_t)(last - first));
649template <
typename sa
idx_t>
653 for (; 0 < n; --n, ++a, ++b)
655 t = *a, *a = *b, *b = t;
659template <
typename sa
idx_t>
660inline void ss_rotate(saidx_t * first, saidx_t * middle, saidx_t * last)
664 l = middle - first, r = last - middle;
665 for (; (0 < l) && (0 < r);)
674 a = last - 1, b = middle - 1;
678 *a-- = *b, *b-- = *a;
683 if ((r -= l + 1) <= l)
687 a -= 1, b = middle - 1;
695 a = first, b = middle;
699 *a++ = *b, *b++ = *a;
704 if ((l -= r + 1) <= r)
717template <
typename sa
idx_t>
719ss_inplacemerge(uint8_t
const * T, saidx_t
const * PA, saidx_t * first, saidx_t * middle, saidx_t * last, saidx_t depth)
732 p = PA + ~*(last - 1);
737 p = PA + *(last - 1);
739 for (a = first, len = middle - first, half = len >> 1, r = -1; 0 < len; len = half, half >>= 1)
742 q =
ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
746 half -= (len & 1) ^ 1;
781template <
typename sa
idx_t>
790 saidx_t *a, *b, *c, *bufend;
794 bufend = buf + (middle - first) - 1;
797 for (t = *(a = first), b = buf, c = middle;;)
818 *a++ = *c, *c++ = *a;
823 *a++ = *b, *b++ = *a;
848 *a++ = *c, *c++ = *a;
853 *a++ = *b, *b++ = *a;
865template <
typename sa
idx_t>
874 saidx_t
const *p1, *p2;
875 saidx_t *a, *b, *c, *bufend;
880 bufend = buf + (last - middle) - 1;
893 if (*(middle - 1) < 0)
895 p2 = PA + ~*(middle - 1);
900 p2 = PA + *(middle - 1);
902 for (t = *(a = last - 1), b = bufend, c = middle - 1;;)
911 *a-- = *b, *b-- = *a;
939 *a-- = *c, *c-- = *a;
944 *a-- = *c, *c-- = *a;
949 *a-- = *b, *b-- = *a;
970 *a-- = *b, *b-- = *a;
986 *a-- = *c, *c-- = *a;
991 *a-- = *c, *c-- = *a;
996 *a-- = *b, *b-- = *a;
1024template <
typename sa
idx_t>
1034# define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
1035# define MERGE_CHECK(a, b, c) \
1038 if (((c)&1) || (((c)&2) && (ss_compare(T, PA + GETIDX(*((a)-1)), PA + *(a), depth) == 0))) \
1042 if (((c)&4) && ((ss_compare(T, PA + GETIDX(*((b)-1)), PA + *(b), depth) == 0))) \
1053 saidx_t *l, *r, *lm, *rm;
1054 saidx_t m, len, half;
1056 int32_t check, next;
1058 for (check = 0, ssize = 0;;)
1060 if ((last - middle) <= bufsize)
1062 if ((first < middle) && (middle < last))
1071 if ((middle - first) <= bufsize)
1082 for (m = 0, len = std::min(middle - first, last - middle), half = len >> 1; 0 < len; len = half, half >>= 1)
1087 half -= (len & 1) ^ 1;
1093 lm = middle - m, rm = middle + m;
1095 l = r = middle, next = 0;
1109 else if (first < lm)
1117 if ((l - first) <= (last - r))
1119 STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
1120 middle = lm, last = l, check = (check & 3) | (next & 4);
1124 if ((next & 2) && (r == middle))
1128 STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
1129 first = r, middle = rm, check = (next & 3) | (check & 4);
1149template <
typename sa
idx_t>
1161#if SS_BLOCKSIZE != 0
1162 saidx_t *b, *middle, *curbuf;
1163 saidx_t j, k, curbufsize, limit;
1167 if (lastsuffix != 0)
1172#if SS_BLOCKSIZE == 0
1175 if ((bufsize <
SS_BLOCKSIZE) && (bufsize < (last - first)) && (bufsize < (limit =
ss_isqrt(last - first))))
1181 buf = middle = last - limit, bufsize = limit;
1185 middle = last, limit = 0;
1189# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1191# elif 1 < SS_BLOCKSIZE
1196 if (curbufsize <= bufsize)
1198 curbufsize = bufsize, curbuf = buf;
1200 for (b = a, k =
SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1)
1202 ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
1205# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1207# elif 1 < SS_BLOCKSIZE
1214 ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
1220# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1222# elif 1 < SS_BLOCKSIZE
1229 if (lastsuffix != 0)
1233 PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
1234 for (a = first, i = *(first - 1); (a < last) && ((*a < 0) || (0 <
ss_compare(T, &(PAi[0]), PA + *a, depth)));
1245 return (n & 0xffff0000) ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
1246 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]);
1251 return (n >> 32) ? ((n >> 48) ? ((n >> 56) ? 56 + lg_table[(n >> 56) & 0xff] : 48 + lg_table[(n >> 48) & 0xff])
1252 : ((n >> 40) ? 40 + lg_table[(n >> 40) & 0xff] : 32 + lg_table[(n >> 32) & 0xff]))
1254 ? ((n & 0xff000000) ? 24 + lg_table[(n >> 24) & 0xff] : 16 + lg_table[(n >> 16) & 0xff])
1255 : ((n & 0x0000ff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff]));
1259template <
typename sa
idx_t>
1265 for (a = first + 1; a < last; ++a)
1267 for (t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);)
1273 while ((first <= --b) && (*b < 0));
1287template <
typename sa
idx_t>
1294 for (v = SA[i], c = ISAd[v]; (j = 2 * i + 1) <
size; SA[i] = SA[k], i = k)
1296 d = ISAd[SA[k = j++]];
1297 if (d < (e = ISAd[SA[j]]))
1311template <
typename sa
idx_t>
1318 if ((
size % 2) == 0)
1321 if (ISAd[SA[m / 2]] < ISAd[SA[m]])
1323 std::swap(SA[m], SA[m / 2]);
1327 for (i = m / 2 - 1; 0 <= i; --i)
1331 if ((
size % 2) == 0)
1333 std::swap(SA[0], SA[m]);
1336 for (i = m - 1; 0 < i; --i)
1338 t = SA[0], SA[0] = SA[i];
1345template <
typename sa
idx_t>
1346inline saidx_t *
tr_median3(saidx_t
const * ISAd, saidx_t * v1, saidx_t * v2, saidx_t * v3)
1348 if (ISAd[*v1] > ISAd[*v2])
1352 if (ISAd[*v2] > ISAd[*v3])
1354 if (ISAd[*v1] > ISAd[*v3])
1367template <
typename sa
idx_t>
1368inline saidx_t *
tr_median5(saidx_t
const * ISAd, saidx_t * v1, saidx_t * v2, saidx_t * v3, saidx_t * v4, saidx_t * v5)
1370 if (ISAd[*v2] > ISAd[*v3])
1374 if (ISAd[*v4] > ISAd[*v5])
1378 if (ISAd[*v2] > ISAd[*v4])
1383 if (ISAd[*v1] > ISAd[*v3])
1387 if (ISAd[*v1] > ISAd[*v4])
1392 if (ISAd[*v3] > ISAd[*v4])
1400template <
typename sa
idx_t>
1401inline saidx_t *
tr_pivot(saidx_t
const * ISAd, saidx_t * first, saidx_t * last)
1407 middle = first + t / 2;
1413 return tr_median3(ISAd, first, middle, last - 1);
1418 return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1422 first =
tr_median3(ISAd, first, first + t, first + (t << 1));
1423 middle =
tr_median3(ISAd, middle - t, middle, middle + t);
1424 last =
tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
1425 return tr_median3(ISAd, first, middle, last);
1428template <
typename sa
idx_t>
1437template <
typename sa
idx_t>
1441template <
typename sa
idx_t>
1448template <
typename sa
idx_t>
1466template <
typename sa
idx_t>
1475 saidx_t *a, *b, *c, *d, *e, *f;
1479 for (b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);)
1481 if (((a = b) < last) && (x < v))
1483 for (; (++b < last) && ((x = ISAd[*b]) <= v);)
1492 for (c = last; (b < --c) && ((x = ISAd[*c]) == v);)
1494 if ((b < (d = c)) && (x > v))
1496 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1508 for (; (++b < c) && ((x = ISAd[*b]) <= v);)
1516 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1529 if ((s = a - first) > (t = b - a))
1533 for (e = first, f = b - s; 0 < s; --s, ++e, ++f)
1537 if ((s = d - c) > (t = last - d - 1))
1541 for (e = b, f = last - s; 0 < s; --s, ++e, ++f)
1545 first += (b - a), last -= (d - c);
1547 *pa = first, *pb = last;
1550template <
typename sa
idx_t>
1552tr_copy(saidx_t * ISA, saidx_t
const * SA, saidx_t * first, saidx_t * a, saidx_t * b, saidx_t * last, saidx_t depth)
1560 for (c = first, d = a - 1; c <= d; ++c)
1562 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1568 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1570 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1578template <
typename sa
idx_t>
1589 saidx_t rank, lastrank, newrank = -1;
1593 for (c = first, d = a - 1; c <= d; ++c)
1595 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1598 rank = ISA[s + depth];
1599 if (lastrank != rank)
1609 for (e = d; first <= e; --e)
1612 if (lastrank != rank)
1617 if (newrank != rank)
1624 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1626 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1629 rank = ISA[s + depth];
1630 if (lastrank != rank)
1640template <
typename sa
idx_t>
1642 saidx_t
const * ISAd,
1656 saidx_t incr = ISAd - ISA;
1657 int32_t limit, next;
1658 int32_t ssize, trlink = -1;
1660 for (ssize = 0, limit =
tr_ilg((saidx_t)(last - first));;)
1668 tr_partition(ISAd - incr, first, first, last, &a, &b, (saidx_t)(last - SA - 1));
1673 for (c = first, v = a - SA - 1; c < a; ++c)
1680 for (c = a, v = b - SA - 1; c < b; ++c)
1690 STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1693 if ((a - first) <= (last - b))
1695 if (1 < (a - first))
1698 last = a, limit =
tr_ilg((saidx_t)(a - first));
1700 else if (1 < (last - b))
1702 first = b, limit =
tr_ilg((saidx_t)(last - b));
1706 STACK_POP5(ISAd, first, last, limit, trlink);
1714 first = b, limit =
tr_ilg((saidx_t)(last - b));
1716 else if (1 < (a - first))
1718 last = a, limit =
tr_ilg((saidx_t)(a - first));
1722 STACK_POP5(ISAd, first, last, limit, trlink);
1726 else if (limit == -2)
1729 a = stack[--ssize].b, b = stack[ssize].c;
1730 if (stack[ssize].d == 0)
1732 tr_copy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA));
1738 stack[trlink].d = -1;
1740 tr_partialcopy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA));
1742 STACK_POP5(ISAd, first, last, limit, trlink);
1754 while ((++a < last) && (0 <= *a));
1765 next = (ISA[*a] != ISAd[*a]) ?
tr_ilg((saidx_t)(a - first + 1)) : -1;
1768 for (b = first, v = a - SA - 1; b < a; ++b)
1777 if ((a - first) <= (last - a))
1780 ISAd += incr, last = a, limit = next;
1787 first = a, limit = -3;
1791 ISAd += incr, last = a, limit = next;
1799 stack[trlink].d = -1;
1803 first = a, limit = -3;
1807 STACK_POP5(ISAd, first, last, limit, trlink);
1813 STACK_POP5(ISAd, first, last, limit, trlink);
1828 tr_heapsort(ISAd, first, (saidx_t)(last - first));
1829 for (a = last - 1; first < a; a = b)
1831 for (x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b)
1842 std::swap(*first, *a);
1847 if ((last - first) != (b - a))
1849 next = (ISA[*a] != v) ?
tr_ilg((saidx_t)(b - a)) : -1;
1852 for (c = first, v = a - SA - 1; c < a; ++c)
1858 for (c = a, v = b - SA - 1; c < b; ++c)
1867 if ((a - first) <= (last - b))
1869 if ((last - b) <= (b - a))
1871 if (1 < (a - first))
1877 else if (1 < (last - b))
1884 ISAd += incr, first = a, last = b, limit = next;
1887 else if ((a - first) <= (b - a))
1889 if (1 < (a - first))
1898 ISAd += incr, first = a, last = b, limit = next;
1905 ISAd += incr, first = a, last = b, limit = next;
1910 if ((a - first) <= (b - a))
1918 else if (1 < (a - first))
1925 ISAd += incr, first = a, last = b, limit = next;
1928 else if ((last - b) <= (b - a))
1939 ISAd += incr, first = a, last = b, limit = next;
1946 ISAd += incr, first = a, last = b, limit = next;
1952 if ((1 < (b - a)) && (0 <= trlink))
1954 stack[trlink].d = -1;
1956 if ((a - first) <= (last - b))
1958 if (1 < (a - first))
1963 else if (1 < (last - b))
1969 STACK_POP5(ISAd, first, last, limit, trlink);
1979 else if (1 < (a - first))
1985 STACK_POP5(ISAd, first, last, limit, trlink);
1994 limit =
tr_ilg((saidx_t)(last - first)), ISAd += incr;
2000 stack[trlink].d = -1;
2002 STACK_POP5(ISAd, first, last, limit, trlink);
2011template <
typename sa
idx_t>
2012inline void trsort(saidx_t * ISA, saidx_t * SA, saidx_t n, saidx_t depth)
2015 saidx_t *first, *last;
2017 saidx_t t, skip, unsorted;
2021 for (ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA)
2028 if ((t = *first) < 0)
2037 *(first + skip) = skip;
2040 last = SA + ISA[t] + 1;
2041 if (1 < (last - first))
2045 if (budget.
count != 0)
2047 unsorted += budget.
count;
2051 skip = first - last;
2054 else if ((last - first) == 1)
2061 while (first < (SA + n));
2064 *(first + skip) = skip;
2074template <
typename sa
idx_t>
2075inline saidx_t
sort_typeBstar(uint8_t
const * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n)
2077 saidx_t *PAb, *ISAb, *buf;
2082 saidx_t i, j, k, t, m, bufsize;
2102 for (i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;)
2109 while ((0 <= --i) && ((c0 = T[i]) >= c1));
2116 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0)
2148 for (i = m - 2; 0 <= i; --i)
2150 t = PAb[i], c0 = T[t], c1 = T[t + 1];
2153 t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
2158 tmp = omp_get_max_threads();
2159 buf = SA + m, bufsize = (n - (2 * m)) / tmp;
2161# pragma omp parallel default(shared) private(curbuf, k, l, d0, d1, tmp)
2163 tmp = omp_get_thread_num();
2164 curbuf = buf + tmp * bufsize;
2168# pragma omp critical(sssort_lock)
2185 while (((l - k) <= 1) && (0 < (l = k)));
2186 c0 = d0, c1 = d1, j = k;
2193 sssort(T, PAb, SA + k, SA + l, curbuf, bufsize, (saidx_t)2, n, *(SA + k) == (m - 1));
2197 buf = SA + m, bufsize = n - (2 * m);
2205 sssort(T, PAb, SA + i, SA + j, buf, bufsize, (saidx_t)2, n, *(SA + i) == (m - 1));
2212 for (i = m - 1; 0 <= i; --i)
2221 while ((0 <= --i) && (0 <= SA[i]));
2231 ISAb[SA[i] = ~SA[i]] = j;
2233 while (SA[--i] < 0);
2238 trsort(ISAb, SA, m, (saidx_t)1);
2241 for (i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;)
2243 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0)
2248 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0)
2250 SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
2265 for (i = t, j =
BUCKET_BSTAR(c0, c1); j <= k; --i, --k)
2279template <
typename sa
idx_t>
2280inline void construct_SA(uint8_t
const * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n, saidx_t m)
2293 for (i = SA +
BUCKET_BSTAR(c1, c1 + 1), j = SA +
BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
2298 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
2299 assert(T[s - 1] <= T[s]);
2302 if ((0 < s) && (T[s - 1] > c0))
2319 assert(((s == 0) && (T[s] == c1)) || (s < 0));
2329 *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
2331 for (i = SA, j = SA + n; i < j; ++i)
2335 assert(T[s - 1] >= T[s]);
2337 if ((s == 0) || (T[s - 1] < c0))
2359template <
typename sa
idx_t>
2361construct_BWT(uint8_t
const * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n, saidx_t m)
2363 saidx_t *i, *j, *k, *orig;
2374 for (i = SA +
BUCKET_BSTAR(c1, c1 + 1), j = SA +
BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
2379 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
2380 assert(T[s - 1] <= T[s]);
2382 *j = ~((saidx_t)c0);
2383 if ((0 < s) && (T[s - 1] > c0))
2415 *k++ = (T[n - 2] < c2) ? ~((saidx_t)T[n - 2]) : (n - 1);
2417 for (i = SA, j = SA + n, orig = SA; i < j; ++i)
2421 assert(T[s - 1] >= T[s]);
2424 if ((0 < s) && (T[s - 1] < c0))
2426 s = ~((saidx_t)T[s - 1]);
2451template <
typename sa
idx_t>
2454 saidx_t *bucket_A, *bucket_B;
2459 if ((T == NULL) || (SA == NULL) || (n < 0))
2475 SA[m ^ 1] = 0, SA[m] = 1;
2479 bucket_A = (saidx_t *)malloc(
BUCKET_A_SIZE *
sizeof(saidx_t));
2480 bucket_B = (saidx_t *)malloc(
BUCKET_B_SIZE *
sizeof(saidx_t));
2483 if ((bucket_A != NULL) && (bucket_B != NULL))
2540template <
typename sa
idx_t>
2541inline int _compare(uint8_t
const * T, saidx_t Tsize, uint8_t
const * P, saidx_t Psize, saidx_t suf, saidx_t * match)
2545 for (i = suf + *match, j = *match, r = 0; (i < Tsize) && (j < Psize) && ((r = T[i] - P[j]) == 0); ++i, ++j)
2548 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 sort_typeBstar(uint8_t const *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n)
void tr_partialcopy(saidx_t *ISA, saidx_t const *SA, saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last, saidx_t depth)
void tr_insertionsort(saidx_t const *ISAd, saidx_t *first, saidx_t *last)
void tr_heapsort(saidx_t const *ISAd, saidx_t *SA, saidx_t size)
void ss_mergeforward(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t depth)
void tr_copy(saidx_t *ISA, saidx_t const *SA, saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last, saidx_t depth)
void ss_fixdown(uint8_t const *Td, saidx_t const *PA, saidx_t *SA, saidx_t i, saidx_t size)
saidx_t * ss_pivot(uint8_t const *Td, saidx_t const *PA, saidx_t *first, saidx_t *last)
int32_t divsufsort(uint8_t const *T, saidx_t *SA, saidx_t n)
void ss_rotate(saidx_t *first, saidx_t *middle, saidx_t *last)
void ss_mintrosort(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *last, saidx_t depth)
void sssort(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *last, saidx_t *buf, saidx_t bufsize, saidx_t depth, saidx_t n, int32_t lastsuffix)
int32_t trbudget_check(trbudget_t< saidx_t > *budget, saidx_t size)
saidx_t * tr_median5(saidx_t const *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5)
void tr_introsort(saidx_t *ISA, saidx_t const *ISAd, saidx_t *SA, saidx_t *first, saidx_t *last, trbudget_t< saidx_t > *budget)
void construct_SA(uint8_t const *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n, saidx_t m)
void ss_swapmerge(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t bufsize, saidx_t depth)
void ss_insertionsort(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *last, saidx_t depth)
saidx_t * ss_partition(saidx_t const *PA, saidx_t *first, saidx_t *last, saidx_t depth)
void ss_mergebackward(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t *buf, saidx_t depth)
void trbudget_init(trbudget_t< saidx_t > *budget, saidx_t chance, saidx_t incval)
void ss_heapsort(uint8_t const *Td, saidx_t const *PA, saidx_t *SA, saidx_t size)
int32_t divsufsort64(uint8_t const *T, int64_t *SA, int64_t n)
void ss_inplacemerge(uint8_t const *T, saidx_t const *PA, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t depth)
void tr_partition(saidx_t const *ISAd, saidx_t *first, saidx_t *middle, saidx_t *last, saidx_t **pa, saidx_t **pb, saidx_t v)
saidx_t construct_BWT(uint8_t const *T, saidx_t *SA, saidx_t *bucket_A, saidx_t *bucket_B, saidx_t n, saidx_t m)
saidx_t * ss_median5(uint8_t const *Td, saidx_t const *PA, saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5)
saidx_t * tr_pivot(saidx_t const *ISAd, saidx_t *first, saidx_t *last)
saidx_t * tr_median3(saidx_t const *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3)
int32_t ss_ilg(int32_t n)
int_vector ::size_type size(range_type const &r)
Size of a range.
saidx_t ss_isqrt(saidx_t x)
void tr_fixdown(saidx_t const *ISAd, saidx_t *SA, saidx_t i, saidx_t size)
int32_t ss_compare(uint8_t const *T, saidx_t const *p1, saidx_t const *p2, saidx_t depth)
saidx_t * ss_median3(uint8_t const *Td, saidx_t const *PA, saidx_t *v1, saidx_t *v2, saidx_t *v3)
int _compare(uint8_t const *T, saidx_t Tsize, uint8_t const *P, saidx_t Psize, saidx_t suf, saidx_t *match)
void ss_blockswap(saidx_t *a, saidx_t *b, saidx_t n)
int32_t tr_ilg(int32_t n)