SDSL 3.0.2
Succinct Data Structure Library
Loading...
Searching...
No Matches
divsufsort.hpp
Go to the documentation of this file.
1/*
2 * libdivsufsort
3 * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
4 *
5 * Permission is hereby granted, free of charge, to any person
6 * obtaining a copy of this software and associated documentation
7 * files (the "Software"), to deal in the Software without
8 * restriction, including without limitation the rights to use,
9 * copy, modify, merge, publish, distribute, sublicense, and/or sell
10 * copies of the Software, and to permit persons to whom the
11 * Software is furnished to do so, subject to the following
12 * conditions:
13 *
14 * The above copyright notice and this permission notice shall be
15 * included in all copies or substantial portions of the Software.
16 *
17 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 * OTHER DEALINGS IN THE SOFTWARE.
25 */
26
27#ifndef INCLUDED_SDSL_DIVSUFSORT
28#define INCLUDED_SDSL_DIVSUFSORT
29
30#include <algorithm>
31#include <assert.h>
32#include <inttypes.h>
33#include <stdio.h>
34#include <stdlib.h>
35#include <utility>
36
37#ifdef _OPENMP
38# include <omp.h>
39#endif
40
41namespace sdsl
42{
43
44#if !defined(UINT8_MAX)
45# define UINT8_MAX (255)
46#endif
47
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)
55
56template <typename saidx_t>
58
59template <>
60struct libdivsufsort_config<int32_t>
61{
62 static constexpr uint64_t TR_STACKSIZE = 64;
63 static constexpr uint64_t SS_SMERGE_STACKSIZE = 32;
64};
65
66template <>
67struct libdivsufsort_config<int64_t>
68{
69 static constexpr uint64_t TR_STACKSIZE = 96;
70 static constexpr uint64_t SS_SMERGE_STACKSIZE = 64;
71};
72
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)])
77#else
78# define BUCKET_B(_c0, _c1) (bucket_B[(_c1)*ALPHABET_SIZE + (_c0)])
79# define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0)*ALPHABET_SIZE + (_c1)])
80#endif
81
82#define STACK_PUSH(_a, _b, _c, _d) \
83 do \
84 { \
85 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize++].d = (_d); \
86 } \
87 while (0)
88#define STACK_PUSH5(_a, _b, _c, _d, _e) \
89 do \
90 { \
91 stack[ssize].a = (_a), stack[ssize].b = (_b), stack[ssize].c = (_c), stack[ssize].d = (_d), \
92 stack[ssize++].e = (_e); \
93 } \
94 while (0)
95#define STACK_POP(_a, _b, _c, _d) \
96 do \
97 { \
98 assert(0 <= ssize); \
99 if (ssize == 0) \
100 { \
101 return; \
102 } \
103 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d; \
104 } \
105 while (0)
106#define STACK_POP5(_a, _b, _c, _d, _e) \
107 do \
108 { \
109 assert(0 <= ssize); \
110 if (ssize == 0) \
111 { \
112 return; \
113 } \
114 (_a) = stack[--ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d, \
115 (_e) = stack[ssize].e; \
116 } \
117 while (0)
118
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};
127
128#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
129
130inline int32_t ss_ilg(int32_t n)
131{
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
136 return lg_table[n];
137# else
138 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
139# endif
140}
141
142inline int32_t ss_ilg(int64_t n)
143{
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]))
147 : ((n & 0xffff0000)
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
151 return lg_table[n];
152# else
153 return (n & 0xff00) ? 8 + lg_table[(n >> 8) & 0xff] : 0 + lg_table[(n >> 0) & 0xff];
154# endif
155}
156
157#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
158
159#if SS_BLOCKSIZE != 0
160
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};
174
175template <typename saidx_t>
176inline saidx_t ss_isqrt(saidx_t x)
177{
178 saidx_t y, e;
179
180 if (x >= (SS_BLOCKSIZE * SS_BLOCKSIZE))
181 {
182 return SS_BLOCKSIZE;
183 }
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]);
186
187 if (e >= 16)
188 {
189 y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
190 if (e >= 24)
191 {
192 y = (y + 1 + x / y) >> 1;
193 }
194 y = (y + 1 + x / y) >> 1;
195 }
196 else if (e >= 8)
197 {
198 y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
199 }
200 else
201 {
202 return sqq_table[x] >> 4;
203 }
204
205 return (x < (y * y)) ? y - 1 : y;
206}
207
208#endif /* SS_BLOCKSIZE != 0 */
209
210/* Compares two suffixes. */
211template <typename saidx_t>
212inline int32_t ss_compare(uint8_t const * T, saidx_t const * p1, saidx_t const * p2, saidx_t depth)
213{
214 uint8_t const *U1, *U2, *U1n, *U2n;
215
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);
218 ++U1, ++U2)
219 {}
220
221 return U1 < U1n ? (U2 < U2n ? *U1 - *U2 : 1) : (U2 < U2n ? -1 : 0);
222}
223
224#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
225
226/* Insertionsort for small size groups */
227template <typename saidx_t>
228inline void ss_insertionsort(uint8_t const * T, saidx_t const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
229{
230 saidx_t *i, *j;
231 saidx_t t;
232 int32_t r;
233
234 for (i = last - 2; first <= i; --i)
235 {
236 for (t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));)
237 {
238 do
239 {
240 *(j - 1) = *j;
241 }
242 while ((++j < last) && (*j < 0));
243 if (last <= j)
244 {
245 break;
246 }
247 }
248 if (r == 0)
249 {
250 *j = ~*j;
251 }
252 *(j - 1) = t;
253 }
254}
255
256#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */
257
258#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
259
260template <typename saidx_t>
261inline void ss_fixdown(uint8_t const * Td, saidx_t const * PA, saidx_t * SA, saidx_t i, saidx_t size)
262{
263 saidx_t j, k;
264 saidx_t v;
265 int32_t c, d, e;
266
267 for (v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k)
268 {
269 d = Td[PA[SA[k = j++]]];
270 if (d < (e = Td[PA[SA[j]]]))
271 {
272 k = j;
273 d = e;
274 }
275 if (d <= c)
276 {
277 break;
278 }
279 }
280 SA[i] = v;
281}
282
283/* Simple top-down heapsort. */
284template <typename saidx_t>
285inline void ss_heapsort(uint8_t const * Td, saidx_t const * PA, saidx_t * SA, saidx_t size)
286{
287 saidx_t i, m;
288 saidx_t t;
289
290 m = size;
291 if ((size % 2) == 0)
292 {
293 m--;
294 if (Td[PA[SA[m / 2]]] < Td[PA[SA[m]]])
295 {
296 std::swap(SA[m], SA[m / 2]);
297 }
298 }
299
300 for (i = m / 2 - 1; 0 <= i; --i)
301 {
302 ss_fixdown(Td, PA, SA, i, m);
303 }
304 if ((size % 2) == 0)
305 {
306 std::swap(SA[0], SA[m]);
307 ss_fixdown(Td, PA, SA, (saidx_t)0, m);
308 }
309 for (i = m - 1; 0 < i; --i)
310 {
311 t = SA[0], SA[0] = SA[i];
312 ss_fixdown(Td, PA, SA, (saidx_t)0, i);
313 SA[i] = t;
314 }
315}
316
317/* Returns the median of three elements. */
318template <typename saidx_t>
319inline saidx_t * ss_median3(uint8_t const * Td, saidx_t const * PA, saidx_t * v1, saidx_t * v2, saidx_t * v3)
320{
321 if (Td[PA[*v1]] > Td[PA[*v2]])
322 {
323 std::swap(v1, v2);
324 }
325 if (Td[PA[*v2]] > Td[PA[*v3]])
326 {
327 if (Td[PA[*v1]] > Td[PA[*v3]])
328 {
329 return v1;
330 }
331 else
332 {
333 return v3;
334 }
335 }
336 return v2;
337}
338
339/* Returns the median of five elements. */
340template <typename saidx_t>
341inline saidx_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)
343{
344 if (Td[PA[*v2]] > Td[PA[*v3]])
345 {
346 std::swap(v2, v3);
347 }
348 if (Td[PA[*v4]] > Td[PA[*v5]])
349 {
350 std::swap(v4, v5);
351 }
352 if (Td[PA[*v2]] > Td[PA[*v4]])
353 {
354 std::swap(v2, v4);
355 std::swap(v3, v5);
356 }
357 if (Td[PA[*v1]] > Td[PA[*v3]])
358 {
359 std::swap(v1, v3);
360 }
361 if (Td[PA[*v1]] > Td[PA[*v4]])
362 {
363 std::swap(v1, v4);
364 std::swap(v3, v5);
365 }
366 if (Td[PA[*v3]] > Td[PA[*v4]])
367 {
368 return v4;
369 }
370 return v3;
371}
372
373/* Returns the pivot element. */
374template <typename saidx_t>
375inline saidx_t * ss_pivot(uint8_t const * Td, saidx_t const * PA, saidx_t * first, saidx_t * last)
376{
377 saidx_t * middle;
378 saidx_t t;
379
380 t = last - first;
381 middle = first + t / 2;
382
383 if (t <= 512)
384 {
385 if (t <= 32)
386 {
387 return ss_median3(Td, PA, first, middle, last - 1);
388 }
389 else
390 {
391 t >>= 2;
392 return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
393 }
394 }
395 t >>= 3;
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);
400}
401
402/* Binary partition for substrings. */
403template <typename saidx_t>
404inline saidx_t * ss_partition(saidx_t const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
405{
406 saidx_t *a, *b;
407 saidx_t t;
408 for (a = first - 1, b = last;;)
409 {
410 for (; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));)
411 {
412 *a = ~*a;
413 }
414 for (; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));)
415 {}
416 if (b <= a)
417 {
418 break;
419 }
420 t = ~*b;
421 *b = *a;
422 *a = t;
423 }
424 if (first < a)
425 {
426 *first = ~*first;
427 }
428 return a;
429}
430
431/* Multikey introsort for medium size groups. */
432template <typename saidx_t>
433inline void ss_mintrosort(uint8_t const * T, saidx_t const * PA, saidx_t * first, saidx_t * last, saidx_t depth)
434{
435 struct
436 {
437 saidx_t *a, *b, c;
438 int32_t d;
439 } stack[SS_MISORT_STACKSIZE];
440 uint8_t const * Td;
441 saidx_t *a, *b, *c, *d, *e, *f;
442 saidx_t s, t;
443 int32_t ssize;
444 int32_t limit;
445 int32_t v, x = 0;
446
447 for (ssize = 0, limit = ss_ilg((saidx_t)(last - first));;)
448 {
449
450 if ((last - first) <= SS_INSERTIONSORT_THRESHOLD)
451 {
452# if 1 < SS_INSERTIONSORT_THRESHOLD
453 if (1 < (last - first))
454 {
455 ss_insertionsort(T, PA, first, last, depth);
456 }
457# endif
458 STACK_POP(first, last, depth, limit);
459 continue;
460 }
461
462 Td = T + depth;
463 if (limit-- == 0)
464 {
465 ss_heapsort(Td, PA, first, (saidx_t)(last - first));
466 }
467 if (limit < 0)
468 {
469 for (a = first + 1, v = Td[PA[*first]]; a < last; ++a)
470 {
471 if ((x = Td[PA[*a]]) != v)
472 {
473 if (1 < (a - first))
474 {
475 break;
476 }
477 v = x;
478 first = a;
479 }
480 }
481 if (Td[PA[*first] - 1] < v)
482 {
483 first = ss_partition(PA, first, a, depth);
484 }
485 if ((a - first) <= (last - a))
486 {
487 if (1 < (a - first))
488 {
489 STACK_PUSH(a, last, depth, -1);
490 last = a, depth += 1, limit = ss_ilg((saidx_t)(a - first));
491 }
492 else
493 {
494 first = a, limit = -1;
495 }
496 }
497 else
498 {
499 if (1 < (last - a))
500 {
501 STACK_PUSH(first, a, depth + 1, ss_ilg((saidx_t)(a - first)));
502 first = a, limit = -1;
503 }
504 else
505 {
506 last = a, depth += 1, limit = ss_ilg((saidx_t)(a - first));
507 }
508 }
509 continue;
510 }
511
512 /* choose pivot */
513 a = ss_pivot(Td, PA, first, last);
514 v = Td[PA[*a]];
515 std::swap(*first, *a);
516
517 /* partition */
518 for (b = first; (++b < last) && ((x = Td[PA[*b]]) == v);)
519 {}
520 if (((a = b) < last) && (x < v))
521 {
522 for (; (++b < last) && ((x = Td[PA[*b]]) <= v);)
523 {
524 if (x == v)
525 {
526 std::swap(*b, *a);
527 ++a;
528 }
529 }
530 }
531 for (c = last; (b < --c) && ((x = Td[PA[*c]]) == v);)
532 {}
533 if ((b < (d = c)) && (x > v))
534 {
535 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
536 {
537 if (x == v)
538 {
539 std::swap(*c, *d);
540 --d;
541 }
542 }
543 }
544 for (; b < c;)
545 {
546 std::swap(*b, *c);
547 for (; (++b < c) && ((x = Td[PA[*b]]) <= v);)
548 {
549 if (x == v)
550 {
551 std::swap(*b, *a);
552 ++a;
553 }
554 }
555 for (; (b < --c) && ((x = Td[PA[*c]]) >= v);)
556 {
557 if (x == v)
558 {
559 std::swap(*c, *d);
560 --d;
561 }
562 }
563 }
564
565 if (a <= d)
566 {
567 c = b - 1;
568
569 if ((s = a - first) > (t = b - a))
570 {
571 s = t;
572 }
573 for (e = first, f = b - s; 0 < s; --s, ++e, ++f)
574 {
575 std::swap(*e, *f);
576 }
577 if ((s = d - c) > (t = last - d - 1))
578 {
579 s = t;
580 }
581 for (e = b, f = last - s; 0 < s; --s, ++e, ++f)
582 {
583 std::swap(*e, *f);
584 }
585
586 a = first + (b - a), c = last - (d - c);
587 b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth);
588
589 if ((a - first) <= (last - c))
590 {
591 if ((last - c) <= (c - b))
592 {
593 STACK_PUSH(b, c, depth + 1, ss_ilg((saidx_t)(c - b)));
594 STACK_PUSH(c, last, depth, limit);
595 last = a;
596 }
597 else if ((a - first) <= (c - b))
598 {
599 STACK_PUSH(c, last, depth, limit);
600 STACK_PUSH(b, c, depth + 1, ss_ilg((saidx_t)(c - b)));
601 last = a;
602 }
603 else
604 {
605 STACK_PUSH(c, last, depth, limit);
606 STACK_PUSH(first, a, depth, limit);
607 first = b, last = c, depth += 1, limit = ss_ilg((saidx_t)(c - b));
608 }
609 }
610 else
611 {
612 if ((a - first) <= (c - b))
613 {
614 STACK_PUSH(b, c, depth + 1, ss_ilg((saidx_t)(c - b)));
615 STACK_PUSH(first, a, depth, limit);
616 first = c;
617 }
618 else if ((last - c) <= (c - b))
619 {
620 STACK_PUSH(first, a, depth, limit);
621 STACK_PUSH(b, c, depth + 1, ss_ilg((saidx_t)(c - b)));
622 first = c;
623 }
624 else
625 {
626 STACK_PUSH(first, a, depth, limit);
627 STACK_PUSH(c, last, depth, limit);
628 first = b, last = c, depth += 1, limit = ss_ilg((saidx_t)(c - b));
629 }
630 }
631 }
632 else
633 {
634 limit += 1;
635 if (Td[PA[*first] - 1] < v)
636 {
637 first = ss_partition(PA, first, last, depth);
638 limit = ss_ilg((saidx_t)(last - first));
639 }
640 depth += 1;
641 }
642 }
643}
644
645#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
646
647#if SS_BLOCKSIZE != 0
648
649template <typename saidx_t>
650inline void ss_blockswap(saidx_t * a, saidx_t * b, saidx_t n)
651{
652 saidx_t t;
653 for (; 0 < n; --n, ++a, ++b)
654 {
655 t = *a, *a = *b, *b = t;
656 }
657}
658
659template <typename saidx_t>
660inline void ss_rotate(saidx_t * first, saidx_t * middle, saidx_t * last)
661{
662 saidx_t *a, *b, t;
663 saidx_t l, r;
664 l = middle - first, r = last - middle;
665 for (; (0 < l) && (0 < r);)
666 {
667 if (l == r)
668 {
669 ss_blockswap(first, middle, l);
670 break;
671 }
672 if (l < r)
673 {
674 a = last - 1, b = middle - 1;
675 t = *a;
676 do
677 {
678 *a-- = *b, *b-- = *a;
679 if (b < first)
680 {
681 *a = t;
682 last = a;
683 if ((r -= l + 1) <= l)
684 {
685 break;
686 }
687 a -= 1, b = middle - 1;
688 t = *a;
689 }
690 }
691 while (1);
692 }
693 else
694 {
695 a = first, b = middle;
696 t = *a;
697 do
698 {
699 *a++ = *b, *b++ = *a;
700 if (last <= b)
701 {
702 *a = t;
703 first = a + 1;
704 if ((l -= r + 1) <= r)
705 {
706 break;
707 }
708 a += 1, b = middle;
709 t = *a;
710 }
711 }
712 while (1);
713 }
714 }
715}
716
717template <typename saidx_t>
718inline void
719ss_inplacemerge(uint8_t const * T, saidx_t const * PA, saidx_t * first, saidx_t * middle, saidx_t * last, saidx_t depth)
720{
721 saidx_t const * p;
722 saidx_t *a, *b;
723 saidx_t len, half;
724 int32_t q, r;
725 int32_t x;
726
727 for (;;)
728 {
729 if (*(last - 1) < 0)
730 {
731 x = 1;
732 p = PA + ~*(last - 1);
733 }
734 else
735 {
736 x = 0;
737 p = PA + *(last - 1);
738 }
739 for (a = first, len = middle - first, half = len >> 1, r = -1; 0 < len; len = half, half >>= 1)
740 {
741 b = a + half;
742 q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
743 if (q < 0)
744 {
745 a = b + 1;
746 half -= (len & 1) ^ 1;
747 }
748 else
749 {
750 r = q;
751 }
752 }
753 if (a < middle)
754 {
755 if (r == 0)
756 {
757 *a = ~*a;
758 }
759 ss_rotate(a, middle, last);
760 last -= middle - a;
761 middle = a;
762 if (first == middle)
763 {
764 break;
765 }
766 }
767 --last;
768 if (x != 0)
769 {
770 while (*--last < 0)
771 {}
772 }
773 if (middle == last)
774 {
775 break;
776 }
777 }
778}
779
780/* Merge-forward with internal buffer. */
781template <typename saidx_t>
782inline void ss_mergeforward(uint8_t const * T,
783 saidx_t const * PA,
784 saidx_t * first,
785 saidx_t * middle,
786 saidx_t * last,
787 saidx_t * buf,
788 saidx_t depth)
789{
790 saidx_t *a, *b, *c, *bufend;
791 saidx_t t;
792 int32_t r;
793
794 bufend = buf + (middle - first) - 1;
795 ss_blockswap(buf, first, (saidx_t)(middle - first));
796
797 for (t = *(a = first), b = buf, c = middle;;)
798 {
799 r = ss_compare(T, PA + *b, PA + *c, depth);
800 if (r < 0)
801 {
802 do
803 {
804 *a++ = *b;
805 if (bufend <= b)
806 {
807 *bufend = t;
808 return;
809 }
810 *b++ = *a;
811 }
812 while (*b < 0);
813 }
814 else if (r > 0)
815 {
816 do
817 {
818 *a++ = *c, *c++ = *a;
819 if (last <= c)
820 {
821 while (b < bufend)
822 {
823 *a++ = *b, *b++ = *a;
824 }
825 *a = *b, *b = t;
826 return;
827 }
828 }
829 while (*c < 0);
830 }
831 else
832 {
833 *c = ~*c;
834 do
835 {
836 *a++ = *b;
837 if (bufend <= b)
838 {
839 *bufend = t;
840 return;
841 }
842 *b++ = *a;
843 }
844 while (*b < 0);
845
846 do
847 {
848 *a++ = *c, *c++ = *a;
849 if (last <= c)
850 {
851 while (b < bufend)
852 {
853 *a++ = *b, *b++ = *a;
854 }
855 *a = *b, *b = t;
856 return;
857 }
858 }
859 while (*c < 0);
860 }
861 }
862}
863
864/* Merge-backward with internal buffer. */
865template <typename saidx_t>
866inline void ss_mergebackward(uint8_t const * T,
867 saidx_t const * PA,
868 saidx_t * first,
869 saidx_t * middle,
870 saidx_t * last,
871 saidx_t * buf,
872 saidx_t depth)
873{
874 saidx_t const *p1, *p2;
875 saidx_t *a, *b, *c, *bufend;
876 saidx_t t;
877 int32_t r;
878 int32_t x;
879
880 bufend = buf + (last - middle) - 1;
881 ss_blockswap(buf, middle, (saidx_t)(last - middle));
882
883 x = 0;
884 if (*bufend < 0)
885 {
886 p1 = PA + ~*bufend;
887 x |= 1;
888 }
889 else
890 {
891 p1 = PA + *bufend;
892 }
893 if (*(middle - 1) < 0)
894 {
895 p2 = PA + ~*(middle - 1);
896 x |= 2;
897 }
898 else
899 {
900 p2 = PA + *(middle - 1);
901 }
902 for (t = *(a = last - 1), b = bufend, c = middle - 1;;)
903 {
904 r = ss_compare(T, p1, p2, depth);
905 if (0 < r)
906 {
907 if (x & 1)
908 {
909 do
910 {
911 *a-- = *b, *b-- = *a;
912 }
913 while (*b < 0);
914 x ^= 1;
915 }
916 *a-- = *b;
917 if (b <= buf)
918 {
919 *buf = t;
920 break;
921 }
922 *b-- = *a;
923 if (*b < 0)
924 {
925 p1 = PA + ~*b;
926 x |= 1;
927 }
928 else
929 {
930 p1 = PA + *b;
931 }
932 }
933 else if (r < 0)
934 {
935 if (x & 2)
936 {
937 do
938 {
939 *a-- = *c, *c-- = *a;
940 }
941 while (*c < 0);
942 x ^= 2;
943 }
944 *a-- = *c, *c-- = *a;
945 if (c < first)
946 {
947 while (buf < b)
948 {
949 *a-- = *b, *b-- = *a;
950 }
951 *a = *b, *b = t;
952 break;
953 }
954 if (*c < 0)
955 {
956 p2 = PA + ~*c;
957 x |= 2;
958 }
959 else
960 {
961 p2 = PA + *c;
962 }
963 }
964 else
965 {
966 if (x & 1)
967 {
968 do
969 {
970 *a-- = *b, *b-- = *a;
971 }
972 while (*b < 0);
973 x ^= 1;
974 }
975 *a-- = ~*b;
976 if (b <= buf)
977 {
978 *buf = t;
979 break;
980 }
981 *b-- = *a;
982 if (x & 2)
983 {
984 do
985 {
986 *a-- = *c, *c-- = *a;
987 }
988 while (*c < 0);
989 x ^= 2;
990 }
991 *a-- = *c, *c-- = *a;
992 if (c < first)
993 {
994 while (buf < b)
995 {
996 *a-- = *b, *b-- = *a;
997 }
998 *a = *b, *b = t;
999 break;
1000 }
1001 if (*b < 0)
1002 {
1003 p1 = PA + ~*b;
1004 x |= 1;
1005 }
1006 else
1007 {
1008 p1 = PA + *b;
1009 }
1010 if (*c < 0)
1011 {
1012 p2 = PA + ~*c;
1013 x |= 2;
1014 }
1015 else
1016 {
1017 p2 = PA + *c;
1018 }
1019 }
1020 }
1021}
1022
1023/* D&C based merge. */
1024template <typename saidx_t>
1025inline void ss_swapmerge(uint8_t const * T,
1026 saidx_t const * PA,
1027 saidx_t * first,
1028 saidx_t * middle,
1029 saidx_t * last,
1030 saidx_t * buf,
1031 saidx_t bufsize,
1032 saidx_t depth)
1033{
1034# define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
1035# define MERGE_CHECK(a, b, c) \
1036 do \
1037 { \
1038 if (((c)&1) || (((c)&2) && (ss_compare(T, PA + GETIDX(*((a)-1)), PA + *(a), depth) == 0))) \
1039 { \
1040 *(a) = ~*(a); \
1041 } \
1042 if (((c)&4) && ((ss_compare(T, PA + GETIDX(*((b)-1)), PA + *(b), depth) == 0))) \
1043 { \
1044 *(b) = ~*(b); \
1045 } \
1046 } \
1047 while (0)
1048 struct
1049 {
1050 saidx_t *a, *b, *c;
1051 int32_t d;
1053 saidx_t *l, *r, *lm, *rm;
1054 saidx_t m, len, half;
1055 int32_t ssize;
1056 int32_t check, next;
1057
1058 for (check = 0, ssize = 0;;)
1059 {
1060 if ((last - middle) <= bufsize)
1061 {
1062 if ((first < middle) && (middle < last))
1063 {
1064 ss_mergebackward(T, PA, first, middle, last, buf, depth);
1065 }
1066 MERGE_CHECK(first, last, check);
1067 STACK_POP(first, middle, last, check);
1068 continue;
1069 }
1070
1071 if ((middle - first) <= bufsize)
1072 {
1073 if (first < middle)
1074 {
1075 ss_mergeforward(T, PA, first, middle, last, buf, depth);
1076 }
1077 MERGE_CHECK(first, last, check);
1078 STACK_POP(first, middle, last, check);
1079 continue;
1080 }
1081
1082 for (m = 0, len = std::min(middle - first, last - middle), half = len >> 1; 0 < len; len = half, half >>= 1)
1083 {
1084 if (ss_compare(T, PA + GETIDX(*(middle + m + half)), PA + GETIDX(*(middle - m - half - 1)), depth) < 0)
1085 {
1086 m += half + 1;
1087 half -= (len & 1) ^ 1;
1088 }
1089 }
1090
1091 if (0 < m)
1092 {
1093 lm = middle - m, rm = middle + m;
1094 ss_blockswap(lm, middle, m);
1095 l = r = middle, next = 0;
1096 if (rm < last)
1097 {
1098 if (*rm < 0)
1099 {
1100 *rm = ~*rm;
1101 if (first < lm)
1102 {
1103 for (; *--l < 0;)
1104 {}
1105 next |= 4;
1106 }
1107 next |= 1;
1108 }
1109 else if (first < lm)
1110 {
1111 for (; *r < 0; ++r)
1112 {}
1113 next |= 2;
1114 }
1115 }
1116
1117 if ((l - first) <= (last - r))
1118 {
1119 STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
1120 middle = lm, last = l, check = (check & 3) | (next & 4);
1121 }
1122 else
1123 {
1124 if ((next & 2) && (r == middle))
1125 {
1126 next ^= 6;
1127 }
1128 STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
1129 first = r, middle = rm, check = (next & 3) | (check & 4);
1130 }
1131 }
1132 else
1133 {
1134 if (ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0)
1135 {
1136 *middle = ~*middle;
1137 }
1138 MERGE_CHECK(first, last, check);
1139 STACK_POP(first, middle, last, check);
1140 }
1141 }
1142}
1143
1144#endif /* SS_BLOCKSIZE != 0 */
1145
1146/*- Function -*/
1147
1148/* Substring sort */
1149template <typename saidx_t>
1150void sssort(uint8_t const * T,
1151 saidx_t const * PA,
1152 saidx_t * first,
1153 saidx_t * last,
1154 saidx_t * buf,
1155 saidx_t bufsize,
1156 saidx_t depth,
1157 saidx_t n,
1158 int32_t lastsuffix)
1159{
1160 saidx_t * a;
1161#if SS_BLOCKSIZE != 0
1162 saidx_t *b, *middle, *curbuf;
1163 saidx_t j, k, curbufsize, limit;
1164#endif
1165 saidx_t i;
1166
1167 if (lastsuffix != 0)
1168 {
1169 ++first;
1170 }
1171
1172#if SS_BLOCKSIZE == 0
1173 ss_mintrosort(T, PA, first, last, depth);
1174#else
1175 if ((bufsize < SS_BLOCKSIZE) && (bufsize < (last - first)) && (bufsize < (limit = ss_isqrt(last - first))))
1176 {
1177 if (SS_BLOCKSIZE < limit)
1178 {
1179 limit = SS_BLOCKSIZE;
1180 }
1181 buf = middle = last - limit, bufsize = limit;
1182 }
1183 else
1184 {
1185 middle = last, limit = 0;
1186 }
1187 for (a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i)
1188 {
1189# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1190 ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth);
1191# elif 1 < SS_BLOCKSIZE
1192 ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth);
1193# endif
1194 curbufsize = last - (a + SS_BLOCKSIZE);
1195 curbuf = a + SS_BLOCKSIZE;
1196 if (curbufsize <= bufsize)
1197 {
1198 curbufsize = bufsize, curbuf = buf;
1199 }
1200 for (b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1)
1201 {
1202 ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
1203 }
1204 }
1205# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1206 ss_mintrosort(T, PA, a, middle, depth);
1207# elif 1 < SS_BLOCKSIZE
1208 ss_insertionsort(T, PA, a, middle, depth);
1209# endif
1210 for (k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1)
1211 {
1212 if (i & 1)
1213 {
1214 ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
1215 a -= k;
1216 }
1217 }
1218 if (limit != 0)
1219 {
1220# if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
1221 ss_mintrosort(T, PA, middle, last, depth);
1222# elif 1 < SS_BLOCKSIZE
1223 ss_insertionsort(T, PA, middle, last, depth);
1224# endif
1225 ss_inplacemerge(T, PA, first, middle, last, depth);
1226 }
1227#endif
1228
1229 if (lastsuffix != 0)
1230 {
1231 /* Insert last type B* suffix. */
1232 saidx_t PAi[2];
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)));
1235 ++a)
1236 {
1237 *(a - 1) = *a;
1238 }
1239 *(a - 1) = i;
1240 }
1241}
1242
1243inline int32_t tr_ilg(int32_t n)
1244{
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]);
1247}
1248
1249inline int32_t tr_ilg(int64_t n)
1250{
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]))
1253 : ((n & 0xffff0000)
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]));
1256}
1257
1258/* Simple insertionsort for small size groups. */
1259template <typename saidx_t>
1260inline void tr_insertionsort(saidx_t const * ISAd, saidx_t * first, saidx_t * last)
1261{
1262 saidx_t *a, *b;
1263 saidx_t t, r;
1264
1265 for (a = first + 1; a < last; ++a)
1266 {
1267 for (t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);)
1268 {
1269 do
1270 {
1271 *(b + 1) = *b;
1272 }
1273 while ((first <= --b) && (*b < 0));
1274 if (b < first)
1275 {
1276 break;
1277 }
1278 }
1279 if (r == 0)
1280 {
1281 *b = ~*b;
1282 }
1283 *(b + 1) = t;
1284 }
1285}
1286
1287template <typename saidx_t>
1288inline void tr_fixdown(saidx_t const * ISAd, saidx_t * SA, saidx_t i, saidx_t size)
1289{
1290 saidx_t j, k;
1291 saidx_t v;
1292 saidx_t c, d, e;
1293
1294 for (v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k)
1295 {
1296 d = ISAd[SA[k = j++]];
1297 if (d < (e = ISAd[SA[j]]))
1298 {
1299 k = j;
1300 d = e;
1301 }
1302 if (d <= c)
1303 {
1304 break;
1305 }
1306 }
1307 SA[i] = v;
1308}
1309
1310/* Simple top-down heapsort. */
1311template <typename saidx_t>
1312inline void tr_heapsort(saidx_t const * ISAd, saidx_t * SA, saidx_t size)
1313{
1314 saidx_t i, m;
1315 saidx_t t;
1316
1317 m = size;
1318 if ((size % 2) == 0)
1319 {
1320 m--;
1321 if (ISAd[SA[m / 2]] < ISAd[SA[m]])
1322 {
1323 std::swap(SA[m], SA[m / 2]);
1324 }
1325 }
1326
1327 for (i = m / 2 - 1; 0 <= i; --i)
1328 {
1329 tr_fixdown(ISAd, SA, i, m);
1330 }
1331 if ((size % 2) == 0)
1332 {
1333 std::swap(SA[0], SA[m]);
1334 tr_fixdown(ISAd, SA, (saidx_t)0, m);
1335 }
1336 for (i = m - 1; 0 < i; --i)
1337 {
1338 t = SA[0], SA[0] = SA[i];
1339 tr_fixdown(ISAd, SA, (saidx_t)0, i);
1340 SA[i] = t;
1341 }
1342}
1343
1344/* Returns the median of three elements. */
1345template <typename saidx_t>
1346inline saidx_t * tr_median3(saidx_t const * ISAd, saidx_t * v1, saidx_t * v2, saidx_t * v3)
1347{
1348 if (ISAd[*v1] > ISAd[*v2])
1349 {
1350 std::swap(v1, v2);
1351 }
1352 if (ISAd[*v2] > ISAd[*v3])
1353 {
1354 if (ISAd[*v1] > ISAd[*v3])
1355 {
1356 return v1;
1357 }
1358 else
1359 {
1360 return v3;
1361 }
1362 }
1363 return v2;
1364}
1365
1366/* Returns the median of five elements. */
1367template <typename saidx_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)
1369{
1370 if (ISAd[*v2] > ISAd[*v3])
1371 {
1372 std::swap(v2, v3);
1373 }
1374 if (ISAd[*v4] > ISAd[*v5])
1375 {
1376 std::swap(v4, v5);
1377 }
1378 if (ISAd[*v2] > ISAd[*v4])
1379 {
1380 std::swap(v2, v4);
1381 std::swap(v3, v5);
1382 }
1383 if (ISAd[*v1] > ISAd[*v3])
1384 {
1385 std::swap(v1, v3);
1386 }
1387 if (ISAd[*v1] > ISAd[*v4])
1388 {
1389 std::swap(v1, v4);
1390 std::swap(v3, v5);
1391 }
1392 if (ISAd[*v3] > ISAd[*v4])
1393 {
1394 return v4;
1395 }
1396 return v3;
1397}
1398
1399/* Returns the pivot element. */
1400template <typename saidx_t>
1401inline saidx_t * tr_pivot(saidx_t const * ISAd, saidx_t * first, saidx_t * last)
1402{
1403 saidx_t * middle;
1404 saidx_t t;
1405
1406 t = last - first;
1407 middle = first + t / 2;
1408
1409 if (t <= 512)
1410 {
1411 if (t <= 32)
1412 {
1413 return tr_median3(ISAd, first, middle, last - 1);
1414 }
1415 else
1416 {
1417 t >>= 2;
1418 return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1419 }
1420 }
1421 t >>= 3;
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);
1426}
1427
1428template <typename saidx_t>
1430{
1431 saidx_t chance;
1432 saidx_t remain;
1433 saidx_t incval;
1434 saidx_t count;
1435};
1436
1437template <typename saidx_t>
1438using trbudget_t = struct _trbudget_t<saidx_t>;
1439// typedef struct _trbudget_t trbudget_t;
1440
1441template <typename saidx_t>
1442inline void trbudget_init(trbudget_t<saidx_t> * budget, saidx_t chance, saidx_t incval)
1443{
1444 budget->chance = chance;
1445 budget->remain = budget->incval = incval;
1446}
1447
1448template <typename saidx_t>
1449inline int32_t trbudget_check(trbudget_t<saidx_t> * budget, saidx_t size)
1450{
1451 if (size <= budget->remain)
1452 {
1453 budget->remain -= size;
1454 return 1;
1455 }
1456 if (budget->chance == 0)
1457 {
1458 budget->count += size;
1459 return 0;
1460 }
1461 budget->remain += budget->incval - size;
1462 budget->chance -= 1;
1463 return 1;
1464}
1465
1466template <typename saidx_t>
1467inline void tr_partition(saidx_t const * ISAd,
1468 saidx_t * first,
1469 saidx_t * middle,
1470 saidx_t * last,
1471 saidx_t ** pa,
1472 saidx_t ** pb,
1473 saidx_t v)
1474{
1475 saidx_t *a, *b, *c, *d, *e, *f;
1476 saidx_t t, s;
1477 saidx_t x = 0;
1478
1479 for (b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);)
1480 {}
1481 if (((a = b) < last) && (x < v))
1482 {
1483 for (; (++b < last) && ((x = ISAd[*b]) <= v);)
1484 {
1485 if (x == v)
1486 {
1487 std::swap(*b, *a);
1488 ++a;
1489 }
1490 }
1491 }
1492 for (c = last; (b < --c) && ((x = ISAd[*c]) == v);)
1493 {}
1494 if ((b < (d = c)) && (x > v))
1495 {
1496 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1497 {
1498 if (x == v)
1499 {
1500 std::swap(*c, *d);
1501 --d;
1502 }
1503 }
1504 }
1505 for (; b < c;)
1506 {
1507 std::swap(*b, *c);
1508 for (; (++b < c) && ((x = ISAd[*b]) <= v);)
1509 {
1510 if (x == v)
1511 {
1512 std::swap(*b, *a);
1513 ++a;
1514 }
1515 }
1516 for (; (b < --c) && ((x = ISAd[*c]) >= v);)
1517 {
1518 if (x == v)
1519 {
1520 std::swap(*c, *d);
1521 --d;
1522 }
1523 }
1524 }
1525
1526 if (a <= d)
1527 {
1528 c = b - 1;
1529 if ((s = a - first) > (t = b - a))
1530 {
1531 s = t;
1532 }
1533 for (e = first, f = b - s; 0 < s; --s, ++e, ++f)
1534 {
1535 std::swap(*e, *f);
1536 }
1537 if ((s = d - c) > (t = last - d - 1))
1538 {
1539 s = t;
1540 }
1541 for (e = b, f = last - s; 0 < s; --s, ++e, ++f)
1542 {
1543 std::swap(*e, *f);
1544 }
1545 first += (b - a), last -= (d - c);
1546 }
1547 *pa = first, *pb = last;
1548}
1549
1550template <typename saidx_t>
1551inline void
1552tr_copy(saidx_t * ISA, saidx_t const * SA, saidx_t * first, saidx_t * a, saidx_t * b, saidx_t * last, saidx_t depth)
1553{
1554 /* sort suffixes of middle partition
1555 by using sorted order of suffixes of left and right partition. */
1556 saidx_t *c, *d, *e;
1557 saidx_t s, v;
1558
1559 v = b - SA - 1;
1560 for (c = first, d = a - 1; c <= d; ++c)
1561 {
1562 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1563 {
1564 *++d = s;
1565 ISA[s] = d - SA;
1566 }
1567 }
1568 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1569 {
1570 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1571 {
1572 *--d = s;
1573 ISA[s] = d - SA;
1574 }
1575 }
1576}
1577
1578template <typename saidx_t>
1579inline void tr_partialcopy(saidx_t * ISA,
1580 saidx_t const * SA,
1581 saidx_t * first,
1582 saidx_t * a,
1583 saidx_t * b,
1584 saidx_t * last,
1585 saidx_t depth)
1586{
1587 saidx_t *c, *d, *e;
1588 saidx_t s, v;
1589 saidx_t rank, lastrank, newrank = -1;
1590
1591 v = b - SA - 1;
1592 lastrank = -1;
1593 for (c = first, d = a - 1; c <= d; ++c)
1594 {
1595 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1596 {
1597 *++d = s;
1598 rank = ISA[s + depth];
1599 if (lastrank != rank)
1600 {
1601 lastrank = rank;
1602 newrank = d - SA;
1603 }
1604 ISA[s] = newrank;
1605 }
1606 }
1607
1608 lastrank = -1;
1609 for (e = d; first <= e; --e)
1610 {
1611 rank = ISA[*e];
1612 if (lastrank != rank)
1613 {
1614 lastrank = rank;
1615 newrank = e - SA;
1616 }
1617 if (newrank != rank)
1618 {
1619 ISA[*e] = newrank;
1620 }
1621 }
1622
1623 lastrank = -1;
1624 for (c = last - 1, e = d + 1, d = b; e < d; --c)
1625 {
1626 if ((0 <= (s = *c - depth)) && (ISA[s] == v))
1627 {
1628 *--d = s;
1629 rank = ISA[s + depth];
1630 if (lastrank != rank)
1631 {
1632 lastrank = rank;
1633 newrank = d - SA;
1634 }
1635 ISA[s] = newrank;
1636 }
1637 }
1638}
1639
1640template <typename saidx_t>
1641inline void tr_introsort(saidx_t * ISA,
1642 saidx_t const * ISAd,
1643 saidx_t * SA,
1644 saidx_t * first,
1645 saidx_t * last,
1646 trbudget_t<saidx_t> * budget)
1647{
1648 struct
1649 {
1650 saidx_t const * a;
1651 saidx_t *b, *c;
1652 int32_t d, e;
1654 saidx_t *a, *b, *c;
1655 saidx_t v, x = 0;
1656 saidx_t incr = ISAd - ISA;
1657 int32_t limit, next;
1658 int32_t ssize, trlink = -1;
1659
1660 for (ssize = 0, limit = tr_ilg((saidx_t)(last - first));;)
1661 {
1662
1663 if (limit < 0)
1664 {
1665 if (limit == -1)
1666 {
1667 /* tandem repeat partition */
1668 tr_partition(ISAd - incr, first, first, last, &a, &b, (saidx_t)(last - SA - 1));
1669
1670 /* update ranks */
1671 if (a < last)
1672 {
1673 for (c = first, v = a - SA - 1; c < a; ++c)
1674 {
1675 ISA[*c] = v;
1676 }
1677 }
1678 if (b < last)
1679 {
1680 for (c = a, v = b - SA - 1; c < b; ++c)
1681 {
1682 ISA[*c] = v;
1683 }
1684 }
1685
1686 /* push */
1687 if (1 < (b - a))
1688 {
1689 STACK_PUSH5(NULL, a, b, 0, 0);
1690 STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1691 trlink = ssize - 2;
1692 }
1693 if ((a - first) <= (last - b))
1694 {
1695 if (1 < (a - first))
1696 {
1697 STACK_PUSH5(ISAd, b, last, tr_ilg((saidx_t)(last - b)), trlink);
1698 last = a, limit = tr_ilg((saidx_t)(a - first));
1699 }
1700 else if (1 < (last - b))
1701 {
1702 first = b, limit = tr_ilg((saidx_t)(last - b));
1703 }
1704 else
1705 {
1706 STACK_POP5(ISAd, first, last, limit, trlink);
1707 }
1708 }
1709 else
1710 {
1711 if (1 < (last - b))
1712 {
1713 STACK_PUSH5(ISAd, first, a, tr_ilg((saidx_t)(a - first)), trlink);
1714 first = b, limit = tr_ilg((saidx_t)(last - b));
1715 }
1716 else if (1 < (a - first))
1717 {
1718 last = a, limit = tr_ilg((saidx_t)(a - first));
1719 }
1720 else
1721 {
1722 STACK_POP5(ISAd, first, last, limit, trlink);
1723 }
1724 }
1725 }
1726 else if (limit == -2)
1727 {
1728 /* tandem repeat copy */
1729 a = stack[--ssize].b, b = stack[ssize].c;
1730 if (stack[ssize].d == 0)
1731 {
1732 tr_copy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA));
1733 }
1734 else
1735 {
1736 if (0 <= trlink)
1737 {
1738 stack[trlink].d = -1;
1739 }
1740 tr_partialcopy(ISA, SA, first, a, b, last, (saidx_t)(ISAd - ISA));
1741 }
1742 STACK_POP5(ISAd, first, last, limit, trlink);
1743 }
1744 else
1745 {
1746 /* sorted partition */
1747 if (0 <= *first)
1748 {
1749 a = first;
1750 do
1751 {
1752 ISA[*a] = a - SA;
1753 }
1754 while ((++a < last) && (0 <= *a));
1755 first = a;
1756 }
1757 if (first < last)
1758 {
1759 a = first;
1760 do
1761 {
1762 *a = ~*a;
1763 }
1764 while (*++a < 0);
1765 next = (ISA[*a] != ISAd[*a]) ? tr_ilg((saidx_t)(a - first + 1)) : -1;
1766 if (++a < last)
1767 {
1768 for (b = first, v = a - SA - 1; b < a; ++b)
1769 {
1770 ISA[*b] = v;
1771 }
1772 }
1773
1774 /* push */
1775 if (trbudget_check(budget, (saidx_t)(a - first)))
1776 {
1777 if ((a - first) <= (last - a))
1778 {
1779 STACK_PUSH5(ISAd, a, last, -3, trlink);
1780 ISAd += incr, last = a, limit = next;
1781 }
1782 else
1783 {
1784 if (1 < (last - a))
1785 {
1786 STACK_PUSH5(ISAd + incr, first, a, next, trlink);
1787 first = a, limit = -3;
1788 }
1789 else
1790 {
1791 ISAd += incr, last = a, limit = next;
1792 }
1793 }
1794 }
1795 else
1796 {
1797 if (0 <= trlink)
1798 {
1799 stack[trlink].d = -1;
1800 }
1801 if (1 < (last - a))
1802 {
1803 first = a, limit = -3;
1804 }
1805 else
1806 {
1807 STACK_POP5(ISAd, first, last, limit, trlink);
1808 }
1809 }
1810 }
1811 else
1812 {
1813 STACK_POP5(ISAd, first, last, limit, trlink);
1814 }
1815 }
1816 continue;
1817 }
1818
1819 if ((last - first) <= TR_INSERTIONSORT_THRESHOLD)
1820 {
1821 tr_insertionsort(ISAd, first, last);
1822 limit = -3;
1823 continue;
1824 }
1825
1826 if (limit-- == 0)
1827 {
1828 tr_heapsort(ISAd, first, (saidx_t)(last - first));
1829 for (a = last - 1; first < a; a = b)
1830 {
1831 for (x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b)
1832 {
1833 *b = ~*b;
1834 }
1835 }
1836 limit = -3;
1837 continue;
1838 }
1839
1840 /* choose pivot */
1841 a = tr_pivot(ISAd, first, last);
1842 std::swap(*first, *a);
1843 v = ISAd[*first];
1844
1845 /* partition */
1846 tr_partition(ISAd, first, first + 1, last, &a, &b, v);
1847 if ((last - first) != (b - a))
1848 {
1849 next = (ISA[*a] != v) ? tr_ilg((saidx_t)(b - a)) : -1;
1850
1851 /* update ranks */
1852 for (c = first, v = a - SA - 1; c < a; ++c)
1853 {
1854 ISA[*c] = v;
1855 }
1856 if (b < last)
1857 {
1858 for (c = a, v = b - SA - 1; c < b; ++c)
1859 {
1860 ISA[*c] = v;
1861 }
1862 }
1863
1864 /* push */
1865 if ((1 < (b - a)) && (trbudget_check(budget, (saidx_t)(b - a))))
1866 {
1867 if ((a - first) <= (last - b))
1868 {
1869 if ((last - b) <= (b - a))
1870 {
1871 if (1 < (a - first))
1872 {
1873 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1874 STACK_PUSH5(ISAd, b, last, limit, trlink);
1875 last = a;
1876 }
1877 else if (1 < (last - b))
1878 {
1879 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1880 first = b;
1881 }
1882 else
1883 {
1884 ISAd += incr, first = a, last = b, limit = next;
1885 }
1886 }
1887 else if ((a - first) <= (b - a))
1888 {
1889 if (1 < (a - first))
1890 {
1891 STACK_PUSH5(ISAd, b, last, limit, trlink);
1892 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1893 last = a;
1894 }
1895 else
1896 {
1897 STACK_PUSH5(ISAd, b, last, limit, trlink);
1898 ISAd += incr, first = a, last = b, limit = next;
1899 }
1900 }
1901 else
1902 {
1903 STACK_PUSH5(ISAd, b, last, limit, trlink);
1904 STACK_PUSH5(ISAd, first, a, limit, trlink);
1905 ISAd += incr, first = a, last = b, limit = next;
1906 }
1907 }
1908 else
1909 {
1910 if ((a - first) <= (b - a))
1911 {
1912 if (1 < (last - b))
1913 {
1914 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1915 STACK_PUSH5(ISAd, first, a, limit, trlink);
1916 first = b;
1917 }
1918 else if (1 < (a - first))
1919 {
1920 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1921 last = a;
1922 }
1923 else
1924 {
1925 ISAd += incr, first = a, last = b, limit = next;
1926 }
1927 }
1928 else if ((last - b) <= (b - a))
1929 {
1930 if (1 < (last - b))
1931 {
1932 STACK_PUSH5(ISAd, first, a, limit, trlink);
1933 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1934 first = b;
1935 }
1936 else
1937 {
1938 STACK_PUSH5(ISAd, first, a, limit, trlink);
1939 ISAd += incr, first = a, last = b, limit = next;
1940 }
1941 }
1942 else
1943 {
1944 STACK_PUSH5(ISAd, first, a, limit, trlink);
1945 STACK_PUSH5(ISAd, b, last, limit, trlink);
1946 ISAd += incr, first = a, last = b, limit = next;
1947 }
1948 }
1949 }
1950 else
1951 {
1952 if ((1 < (b - a)) && (0 <= trlink))
1953 {
1954 stack[trlink].d = -1;
1955 }
1956 if ((a - first) <= (last - b))
1957 {
1958 if (1 < (a - first))
1959 {
1960 STACK_PUSH5(ISAd, b, last, limit, trlink);
1961 last = a;
1962 }
1963 else if (1 < (last - b))
1964 {
1965 first = b;
1966 }
1967 else
1968 {
1969 STACK_POP5(ISAd, first, last, limit, trlink);
1970 }
1971 }
1972 else
1973 {
1974 if (1 < (last - b))
1975 {
1976 STACK_PUSH5(ISAd, first, a, limit, trlink);
1977 first = b;
1978 }
1979 else if (1 < (a - first))
1980 {
1981 last = a;
1982 }
1983 else
1984 {
1985 STACK_POP5(ISAd, first, last, limit, trlink);
1986 }
1987 }
1988 }
1989 }
1990 else
1991 {
1992 if (trbudget_check(budget, (saidx_t)(last - first)))
1993 {
1994 limit = tr_ilg((saidx_t)(last - first)), ISAd += incr;
1995 }
1996 else
1997 {
1998 if (0 <= trlink)
1999 {
2000 stack[trlink].d = -1;
2001 }
2002 STACK_POP5(ISAd, first, last, limit, trlink);
2003 }
2004 }
2005 }
2006}
2007
2008/*- Function -*/
2009
2010/* Tandem repeat sort */
2011template <typename saidx_t>
2012inline void trsort(saidx_t * ISA, saidx_t * SA, saidx_t n, saidx_t depth)
2013{
2014 saidx_t * ISAd;
2015 saidx_t *first, *last;
2016 trbudget_t<saidx_t> budget;
2017 saidx_t t, skip, unsorted;
2018
2019 trbudget_init(&budget, (saidx_t)(tr_ilg(n) * 2 / 3), n);
2020 /* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */
2021 for (ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA)
2022 {
2023 first = SA;
2024 skip = 0;
2025 unsorted = 0;
2026 do
2027 {
2028 if ((t = *first) < 0)
2029 {
2030 first -= t;
2031 skip += t;
2032 }
2033 else
2034 {
2035 if (skip != 0)
2036 {
2037 *(first + skip) = skip;
2038 skip = 0;
2039 }
2040 last = SA + ISA[t] + 1;
2041 if (1 < (last - first))
2042 {
2043 budget.count = 0;
2044 tr_introsort(ISA, ISAd, SA, first, last, &budget);
2045 if (budget.count != 0)
2046 {
2047 unsorted += budget.count;
2048 }
2049 else
2050 {
2051 skip = first - last;
2052 }
2053 }
2054 else if ((last - first) == 1)
2055 {
2056 skip = -1;
2057 }
2058 first = last;
2059 }
2060 }
2061 while (first < (SA + n));
2062 if (skip != 0)
2063 {
2064 *(first + skip) = skip;
2065 }
2066 if (unsorted == 0)
2067 {
2068 break;
2069 }
2070 }
2071}
2072
2073/* Sorts suffixes of type B*. */
2074template <typename saidx_t>
2075inline saidx_t sort_typeBstar(uint8_t const * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n)
2076{
2077 saidx_t *PAb, *ISAb, *buf;
2078#ifdef _OPENMP
2079 saidx_t * curbuf;
2080 saidx_t l;
2081#endif
2082 saidx_t i, j, k, t, m, bufsize;
2083 int32_t c0, c1;
2084#ifdef _OPENMP
2085 int32_t d0, d1;
2086 int tmp;
2087#endif
2088
2089 /* Initialize bucket arrays. */
2090 for (i = 0; i < BUCKET_A_SIZE; ++i)
2091 {
2092 bucket_A[i] = 0;
2093 }
2094 for (i = 0; i < BUCKET_B_SIZE; ++i)
2095 {
2096 bucket_B[i] = 0;
2097 }
2098
2099 /* Count the number of occurrences of the first one or two characters of each
2100 * type A, B and B* suffix. Moreover, store the beginning position of all
2101 type B* suffixes into the array SA. */
2102 for (i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;)
2103 {
2104 /* type A suffix. */
2105 do
2106 {
2107 ++BUCKET_A(c1 = c0);
2108 }
2109 while ((0 <= --i) && ((c0 = T[i]) >= c1));
2110 if (0 <= i)
2111 {
2112 /* type B* suffix. */
2113 ++BUCKET_BSTAR(c0, c1);
2114 SA[--m] = i;
2115 /* type B suffix. */
2116 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0)
2117 {
2118 ++BUCKET_B(c0, c1);
2119 }
2120 }
2121 }
2122 m = n - m;
2123 /*
2124 * note:
2125 * A type B* suffix is lexicographically smaller than a type B suffix that
2126 * begins with the same first two characters.
2127 */
2128
2129 /* Calculate the index of start/end point of each bucket. */
2130 for (c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0)
2131 {
2132 t = i + BUCKET_A(c0);
2133 BUCKET_A(c0) = i + j; /* start point */
2134 i = t + BUCKET_B(c0, c0);
2135 for (c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1)
2136 {
2137 j += BUCKET_BSTAR(c0, c1);
2138 BUCKET_BSTAR(c0, c1) = j; /* end point */
2139 i += BUCKET_B(c0, c1);
2140 }
2141 }
2142
2143 if (0 < m)
2144 {
2145 /* Sort the type B* suffixes by their first two characters. */
2146 PAb = SA + n - m;
2147 ISAb = SA + m;
2148 for (i = m - 2; 0 <= i; --i)
2149 {
2150 t = PAb[i], c0 = T[t], c1 = T[t + 1];
2151 SA[--BUCKET_BSTAR(c0, c1)] = i;
2152 }
2153 t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
2154 SA[--BUCKET_BSTAR(c0, c1)] = m - 1;
2155
2156 /* Sort the type B* substrings using sssort. */
2157#ifdef _OPENMP
2158 tmp = omp_get_max_threads();
2159 buf = SA + m, bufsize = (n - (2 * m)) / tmp;
2160 c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m;
2161# pragma omp parallel default(shared) private(curbuf, k, l, d0, d1, tmp)
2162 {
2163 tmp = omp_get_thread_num();
2164 curbuf = buf + tmp * bufsize;
2165 k = 0;
2166 for (;;)
2167 {
2168# pragma omp critical(sssort_lock)
2169 {
2170 if (0 < (l = j))
2171 {
2172 d0 = c0, d1 = c1;
2173 do
2174 {
2175 k = BUCKET_BSTAR(d0, d1);
2176 if (--d1 <= d0)
2177 {
2178 d1 = ALPHABET_SIZE - 1;
2179 if (--d0 < 0)
2180 {
2181 break;
2182 }
2183 }
2184 }
2185 while (((l - k) <= 1) && (0 < (l = k)));
2186 c0 = d0, c1 = d1, j = k;
2187 }
2188 }
2189 if (l == 0)
2190 {
2191 break;
2192 }
2193 sssort(T, PAb, SA + k, SA + l, curbuf, bufsize, (saidx_t)2, n, *(SA + k) == (m - 1));
2194 }
2195 }
2196#else
2197 buf = SA + m, bufsize = n - (2 * m);
2198 for (c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0)
2199 {
2200 for (c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1)
2201 {
2202 i = BUCKET_BSTAR(c0, c1);
2203 if (1 < (j - i))
2204 {
2205 sssort(T, PAb, SA + i, SA + j, buf, bufsize, (saidx_t)2, n, *(SA + i) == (m - 1));
2206 }
2207 }
2208 }
2209#endif
2210
2211 /* Compute ranks of type B* substrings. */
2212 for (i = m - 1; 0 <= i; --i)
2213 {
2214 if (0 <= SA[i])
2215 {
2216 j = i;
2217 do
2218 {
2219 ISAb[SA[i]] = i;
2220 }
2221 while ((0 <= --i) && (0 <= SA[i]));
2222 SA[i + 1] = i - j;
2223 if (i <= 0)
2224 {
2225 break;
2226 }
2227 }
2228 j = i;
2229 do
2230 {
2231 ISAb[SA[i] = ~SA[i]] = j;
2232 }
2233 while (SA[--i] < 0);
2234 ISAb[SA[i]] = j;
2235 }
2236
2237 /* Construct the inverse suffix array of type B* suffixes using trsort. */
2238 trsort(ISAb, SA, m, (saidx_t)1);
2239
2240 /* Set the sorted order of tyoe B* suffixes. */
2241 for (i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;)
2242 {
2243 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0)
2244 {}
2245 if (0 <= i)
2246 {
2247 t = i;
2248 for (--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0)
2249 {}
2250 SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
2251 }
2252 }
2253
2254 /* Calculate the index of start/end point of each bucket. */
2255 BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */
2256 for (c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0)
2257 {
2258 i = BUCKET_A(c0 + 1) - 1;
2259 for (c1 = ALPHABET_SIZE - 1; c0 < c1; --c1)
2260 {
2261 t = i - BUCKET_B(c0, c1);
2262 BUCKET_B(c0, c1) = i; /* end point */
2263
2264 /* Move all type B* suffixes to the correct position. */
2265 for (i = t, j = BUCKET_BSTAR(c0, c1); j <= k; --i, --k)
2266 {
2267 SA[i] = SA[k];
2268 }
2269 }
2270 BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */
2271 BUCKET_B(c0, c0) = i; /* end point */
2272 }
2273 }
2274
2275 return m;
2276}
2277
2278/* Constructs the suffix array by using the sorted order of type B* suffixes. */
2279template <typename saidx_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)
2281{
2282 saidx_t *i, *j, *k;
2283 saidx_t s;
2284 int32_t c0, c1, c2;
2285
2286 if (0 < m)
2287 {
2288 /* Construct the sorted order of type B suffixes by using
2289 the sorted order of type B* suffixes. */
2290 for (c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1)
2291 {
2292 /* Scan the suffix array from right to left. */
2293 for (i = SA + BUCKET_BSTAR(c1, c1 + 1), j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
2294 {
2295 if (0 < (s = *j))
2296 {
2297 assert(T[s] == c1);
2298 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
2299 assert(T[s - 1] <= T[s]);
2300 *j = ~s;
2301 c0 = T[--s];
2302 if ((0 < s) && (T[s - 1] > c0))
2303 {
2304 s = ~s;
2305 }
2306 if (c0 != c2)
2307 {
2308 if (0 <= c2)
2309 {
2310 BUCKET_B(c2, c1) = k - SA;
2311 }
2312 k = SA + BUCKET_B(c2 = c0, c1);
2313 }
2314 assert(k < j);
2315 *k-- = s;
2316 }
2317 else
2318 {
2319 assert(((s == 0) && (T[s] == c1)) || (s < 0));
2320 *j = ~s;
2321 }
2322 }
2323 }
2324 }
2325
2326 /* Construct the suffix array by using
2327 the sorted order of type B suffixes. */
2328 k = SA + BUCKET_A(c2 = T[n - 1]);
2329 *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
2330 /* Scan the suffix array from left to right. */
2331 for (i = SA, j = SA + n; i < j; ++i)
2332 {
2333 if (0 < (s = *i))
2334 {
2335 assert(T[s - 1] >= T[s]);
2336 c0 = T[--s];
2337 if ((s == 0) || (T[s - 1] < c0))
2338 {
2339 s = ~s;
2340 }
2341 if (c0 != c2)
2342 {
2343 BUCKET_A(c2) = k - SA;
2344 k = SA + BUCKET_A(c2 = c0);
2345 }
2346 assert(i < k);
2347 *k++ = s;
2348 }
2349 else
2350 {
2351 assert(s < 0);
2352 *i = ~s;
2353 }
2354 }
2355}
2356
2357/* Constructs the burrows-wheeler transformed string directly
2358 by using the sorted order of type B* suffixes. */
2359template <typename saidx_t>
2360inline saidx_t
2361construct_BWT(uint8_t const * T, saidx_t * SA, saidx_t * bucket_A, saidx_t * bucket_B, saidx_t n, saidx_t m)
2362{
2363 saidx_t *i, *j, *k, *orig;
2364 saidx_t s;
2365 int32_t c0, c1, c2;
2366
2367 if (0 < m)
2368 {
2369 /* Construct the sorted order of type B suffixes by using
2370 the sorted order of type B* suffixes. */
2371 for (c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1)
2372 {
2373 /* Scan the suffix array from right to left. */
2374 for (i = SA + BUCKET_BSTAR(c1, c1 + 1), j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; i <= j; --j)
2375 {
2376 if (0 < (s = *j))
2377 {
2378 assert(T[s] == c1);
2379 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
2380 assert(T[s - 1] <= T[s]);
2381 c0 = T[--s];
2382 *j = ~((saidx_t)c0);
2383 if ((0 < s) && (T[s - 1] > c0))
2384 {
2385 s = ~s;
2386 }
2387 if (c0 != c2)
2388 {
2389 if (0 <= c2)
2390 {
2391 BUCKET_B(c2, c1) = k - SA;
2392 }
2393 k = SA + BUCKET_B(c2 = c0, c1);
2394 }
2395 assert(k < j);
2396 *k-- = s;
2397 }
2398 else if (s != 0)
2399 {
2400 *j = ~s;
2401#ifndef NDEBUG
2402 }
2403 else
2404 {
2405 assert(T[s] == c1);
2406#endif
2407 }
2408 }
2409 }
2410 }
2411
2412 /* Construct the BWTed string by using
2413 the sorted order of type B suffixes. */
2414 k = SA + BUCKET_A(c2 = T[n - 1]);
2415 *k++ = (T[n - 2] < c2) ? ~((saidx_t)T[n - 2]) : (n - 1);
2416 /* Scan the suffix array from left to right. */
2417 for (i = SA, j = SA + n, orig = SA; i < j; ++i)
2418 {
2419 if (0 < (s = *i))
2420 {
2421 assert(T[s - 1] >= T[s]);
2422 c0 = T[--s];
2423 *i = c0;
2424 if ((0 < s) && (T[s - 1] < c0))
2425 {
2426 s = ~((saidx_t)T[s - 1]);
2427 }
2428 if (c0 != c2)
2429 {
2430 BUCKET_A(c2) = k - SA;
2431 k = SA + BUCKET_A(c2 = c0);
2432 }
2433 assert(i < k);
2434 *k++ = s;
2435 }
2436 else if (s != 0)
2437 {
2438 *i = ~s;
2439 }
2440 else
2441 {
2442 orig = i;
2443 }
2444 }
2445
2446 return orig - SA;
2447}
2448
2449/*- Function -*/
2450
2451template <typename saidx_t>
2452int32_t divsufsort(uint8_t const * T, saidx_t * SA, saidx_t n)
2453{
2454 saidx_t *bucket_A, *bucket_B;
2455 saidx_t m;
2456 int32_t err = 0;
2457
2458 /* Check arguments. */
2459 if ((T == NULL) || (SA == NULL) || (n < 0))
2460 {
2461 return -1;
2462 }
2463 else if (n == 0)
2464 {
2465 return 0;
2466 }
2467 else if (n == 1)
2468 {
2469 SA[0] = 0;
2470 return 0;
2471 }
2472 else if (n == 2)
2473 {
2474 m = (T[0] < T[1]);
2475 SA[m ^ 1] = 0, SA[m] = 1;
2476 return 0;
2477 }
2478
2479 bucket_A = (saidx_t *)malloc(BUCKET_A_SIZE * sizeof(saidx_t));
2480 bucket_B = (saidx_t *)malloc(BUCKET_B_SIZE * sizeof(saidx_t));
2481
2482 /* Suffixsort. */
2483 if ((bucket_A != NULL) && (bucket_B != NULL))
2484 {
2485 m = sort_typeBstar(T, SA, bucket_A, bucket_B, n);
2486 construct_SA(T, SA, bucket_A, bucket_B, n, m);
2487 }
2488 else
2489 {
2490 err = -2;
2491 }
2492
2493 free(bucket_B);
2494 free(bucket_A);
2495
2496 return err;
2497}
2498
2499// template <typename saidx_t>
2500// saidx_t
2501// divbwt(const uint8_t *T, uint8_t *U, saidx_t *A, saidx_t n) {
2502// saidx_t *B;
2503// saidx_t *bucket_A, *bucket_B;
2504// saidx_t m, pidx, i;
2505//
2506// /* Check arguments. */
2507// if((T == NULL) || (U == NULL) || (n < 0)) { return -1; }
2508// else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; }
2509//
2510// if((B = A) == NULL) { B = (saidx_t *)malloc((size_t)(n + 1) * sizeof(saidx_t)); }
2511// bucket_A = (saidx_t *)malloc(BUCKET_A_SIZE * sizeof(saidx_t));
2512// bucket_B = (saidx_t *)malloc(BUCKET_B_SIZE * sizeof(saidx_t));
2513//
2514// /* Burrows-Wheeler Transform. */
2515// if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) {
2516// m = sort_typeBstar(T, B, bucket_A, bucket_B, n);
2517// pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m);
2518//
2519// /* Copy to output string. */
2520// U[0] = T[n - 1];
2521// for(i = 0; i < pidx; ++i) { U[i + 1] = (uint8_t)B[i]; }
2522// for(i += 1; i < n; ++i) { U[i] = (uint8_t)B[i]; }
2523// pidx += 1;
2524// } else {
2525// pidx = -2;
2526// }
2527//
2528// free(bucket_B);
2529// free(bucket_A);
2530// if(A == NULL) { free(B); }
2531//
2532// return pidx;
2533// }
2534
2535inline int32_t divsufsort64(uint8_t const * T, int64_t * SA, int64_t n)
2536{
2537 return divsufsort(T, SA, n);
2538}
2539
2540template <typename saidx_t>
2541inline int _compare(uint8_t const * T, saidx_t Tsize, uint8_t const * P, saidx_t Psize, saidx_t suf, saidx_t * match)
2542{
2543 saidx_t i, j;
2544 int32_t r;
2545 for (i = suf + *match, j = *match, r = 0; (i < Tsize) && (j < Psize) && ((r = T[i] - P[j]) == 0); ++i, ++j)
2546 {}
2547 *match = j;
2548 return (r == 0) ? -(j != Psize) : r;
2549}
2550
2551/* Burrows-Wheeler transform. */
2552// TODO: use?
2553// template <typename saidx_t>
2554// int32_t
2555// bw_transform(const uint8_t *T, uint8_t *U, saidx_t *SA,
2556// saidx_t n, saidx_t *idx) {
2557// saidx_t *A, i, j, p, t;
2558// int32_t c;
2559//
2560// /* Check arguments. */
2561// if((T == NULL) || (U == NULL) || (n < 0) || (idx == NULL)) { return -1; }
2562// if(n <= 1) {
2563// if(n == 1) { U[0] = T[0]; }
2564// *idx = n;
2565// return 0;
2566// }
2567//
2568// if((A = SA) == NULL) {
2569// i = divbwt(T, U, NULL, n);
2570// if(0 <= i) { *idx = i; i = 0; }
2571// return (int32_t)i;
2572// }
2573//
2574// /* BW transform. */
2575// if(T == U) {
2576// t = n;
2577// for(i = 0, j = 0; i < n; ++i) {
2578// p = t - 1;
2579// t = A[i];
2580// if(0 <= p) {
2581// c = T[j];
2582// U[j] = (j <= p) ? T[p] : (uint8_t)A[p];
2583// A[j] = c;
2584// j++;
2585// } else {
2586// *idx = i;
2587// }
2588// }
2589// p = t - 1;
2590// if(0 <= p) {
2591// c = T[j];
2592// U[j] = (j <= p) ? T[p] : (uint8_t)A[p];
2593// A[j] = c;
2594// } else {
2595// *idx = i;
2596// }
2597// } else {
2598// U[0] = T[n - 1];
2599// for(i = 0; A[i] != 0; ++i) { U[i + 1] = T[A[i] - 1]; }
2600// *idx = i + 1;
2601// for(++i; i < n; ++i) { U[i] = T[A[i] - 1]; }
2602// }
2603//
2604// if(SA == NULL) {
2605// /* Deallocate memory. */
2606// free(A);
2607// }
2608//
2609// return 0;
2610// }
2611
2612// TODO: use?
2613/* Inverse Burrows-Wheeler transform. */
2614// template <typename saidx_t>
2615// int32_t
2616// inverse_bw_transform(const uint8_t *T, uint8_t *U, saidx_t *A,
2617// saidx_t n, saidx_t idx) {
2618// saidx_t C[ALPHABET_SIZE];
2619// uint8_t D[ALPHABET_SIZE];
2620// saidx_t *B;
2621// saidx_t i, p;
2622// int32_t c, d;
2623//
2624// /* Check arguments. */
2625// if((T == NULL) || (U == NULL) || (n < 0) || (idx < 0) ||
2626// (n < idx) || ((0 < n) && (idx == 0))) {
2627// return -1;
2628// }
2629// if(n <= 1) { return 0; }
2630//
2631// if((B = A) == NULL) {
2632// /* Allocate n*sizeof(saidx_t) bytes of memory. */
2633// if((B = (saidx_t *)malloc((size_t)n * sizeof(saidx_t))) == NULL) { return -2; }
2634// }
2635//
2636// /* Inverse BW transform. */
2637// for(c = 0; c < ALPHABET_SIZE; ++c) { C[c] = 0; }
2638// for(i = 0; i < n; ++i) { ++C[T[i]]; }
2639// for(c = 0, d = 0, i = 0; c < ALPHABET_SIZE; ++c) {
2640// p = C[c];
2641// if(0 < p) {
2642// C[c] = i;
2643// D[d++] = (uint8_t)c;
2644// i += p;
2645// }
2646// }
2647// for(i = 0; i < idx; ++i) { B[C[T[i]]++] = i; }
2648// for( ; i < n; ++i) { B[C[T[i]]++] = i + 1; }
2649// for(c = 0; c < d; ++c) { C[c] = C[D[c]]; }
2650// for(i = 0, p = idx; i < n; ++i) {
2651// U[i] = D[binarysearch_lower(C, d, p)];
2652// p = B[p - 1];
2653// }
2654//
2655// if(A == NULL) {
2656// /* Deallocate memory. */
2657// free(B);
2658// }
2659//
2660// return 0;
2661// }
2662
2663} // namespace sdsl
2664
2665#endif
#define STACK_POP5(_a, _b, _c, _d, _e)
#define SS_BLOCKSIZE
#define STACK_PUSH5(_a, _b, _c, _d, _e)
#define TR_INSERTIONSORT_THRESHOLD
#define GETIDX(a)
#define BUCKET_BSTAR(_c0, _c1)
#define SS_INSERTIONSORT_THRESHOLD
#define BUCKET_A_SIZE
#define ALPHABET_SIZE
#define BUCKET_B_SIZE
#define SS_MISORT_STACKSIZE
#define STACK_POP(_a, _b, _c, _d)
#define BUCKET_A(_c0)
#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)