1/*
2 * divsufsort.c for libdivsufsort-lite
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/*- Compiler specifics -*/
28#ifdef __clang__
29#pragma clang diagnostic ignored "-Wshorten-64-to-32"
30#endif
31
32#if defined(_MSC_VER)
33# pragma warning(disable : 4244)
34# pragma warning(disable : 4127) /* C4127 : Condition expression is constant */
35#endif
36
37
38/*- Dependencies -*/
39#include <assert.h>
40#include <stdio.h>
41#include <stdlib.h>
42
43#include "divsufsort.h"
44
45/*- Constants -*/
46#if defined(INLINE)
47# undef INLINE
48#endif
49#if !defined(INLINE)
50# define INLINE __inline
51#endif
52#if defined(ALPHABET_SIZE) && (ALPHABET_SIZE < 1)
53# undef ALPHABET_SIZE
54#endif
55#if !defined(ALPHABET_SIZE)
56# define ALPHABET_SIZE (256)
57#endif
58#define BUCKET_A_SIZE (ALPHABET_SIZE)
59#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
60#if defined(SS_INSERTIONSORT_THRESHOLD)
61# if SS_INSERTIONSORT_THRESHOLD < 1
62# undef SS_INSERTIONSORT_THRESHOLD
63# define SS_INSERTIONSORT_THRESHOLD (1)
64# endif
65#else
66# define SS_INSERTIONSORT_THRESHOLD (8)
67#endif
68#if defined(SS_BLOCKSIZE)
69# if SS_BLOCKSIZE < 0
70# undef SS_BLOCKSIZE
71# define SS_BLOCKSIZE (0)
72# elif 32768 <= SS_BLOCKSIZE
73# undef SS_BLOCKSIZE
74# define SS_BLOCKSIZE (32767)
75# endif
76#else
77# define SS_BLOCKSIZE (1024)
78#endif
79/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */
80#if SS_BLOCKSIZE == 0
81# define SS_MISORT_STACKSIZE (96)
82#elif SS_BLOCKSIZE <= 4096
83# define SS_MISORT_STACKSIZE (16)
84#else
85# define SS_MISORT_STACKSIZE (24)
86#endif
87#define SS_SMERGE_STACKSIZE (32)
88#define TR_INSERTIONSORT_THRESHOLD (8)
89#define TR_STACKSIZE (64)
90
91
92/*- Macros -*/
93#ifndef SWAP
94# define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0)
95#endif /* SWAP */
96#ifndef MIN
97# define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b))
98#endif /* MIN */
99#ifndef MAX
100# define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b))
101#endif /* MAX */
102#define STACK_PUSH(_a, _b, _c, _d)\
103 do {\
104 assert(ssize < STACK_SIZE);\
105 stack[ssize].a = (_a), stack[ssize].b = (_b),\
106 stack[ssize].c = (_c), stack[ssize++].d = (_d);\
107 } while(0)
108#define STACK_PUSH5(_a, _b, _c, _d, _e)\
109 do {\
110 assert(ssize < STACK_SIZE);\
111 stack[ssize].a = (_a), stack[ssize].b = (_b),\
112 stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\
113 } while(0)
114#define STACK_POP(_a, _b, _c, _d)\
115 do {\
116 assert(0 <= ssize);\
117 if(ssize == 0) { return; }\
118 (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
119 (_c) = stack[ssize].c, (_d) = stack[ssize].d;\
120 } while(0)
121#define STACK_POP5(_a, _b, _c, _d, _e)\
122 do {\
123 assert(0 <= ssize);\
124 if(ssize == 0) { return; }\
125 (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
126 (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\
127 } while(0)
128#define BUCKET_A(_c0) bucket_A[(_c0)]
129#if ALPHABET_SIZE == 256
130#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
131#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
132#else
133#define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)])
134#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)])
135#endif
136
137
138/*- Private Functions -*/
139
140static const int lg_table[256]= {
141 -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,
142 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
143 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,
144 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,
145 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
146 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
147 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
148 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
149};
150
151#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
152
153static INLINE
154int
155ss_ilg(int n) {
156#if SS_BLOCKSIZE == 0
157 return (n & 0xffff0000) ?
158 ((n & 0xff000000) ?
159 24 + lg_table[(n >> 24) & 0xff] :
160 16 + lg_table[(n >> 16) & 0xff]) :
161 ((n & 0x0000ff00) ?
162 8 + lg_table[(n >> 8) & 0xff] :
163 0 + lg_table[(n >> 0) & 0xff]);
164#elif SS_BLOCKSIZE < 256
165 return lg_table[n];
166#else
167 return (n & 0xff00) ?
168 8 + lg_table[(n >> 8) & 0xff] :
169 0 + lg_table[(n >> 0) & 0xff];
170#endif
171}
172
173#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
174
175#if SS_BLOCKSIZE != 0
176
177static const int sqq_table[256] = {
178 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61,
179 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89,
180 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109,
181110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
182128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
183143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155,
184156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
185169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180,
186181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191,
187192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201,
188202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211,
189212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221,
190221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230,
191230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
192239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247,
193247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255
194};
195
196static INLINE
197int
198ss_isqrt(int x) {
199 int y, e;
200
201 if(x >= (SS_BLOCKSIZE * SS_BLOCKSIZE)) { return SS_BLOCKSIZE; }
202 e = (x & 0xffff0000) ?
203 ((x & 0xff000000) ?
204 24 + lg_table[(x >> 24) & 0xff] :
205 16 + lg_table[(x >> 16) & 0xff]) :
206 ((x & 0x0000ff00) ?
207 8 + lg_table[(x >> 8) & 0xff] :
208 0 + lg_table[(x >> 0) & 0xff]);
209
210 if(e >= 16) {
211 y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
212 if(e >= 24) { y = (y + 1 + x / y) >> 1; }
213 y = (y + 1 + x / y) >> 1;
214 } else if(e >= 8) {
215 y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
216 } else {
217 return sqq_table[x] >> 4;
218 }
219
220 return (x < (y * y)) ? y - 1 : y;
221}
222
223#endif /* SS_BLOCKSIZE != 0 */
224
225
226/*---------------------------------------------------------------------------*/
227
228/* Compares two suffixes. */
229static INLINE
230int
231ss_compare(const unsigned char *T,
232 const int *p1, const int *p2,
233 int depth) {
234 const unsigned char *U1, *U2, *U1n, *U2n;
235
236 for(U1 = T + depth + *p1,
237 U2 = T + depth + *p2,
238 U1n = T + *(p1 + 1) + 2,
239 U2n = T + *(p2 + 1) + 2;
240 (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
241 ++U1, ++U2) {
242 }
243
244 return U1 < U1n ?
245 (U2 < U2n ? *U1 - *U2 : 1) :
246 (U2 < U2n ? -1 : 0);
247}
248
249
250/*---------------------------------------------------------------------------*/
251
252#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
253
254/* Insertionsort for small size groups */
255static
256void
257ss_insertionsort(const unsigned char *T, const int *PA,
258 int *first, int *last, int depth) {
259 int *i, *j;
260 int t;
261 int r;
262
263 for(i = last - 2; first <= i; --i) {
264 for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) {
265 do { *(j - 1) = *j; } while((++j < last) && (*j < 0));
266 if(last <= j) { break; }
267 }
268 if(r == 0) { *j = ~*j; }
269 *(j - 1) = t;
270 }
271}
272
273#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */
274
275
276/*---------------------------------------------------------------------------*/
277
278#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
279
280static INLINE
281void
282ss_fixdown(const unsigned char *Td, const int *PA,
283 int *SA, int i, int size) {
284 int j, k;
285 int v;
286 int c, d, e;
287
288 for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
289 d = Td[PA[SA[k = j++]]];
290 if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; }
291 if(d <= c) { break; }
292 }
293 SA[i] = v;
294}
295
296/* Simple top-down heapsort. */
297static
298void
299ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) {
300 int i, m;
301 int t;
302
303 m = size;
304 if((size % 2) == 0) {
305 m--;
306 if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); }
307 }
308
309 for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); }
310 if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); }
311 for(i = m - 1; 0 < i; --i) {
312 t = SA[0], SA[0] = SA[i];
313 ss_fixdown(Td, PA, SA, 0, i);
314 SA[i] = t;
315 }
316}
317
318
319/*---------------------------------------------------------------------------*/
320
321/* Returns the median of three elements. */
322static INLINE
323int *
324ss_median3(const unsigned char *Td, const int *PA,
325 int *v1, int *v2, int *v3) {
326 int *t;
327 if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); }
328 if(Td[PA[*v2]] > Td[PA[*v3]]) {
329 if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; }
330 else { return v3; }
331 }
332 return v2;
333}
334
335/* Returns the median of five elements. */
336static INLINE
337int *
338ss_median5(const unsigned char *Td, const int *PA,
339 int *v1, int *v2, int *v3, int *v4, int *v5) {
340 int *t;
341 if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); }
342 if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); }
343 if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); }
344 if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); }
345 if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); }
346 if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; }
347 return v3;
348}
349
350/* Returns the pivot element. */
351static INLINE
352int *
353ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) {
354 int *middle;
355 int t;
356
357 t = last - first;
358 middle = first + t / 2;
359
360 if(t <= 512) {
361 if(t <= 32) {
362 return ss_median3(Td, PA, first, middle, last - 1);
363 } else {
364 t >>= 2;
365 return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
366 }
367 }
368 t >>= 3;
369 first = ss_median3(Td, PA, first, first + t, first + (t << 1));
370 middle = ss_median3(Td, PA, middle - t, middle, middle + t);
371 last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
372 return ss_median3(Td, PA, first, middle, last);
373}
374
375
376/*---------------------------------------------------------------------------*/
377
378/* Binary partition for substrings. */
379static INLINE
380int *
381ss_partition(const int *PA,
382 int *first, int *last, int depth) {
383 int *a, *b;
384 int t;
385 for(a = first - 1, b = last;;) {
386 for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
387 for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { }
388 if(b <= a) { break; }
389 t = ~*b;
390 *b = *a;
391 *a = t;
392 }
393 if(first < a) { *first = ~*first; }
394 return a;
395}
396
397/* Multikey introsort for medium size groups. */
398static
399void
400ss_mintrosort(const unsigned char *T, const int *PA,
401 int *first, int *last,
402 int depth) {
403#define STACK_SIZE SS_MISORT_STACKSIZE
404 struct { int *a, *b, c; int d; } stack[STACK_SIZE];
405 const unsigned char *Td;
406 int *a, *b, *c, *d, *e, *f;
407 int s, t;
408 int ssize;
409 int limit;
410 int v, x = 0;
411
412 for(ssize = 0, limit = ss_ilg(last - first);;) {
413
414 if((last - first) <= SS_INSERTIONSORT_THRESHOLD) {
415#if 1 < SS_INSERTIONSORT_THRESHOLD
416 if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); }
417#endif
418 STACK_POP(first, last, depth, limit);
419 continue;
420 }
421
422 Td = T + depth;
423 if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); }
424 if(limit < 0) {
425 for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) {
426 if((x = Td[PA[*a]]) != v) {
427 if(1 < (a - first)) { break; }
428 v = x;
429 first = a;
430 }
431 }
432 if(Td[PA[*first] - 1] < v) {
433 first = ss_partition(PA, first, a, depth);
434 }
435 if((a - first) <= (last - a)) {
436 if(1 < (a - first)) {
437 STACK_PUSH(a, last, depth, -1);
438 last = a, depth += 1, limit = ss_ilg(a - first);
439 } else {
440 first = a, limit = -1;
441 }
442 } else {
443 if(1 < (last - a)) {
444 STACK_PUSH(first, a, depth + 1, ss_ilg(a - first));
445 first = a, limit = -1;
446 } else {
447 last = a, depth += 1, limit = ss_ilg(a - first);
448 }
449 }
450 continue;
451 }
452
453 /* choose pivot */
454 a = ss_pivot(Td, PA, first, last);
455 v = Td[PA[*a]];
456 SWAP(*first, *a);
457
458 /* partition */
459 for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { }
460 if(((a = b) < last) && (x < v)) {
461 for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) {
462 if(x == v) { SWAP(*b, *a); ++a; }
463 }
464 }
465 for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { }
466 if((b < (d = c)) && (x > v)) {
467 for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
468 if(x == v) { SWAP(*c, *d); --d; }
469 }
470 }
471 for(; b < c;) {
472 SWAP(*b, *c);
473 for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) {
474 if(x == v) { SWAP(*b, *a); ++a; }
475 }
476 for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
477 if(x == v) { SWAP(*c, *d); --d; }
478 }
479 }
480
481 if(a <= d) {
482 c = b - 1;
483
484 if((s = a - first) > (t = b - a)) { s = t; }
485 for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
486 if((s = d - c) > (t = last - d - 1)) { s = t; }
487 for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
488
489 a = first + (b - a), c = last - (d - c);
490 b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth);
491
492 if((a - first) <= (last - c)) {
493 if((last - c) <= (c - b)) {
494 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
495 STACK_PUSH(c, last, depth, limit);
496 last = a;
497 } else if((a - first) <= (c - b)) {
498 STACK_PUSH(c, last, depth, limit);
499 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
500 last = a;
501 } else {
502 STACK_PUSH(c, last, depth, limit);
503 STACK_PUSH(first, a, depth, limit);
504 first = b, last = c, depth += 1, limit = ss_ilg(c - b);
505 }
506 } else {
507 if((a - first) <= (c - b)) {
508 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
509 STACK_PUSH(first, a, depth, limit);
510 first = c;
511 } else if((last - c) <= (c - b)) {
512 STACK_PUSH(first, a, depth, limit);
513 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
514 first = c;
515 } else {
516 STACK_PUSH(first, a, depth, limit);
517 STACK_PUSH(c, last, depth, limit);
518 first = b, last = c, depth += 1, limit = ss_ilg(c - b);
519 }
520 }
521 } else {
522 limit += 1;
523 if(Td[PA[*first] - 1] < v) {
524 first = ss_partition(PA, first, last, depth);
525 limit = ss_ilg(last - first);
526 }
527 depth += 1;
528 }
529 }
530#undef STACK_SIZE
531}
532
533#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
534
535
536/*---------------------------------------------------------------------------*/
537
538#if SS_BLOCKSIZE != 0
539
540static INLINE
541void
542ss_blockswap(int *a, int *b, int n) {
543 int t;
544 for(; 0 < n; --n, ++a, ++b) {
545 t = *a, *a = *b, *b = t;
546 }
547}
548
549static INLINE
550void
551ss_rotate(int *first, int *middle, int *last) {
552 int *a, *b, t;
553 int l, r;
554 l = middle - first, r = last - middle;
555 for(; (0 < l) && (0 < r);) {
556 if(l == r) { ss_blockswap(first, middle, l); break; }
557 if(l < r) {
558 a = last - 1, b = middle - 1;
559 t = *a;
560 do {
561 *a-- = *b, *b-- = *a;
562 if(b < first) {
563 *a = t;
564 last = a;
565 if((r -= l + 1) <= l) { break; }
566 a -= 1, b = middle - 1;
567 t = *a;
568 }
569 } while(1);
570 } else {
571 a = first, b = middle;
572 t = *a;
573 do {
574 *a++ = *b, *b++ = *a;
575 if(last <= b) {
576 *a = t;
577 first = a + 1;
578 if((l -= r + 1) <= r) { break; }
579 a += 1, b = middle;
580 t = *a;
581 }
582 } while(1);
583 }
584 }
585}
586
587
588/*---------------------------------------------------------------------------*/
589
590static
591void
592ss_inplacemerge(const unsigned char *T, const int *PA,
593 int *first, int *middle, int *last,
594 int depth) {
595 const int *p;
596 int *a, *b;
597 int len, half;
598 int q, r;
599 int x;
600
601 for(;;) {
602 if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); }
603 else { x = 0; p = PA + *(last - 1); }
604 for(a = first, len = middle - first, half = len >> 1, r = -1;
605 0 < len;
606 len = half, half >>= 1) {
607 b = a + half;
608 q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
609 if(q < 0) {
610 a = b + 1;
611 half -= (len & 1) ^ 1;
612 } else {
613 r = q;
614 }
615 }
616 if(a < middle) {
617 if(r == 0) { *a = ~*a; }
618 ss_rotate(a, middle, last);
619 last -= middle - a;
620 middle = a;
621 if(first == middle) { break; }
622 }
623 --last;
624 if(x != 0) { while(*--last < 0) { } }
625 if(middle == last) { break; }
626 }
627}
628
629
630/*---------------------------------------------------------------------------*/
631
632/* Merge-forward with internal buffer. */
633static
634void
635ss_mergeforward(const unsigned char *T, const int *PA,
636 int *first, int *middle, int *last,
637 int *buf, int depth) {
638 int *a, *b, *c, *bufend;
639 int t;
640 int r;
641
642 bufend = buf + (middle - first) - 1;
643 ss_blockswap(buf, first, middle - first);
644
645 for(t = *(a = first), b = buf, c = middle;;) {
646 r = ss_compare(T, PA + *b, PA + *c, depth);
647 if(r < 0) {
648 do {
649 *a++ = *b;
650 if(bufend <= b) { *bufend = t; return; }
651 *b++ = *a;
652 } while(*b < 0);
653 } else if(r > 0) {
654 do {
655 *a++ = *c, *c++ = *a;
656 if(last <= c) {
657 while(b < bufend) { *a++ = *b, *b++ = *a; }
658 *a = *b, *b = t;
659 return;
660 }
661 } while(*c < 0);
662 } else {
663 *c = ~*c;
664 do {
665 *a++ = *b;
666 if(bufend <= b) { *bufend = t; return; }
667 *b++ = *a;
668 } while(*b < 0);
669
670 do {
671 *a++ = *c, *c++ = *a;
672 if(last <= c) {
673 while(b < bufend) { *a++ = *b, *b++ = *a; }
674 *a = *b, *b = t;
675 return;
676 }
677 } while(*c < 0);
678 }
679 }
680}
681
682/* Merge-backward with internal buffer. */
683static
684void
685ss_mergebackward(const unsigned char *T, const int *PA,
686 int *first, int *middle, int *last,
687 int *buf, int depth) {
688 const int *p1, *p2;
689 int *a, *b, *c, *bufend;
690 int t;
691 int r;
692 int x;
693
694 bufend = buf + (last - middle) - 1;
695 ss_blockswap(buf, middle, last - middle);
696
697 x = 0;
698 if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; }
699 else { p1 = PA + *bufend; }
700 if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; }
701 else { p2 = PA + *(middle - 1); }
702 for(t = *(a = last - 1), b = bufend, c = middle - 1;;) {
703 r = ss_compare(T, p1, p2, depth);
704 if(0 < r) {
705 if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
706 *a-- = *b;
707 if(b <= buf) { *buf = t; break; }
708 *b-- = *a;
709 if(*b < 0) { p1 = PA + ~*b; x |= 1; }
710 else { p1 = PA + *b; }
711 } else if(r < 0) {
712 if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
713 *a-- = *c, *c-- = *a;
714 if(c < first) {
715 while(buf < b) { *a-- = *b, *b-- = *a; }
716 *a = *b, *b = t;
717 break;
718 }
719 if(*c < 0) { p2 = PA + ~*c; x |= 2; }
720 else { p2 = PA + *c; }
721 } else {
722 if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
723 *a-- = ~*b;
724 if(b <= buf) { *buf = t; break; }
725 *b-- = *a;
726 if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
727 *a-- = *c, *c-- = *a;
728 if(c < first) {
729 while(buf < b) { *a-- = *b, *b-- = *a; }
730 *a = *b, *b = t;
731 break;
732 }
733 if(*b < 0) { p1 = PA + ~*b; x |= 1; }
734 else { p1 = PA + *b; }
735 if(*c < 0) { p2 = PA + ~*c; x |= 2; }
736 else { p2 = PA + *c; }
737 }
738 }
739}
740
741/* D&C based merge. */
742static
743void
744ss_swapmerge(const unsigned char *T, const int *PA,
745 int *first, int *middle, int *last,
746 int *buf, int bufsize, int depth) {
747#define STACK_SIZE SS_SMERGE_STACKSIZE
748#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
749#define MERGE_CHECK(a, b, c)\
750 do {\
751 if(((c) & 1) ||\
752 (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\
753 *(a) = ~*(a);\
754 }\
755 if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\
756 *(b) = ~*(b);\
757 }\
758 } while(0)
759 struct { int *a, *b, *c; int d; } stack[STACK_SIZE];
760 int *l, *r, *lm, *rm;
761 int m, len, half;
762 int ssize;
763 int check, next;
764
765 for(check = 0, ssize = 0;;) {
766 if((last - middle) <= bufsize) {
767 if((first < middle) && (middle < last)) {
768 ss_mergebackward(T, PA, first, middle, last, buf, depth);
769 }
770 MERGE_CHECK(first, last, check);
771 STACK_POP(first, middle, last, check);
772 continue;
773 }
774
775 if((middle - first) <= bufsize) {
776 if(first < middle) {
777 ss_mergeforward(T, PA, first, middle, last, buf, depth);
778 }
779 MERGE_CHECK(first, last, check);
780 STACK_POP(first, middle, last, check);
781 continue;
782 }
783
784 for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1;
785 0 < len;
786 len = half, half >>= 1) {
787 if(ss_compare(T, PA + GETIDX(*(middle + m + half)),
788 PA + GETIDX(*(middle - m - half - 1)), depth) < 0) {
789 m += half + 1;
790 half -= (len & 1) ^ 1;
791 }
792 }
793
794 if(0 < m) {
795 lm = middle - m, rm = middle + m;
796 ss_blockswap(lm, middle, m);
797 l = r = middle, next = 0;
798 if(rm < last) {
799 if(*rm < 0) {
800 *rm = ~*rm;
801 if(first < lm) { for(; *--l < 0;) { } next |= 4; }
802 next |= 1;
803 } else if(first < lm) {
804 for(; *r < 0; ++r) { }
805 next |= 2;
806 }
807 }
808
809 if((l - first) <= (last - r)) {
810 STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
811 middle = lm, last = l, check = (check & 3) | (next & 4);
812 } else {
813 if((next & 2) && (r == middle)) { next ^= 6; }
814 STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
815 first = r, middle = rm, check = (next & 3) | (check & 4);
816 }
817 } else {
818 if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) {
819 *middle = ~*middle;
820 }
821 MERGE_CHECK(first, last, check);
822 STACK_POP(first, middle, last, check);
823 }
824 }
825#undef STACK_SIZE
826}
827
828#endif /* SS_BLOCKSIZE != 0 */
829
830
831/*---------------------------------------------------------------------------*/
832
833/* Substring sort */
834static
835void
836sssort(const unsigned char *T, const int *PA,
837 int *first, int *last,
838 int *buf, int bufsize,
839 int depth, int n, int lastsuffix) {
840 int *a;
841#if SS_BLOCKSIZE != 0
842 int *b, *middle, *curbuf;
843 int j, k, curbufsize, limit;
844#endif
845 int i;
846
847 if(lastsuffix != 0) { ++first; }
848
849#if SS_BLOCKSIZE == 0
850 ss_mintrosort(T, PA, first, last, depth);
851#else
852 if((bufsize < SS_BLOCKSIZE) &&
853 (bufsize < (last - first)) &&
854 (bufsize < (limit = ss_isqrt(last - first)))) {
855 if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; }
856 buf = middle = last - limit, bufsize = limit;
857 } else {
858 middle = last, limit = 0;
859 }
860 for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) {
861#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
862 ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth);
863#elif 1 < SS_BLOCKSIZE
864 ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth);
865#endif
866 curbufsize = last - (a + SS_BLOCKSIZE);
867 curbuf = a + SS_BLOCKSIZE;
868 if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
869 for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) {
870 ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
871 }
872 }
873#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
874 ss_mintrosort(T, PA, a, middle, depth);
875#elif 1 < SS_BLOCKSIZE
876 ss_insertionsort(T, PA, a, middle, depth);
877#endif
878 for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) {
879 if(i & 1) {
880 ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
881 a -= k;
882 }
883 }
884 if(limit != 0) {
885#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
886 ss_mintrosort(T, PA, middle, last, depth);
887#elif 1 < SS_BLOCKSIZE
888 ss_insertionsort(T, PA, middle, last, depth);
889#endif
890 ss_inplacemerge(T, PA, first, middle, last, depth);
891 }
892#endif
893
894 if(lastsuffix != 0) {
895 /* Insert last type B* suffix. */
896 int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
897 for(a = first, i = *(first - 1);
898 (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth)));
899 ++a) {
900 *(a - 1) = *a;
901 }
902 *(a - 1) = i;
903 }
904}
905
906
907/*---------------------------------------------------------------------------*/
908
909static INLINE
910int
911tr_ilg(int n) {
912 return (n & 0xffff0000) ?
913 ((n & 0xff000000) ?
914 24 + lg_table[(n >> 24) & 0xff] :
915 16 + lg_table[(n >> 16) & 0xff]) :
916 ((n & 0x0000ff00) ?
917 8 + lg_table[(n >> 8) & 0xff] :
918 0 + lg_table[(n >> 0) & 0xff]);
919}
920
921
922/*---------------------------------------------------------------------------*/
923
924/* Simple insertionsort for small size groups. */
925static
926void
927tr_insertionsort(const int *ISAd, int *first, int *last) {
928 int *a, *b;
929 int t, r;
930
931 for(a = first + 1; a < last; ++a) {
932 for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) {
933 do { *(b + 1) = *b; } while((first <= --b) && (*b < 0));
934 if(b < first) { break; }
935 }
936 if(r == 0) { *b = ~*b; }
937 *(b + 1) = t;
938 }
939}
940
941
942/*---------------------------------------------------------------------------*/
943
944static INLINE
945void
946tr_fixdown(const int *ISAd, int *SA, int i, int size) {
947 int j, k;
948 int v;
949 int c, d, e;
950
951 for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
952 d = ISAd[SA[k = j++]];
953 if(d < (e = ISAd[SA[j]])) { k = j; d = e; }
954 if(d <= c) { break; }
955 }
956 SA[i] = v;
957}
958
959/* Simple top-down heapsort. */
960static
961void
962tr_heapsort(const int *ISAd, int *SA, int size) {
963 int i, m;
964 int t;
965
966 m = size;
967 if((size % 2) == 0) {
968 m--;
969 if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); }
970 }
971
972 for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); }
973 if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); }
974 for(i = m - 1; 0 < i; --i) {
975 t = SA[0], SA[0] = SA[i];
976 tr_fixdown(ISAd, SA, 0, i);
977 SA[i] = t;
978 }
979}
980
981
982/*---------------------------------------------------------------------------*/
983
984/* Returns the median of three elements. */
985static INLINE
986int *
987tr_median3(const int *ISAd, int *v1, int *v2, int *v3) {
988 int *t;
989 if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); }
990 if(ISAd[*v2] > ISAd[*v3]) {
991 if(ISAd[*v1] > ISAd[*v3]) { return v1; }
992 else { return v3; }
993 }
994 return v2;
995}
996
997/* Returns the median of five elements. */
998static INLINE
999int *
1000tr_median5(const int *ISAd,
1001 int *v1, int *v2, int *v3, int *v4, int *v5) {
1002 int *t;
1003 if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); }
1004 if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); }
1005 if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4); SWAP(v3, v5); }
1006 if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3); }
1007 if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4); SWAP(v3, v5); }
1008 if(ISAd[*v3] > ISAd[*v4]) { return v4; }
1009 return v3;
1010}
1011
1012/* Returns the pivot element. */
1013static INLINE
1014int *
1015tr_pivot(const int *ISAd, int *first, int *last) {
1016 int *middle;
1017 int t;
1018
1019 t = last - first;
1020 middle = first + t / 2;
1021
1022 if(t <= 512) {
1023 if(t <= 32) {
1024 return tr_median3(ISAd, first, middle, last - 1);
1025 } else {
1026 t >>= 2;
1027 return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1028 }
1029 }
1030 t >>= 3;
1031 first = tr_median3(ISAd, first, first + t, first + (t << 1));
1032 middle = tr_median3(ISAd, middle - t, middle, middle + t);
1033 last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
1034 return tr_median3(ISAd, first, middle, last);
1035}
1036
1037
1038/*---------------------------------------------------------------------------*/
1039
1040typedef struct _trbudget_t trbudget_t;
1041struct _trbudget_t {
1042 int chance;
1043 int remain;
1044 int incval;
1045 int count;
1046};
1047
1048static INLINE
1049void
1050trbudget_init(trbudget_t *budget, int chance, int incval) {
1051 budget->chance = chance;
1052 budget->remain = budget->incval = incval;
1053}
1054
1055static INLINE
1056int
1057trbudget_check(trbudget_t *budget, int size) {
1058 if(size <= budget->remain) { budget->remain -= size; return 1; }
1059 if(budget->chance == 0) { budget->count += size; return 0; }
1060 budget->remain += budget->incval - size;
1061 budget->chance -= 1;
1062 return 1;
1063}
1064
1065
1066/*---------------------------------------------------------------------------*/
1067
1068static INLINE
1069void
1070tr_partition(const int *ISAd,
1071 int *first, int *middle, int *last,
1072 int **pa, int **pb, int v) {
1073 int *a, *b, *c, *d, *e, *f;
1074 int t, s;
1075 int x = 0;
1076
1077 for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { }
1078 if(((a = b) < last) && (x < v)) {
1079 for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
1080 if(x == v) { SWAP(*b, *a); ++a; }
1081 }
1082 }
1083 for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { }
1084 if((b < (d = c)) && (x > v)) {
1085 for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1086 if(x == v) { SWAP(*c, *d); --d; }
1087 }
1088 }
1089 for(; b < c;) {
1090 SWAP(*b, *c);
1091 for(; (++b < c) && ((x = ISAd[*b]) <= v);) {
1092 if(x == v) { SWAP(*b, *a); ++a; }
1093 }
1094 for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1095 if(x == v) { SWAP(*c, *d); --d; }
1096 }
1097 }
1098
1099 if(a <= d) {
1100 c = b - 1;
1101 if((s = a - first) > (t = b - a)) { s = t; }
1102 for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1103 if((s = d - c) > (t = last - d - 1)) { s = t; }
1104 for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1105 first += (b - a), last -= (d - c);
1106 }
1107 *pa = first, *pb = last;
1108}
1109
1110static
1111void
1112tr_copy(int *ISA, const int *SA,
1113 int *first, int *a, int *b, int *last,
1114 int depth) {
1115 /* sort suffixes of middle partition
1116 by using sorted order of suffixes of left and right partition. */
1117 int *c, *d, *e;
1118 int s, v;
1119
1120 v = b - SA - 1;
1121 for(c = first, d = a - 1; c <= d; ++c) {
1122 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1123 *++d = s;
1124 ISA[s] = d - SA;
1125 }
1126 }
1127 for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1128 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1129 *--d = s;
1130 ISA[s] = d - SA;
1131 }
1132 }
1133}
1134
1135static
1136void
1137tr_partialcopy(int *ISA, const int *SA,
1138 int *first, int *a, int *b, int *last,
1139 int depth) {
1140 int *c, *d, *e;
1141 int s, v;
1142 int rank, lastrank, newrank = -1;
1143
1144 v = b - SA - 1;
1145 lastrank = -1;
1146 for(c = first, d = a - 1; c <= d; ++c) {
1147 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1148 *++d = s;
1149 rank = ISA[s + depth];
1150 if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1151 ISA[s] = newrank;
1152 }
1153 }
1154
1155 lastrank = -1;
1156 for(e = d; first <= e; --e) {
1157 rank = ISA[*e];
1158 if(lastrank != rank) { lastrank = rank; newrank = e - SA; }
1159 if(newrank != rank) { ISA[*e] = newrank; }
1160 }
1161
1162 lastrank = -1;
1163 for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1164 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1165 *--d = s;
1166 rank = ISA[s + depth];
1167 if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1168 ISA[s] = newrank;
1169 }
1170 }
1171}
1172
1173static
1174void
1175tr_introsort(int *ISA, const int *ISAd,
1176 int *SA, int *first, int *last,
1177 trbudget_t *budget) {
1178#define STACK_SIZE TR_STACKSIZE
1179 struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE];
1180 int *a, *b, *c;
1181 int t;
1182 int v, x = 0;
1183 int incr = ISAd - ISA;
1184 int limit, next;
1185 int ssize, trlink = -1;
1186
1187 for(ssize = 0, limit = tr_ilg(last - first);;) {
1188
1189 if(limit < 0) {
1190 if(limit == -1) {
1191 /* tandem repeat partition */
1192 tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1);
1193
1194 /* update ranks */
1195 if(a < last) {
1196 for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1197 }
1198 if(b < last) {
1199 for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
1200 }
1201
1202 /* push */
1203 if(1 < (b - a)) {
1204 STACK_PUSH5(NULL, a, b, 0, 0);
1205 STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1206 trlink = ssize - 2;
1207 }
1208 if((a - first) <= (last - b)) {
1209 if(1 < (a - first)) {
1210 STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink);
1211 last = a, limit = tr_ilg(a - first);
1212 } else if(1 < (last - b)) {
1213 first = b, limit = tr_ilg(last - b);
1214 } else {
1215 STACK_POP5(ISAd, first, last, limit, trlink);
1216 }
1217 } else {
1218 if(1 < (last - b)) {
1219 STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink);
1220 first = b, limit = tr_ilg(last - b);
1221 } else if(1 < (a - first)) {
1222 last = a, limit = tr_ilg(a - first);
1223 } else {
1224 STACK_POP5(ISAd, first, last, limit, trlink);
1225 }
1226 }
1227 } else if(limit == -2) {
1228 /* tandem repeat copy */
1229 a = stack[--ssize].b, b = stack[ssize].c;
1230 if(stack[ssize].d == 0) {
1231 tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
1232 } else {
1233 if(0 <= trlink) { stack[trlink].d = -1; }
1234 tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
1235 }
1236 STACK_POP5(ISAd, first, last, limit, trlink);
1237 } else {
1238 /* sorted partition */
1239 if(0 <= *first) {
1240 a = first;
1241 do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a));
1242 first = a;
1243 }
1244 if(first < last) {
1245 a = first; do { *a = ~*a; } while(*++a < 0);
1246 next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1;
1247 if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } }
1248
1249 /* push */
1250 if(trbudget_check(budget, a - first)) {
1251 if((a - first) <= (last - a)) {
1252 STACK_PUSH5(ISAd, a, last, -3, trlink);
1253 ISAd += incr, last = a, limit = next;
1254 } else {
1255 if(1 < (last - a)) {
1256 STACK_PUSH5(ISAd + incr, first, a, next, trlink);
1257 first = a, limit = -3;
1258 } else {
1259 ISAd += incr, last = a, limit = next;
1260 }
1261 }
1262 } else {
1263 if(0 <= trlink) { stack[trlink].d = -1; }
1264 if(1 < (last - a)) {
1265 first = a, limit = -3;
1266 } else {
1267 STACK_POP5(ISAd, first, last, limit, trlink);
1268 }
1269 }
1270 } else {
1271 STACK_POP5(ISAd, first, last, limit, trlink);
1272 }
1273 }
1274 continue;
1275 }
1276
1277 if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
1278 tr_insertionsort(ISAd, first, last);
1279 limit = -3;
1280 continue;
1281 }
1282
1283 if(limit-- == 0) {
1284 tr_heapsort(ISAd, first, last - first);
1285 for(a = last - 1; first < a; a = b) {
1286 for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
1287 }
1288 limit = -3;
1289 continue;
1290 }
1291
1292 /* choose pivot */
1293 a = tr_pivot(ISAd, first, last);
1294 SWAP(*first, *a);
1295 v = ISAd[*first];
1296
1297 /* partition */
1298 tr_partition(ISAd, first, first + 1, last, &a, &b, v);
1299 if((last - first) != (b - a)) {
1300 next = (ISA[*a] != v) ? tr_ilg(b - a) : -1;
1301
1302 /* update ranks */
1303 for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1304 if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } }
1305
1306 /* push */
1307 if((1 < (b - a)) && (trbudget_check(budget, b - a))) {
1308 if((a - first) <= (last - b)) {
1309 if((last - b) <= (b - a)) {
1310 if(1 < (a - first)) {
1311 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1312 STACK_PUSH5(ISAd, b, last, limit, trlink);
1313 last = a;
1314 } else if(1 < (last - b)) {
1315 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1316 first = b;
1317 } else {
1318 ISAd += incr, first = a, last = b, limit = next;
1319 }
1320 } else if((a - first) <= (b - a)) {
1321 if(1 < (a - first)) {
1322 STACK_PUSH5(ISAd, b, last, limit, trlink);
1323 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1324 last = a;
1325 } else {
1326 STACK_PUSH5(ISAd, b, last, limit, trlink);
1327 ISAd += incr, first = a, last = b, limit = next;
1328 }
1329 } else {
1330 STACK_PUSH5(ISAd, b, last, limit, trlink);
1331 STACK_PUSH5(ISAd, first, a, limit, trlink);
1332 ISAd += incr, first = a, last = b, limit = next;
1333 }
1334 } else {
1335 if((a - first) <= (b - a)) {
1336 if(1 < (last - b)) {
1337 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1338 STACK_PUSH5(ISAd, first, a, limit, trlink);
1339 first = b;
1340 } else if(1 < (a - first)) {
1341 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1342 last = a;
1343 } else {
1344 ISAd += incr, first = a, last = b, limit = next;
1345 }
1346 } else if((last - b) <= (b - a)) {
1347 if(1 < (last - b)) {
1348 STACK_PUSH5(ISAd, first, a, limit, trlink);
1349 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1350 first = b;
1351 } else {
1352 STACK_PUSH5(ISAd, first, a, limit, trlink);
1353 ISAd += incr, first = a, last = b, limit = next;
1354 }
1355 } else {
1356 STACK_PUSH5(ISAd, first, a, limit, trlink);
1357 STACK_PUSH5(ISAd, b, last, limit, trlink);
1358 ISAd += incr, first = a, last = b, limit = next;
1359 }
1360 }
1361 } else {
1362 if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
1363 if((a - first) <= (last - b)) {
1364 if(1 < (a - first)) {
1365 STACK_PUSH5(ISAd, b, last, limit, trlink);
1366 last = a;
1367 } else if(1 < (last - b)) {
1368 first = b;
1369 } else {
1370 STACK_POP5(ISAd, first, last, limit, trlink);
1371 }
1372 } else {
1373 if(1 < (last - b)) {
1374 STACK_PUSH5(ISAd, first, a, limit, trlink);
1375 first = b;
1376 } else if(1 < (a - first)) {
1377 last = a;
1378 } else {
1379 STACK_POP5(ISAd, first, last, limit, trlink);
1380 }
1381 }
1382 }
1383 } else {
1384 if(trbudget_check(budget, last - first)) {
1385 limit = tr_ilg(last - first), ISAd += incr;
1386 } else {
1387 if(0 <= trlink) { stack[trlink].d = -1; }
1388 STACK_POP5(ISAd, first, last, limit, trlink);
1389 }
1390 }
1391 }
1392#undef STACK_SIZE
1393}
1394
1395
1396
1397/*---------------------------------------------------------------------------*/
1398
1399/* Tandem repeat sort */
1400static
1401void
1402trsort(int *ISA, int *SA, int n, int depth) {
1403 int *ISAd;
1404 int *first, *last;
1405 trbudget_t budget;
1406 int t, skip, unsorted;
1407
1408 trbudget_init(&budget, tr_ilg(n) * 2 / 3, n);
1409/* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */
1410 for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) {
1411 first = SA;
1412 skip = 0;
1413 unsorted = 0;
1414 do {
1415 if((t = *first) < 0) { first -= t; skip += t; }
1416 else {
1417 if(skip != 0) { *(first + skip) = skip; skip = 0; }
1418 last = SA + ISA[t] + 1;
1419 if(1 < (last - first)) {
1420 budget.count = 0;
1421 tr_introsort(ISA, ISAd, SA, first, last, &budget);
1422 if(budget.count != 0) { unsorted += budget.count; }
1423 else { skip = first - last; }
1424 } else if((last - first) == 1) {
1425 skip = -1;
1426 }
1427 first = last;
1428 }
1429 } while(first < (SA + n));
1430 if(skip != 0) { *(first + skip) = skip; }
1431 if(unsorted == 0) { break; }
1432 }
1433}
1434
1435
1436/*---------------------------------------------------------------------------*/
1437
1438/* Sorts suffixes of type B*. */
1439static
1440int
1441sort_typeBstar(const unsigned char *T, int *SA,
1442 int *bucket_A, int *bucket_B,
1443 int n, int openMP) {
1444 int *PAb, *ISAb, *buf;
1445#ifdef LIBBSC_OPENMP
1446 int *curbuf;
1447 int l;
1448#endif
1449 int i, j, k, t, m, bufsize;
1450 int c0, c1;
1451#ifdef LIBBSC_OPENMP
1452 int d0, d1;
1453#endif
1454 (void)openMP;
1455
1456 /* Initialize bucket arrays. */
1457 for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; }
1458 for(i = 0; i < BUCKET_B_SIZE; ++i) { bucket_B[i] = 0; }
1459
1460 /* Count the number of occurrences of the first one or two characters of each
1461 type A, B and B* suffix. Moreover, store the beginning position of all
1462 type B* suffixes into the array SA. */
1463 for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) {
1464 /* type A suffix. */
1465 do { ++BUCKET_A(c1 = c0); } while((0 <= --i) && ((c0 = T[i]) >= c1));
1466 if(0 <= i) {
1467 /* type B* suffix. */
1468 ++BUCKET_BSTAR(c0, c1);
1469 SA[--m] = i;
1470 /* type B suffix. */
1471 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) {
1472 ++BUCKET_B(c0, c1);
1473 }
1474 }
1475 }
1476 m = n - m;
1477/*
1478note:
1479 A type B* suffix is lexicographically smaller than a type B suffix that
1480 begins with the same first two characters.
1481*/
1482
1483 /* Calculate the index of start/end point of each bucket. */
1484 for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0) {
1485 t = i + BUCKET_A(c0);
1486 BUCKET_A(c0) = i + j; /* start point */
1487 i = t + BUCKET_B(c0, c0);
1488 for(c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1) {
1489 j += BUCKET_BSTAR(c0, c1);
1490 BUCKET_BSTAR(c0, c1) = j; /* end point */
1491 i += BUCKET_B(c0, c1);
1492 }
1493 }
1494
1495 if(0 < m) {
1496 /* Sort the type B* suffixes by their first two characters. */
1497 PAb = SA + n - m; ISAb = SA + m;
1498 for(i = m - 2; 0 <= i; --i) {
1499 t = PAb[i], c0 = T[t], c1 = T[t + 1];
1500 SA[--BUCKET_BSTAR(c0, c1)] = i;
1501 }
1502 t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
1503 SA[--BUCKET_BSTAR(c0, c1)] = m - 1;
1504
1505 /* Sort the type B* substrings using sssort. */
1506#ifdef LIBBSC_OPENMP
1507 if (openMP)
1508 {
1509 buf = SA + m;
1510 c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m;
1511#pragma omp parallel default(shared) private(bufsize, curbuf, k, l, d0, d1)
1512 {
1513 bufsize = (n - (2 * m)) / omp_get_num_threads();
1514 curbuf = buf + omp_get_thread_num() * bufsize;
1515 k = 0;
1516 for(;;) {
1517 #pragma omp critical(sssort_lock)
1518 {
1519 if(0 < (l = j)) {
1520 d0 = c0, d1 = c1;
1521 do {
1522 k = BUCKET_BSTAR(d0, d1);
1523 if(--d1 <= d0) {
1524 d1 = ALPHABET_SIZE - 1;
1525 if(--d0 < 0) { break; }
1526 }
1527 } while(((l - k) <= 1) && (0 < (l = k)));
1528 c0 = d0, c1 = d1, j = k;
1529 }
1530 }
1531 if(l == 0) { break; }
1532 sssort(T, PAb, SA + k, SA + l,
1533 curbuf, bufsize, 2, n, *(SA + k) == (m - 1));
1534 }
1535 }
1536 }
1537 else
1538 {
1539 buf = SA + m, bufsize = n - (2 * m);
1540 for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
1541 for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
1542 i = BUCKET_BSTAR(c0, c1);
1543 if(1 < (j - i)) {
1544 sssort(T, PAb, SA + i, SA + j,
1545 buf, bufsize, 2, n, *(SA + i) == (m - 1));
1546 }
1547 }
1548 }
1549 }
1550#else
1551 buf = SA + m, bufsize = n - (2 * m);
1552 for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
1553 for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
1554 i = BUCKET_BSTAR(c0, c1);
1555 if(1 < (j - i)) {
1556 sssort(T, PAb, SA + i, SA + j,
1557 buf, bufsize, 2, n, *(SA + i) == (m - 1));
1558 }
1559 }
1560 }
1561#endif
1562
1563 /* Compute ranks of type B* substrings. */
1564 for(i = m - 1; 0 <= i; --i) {
1565 if(0 <= SA[i]) {
1566 j = i;
1567 do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i]));
1568 SA[i + 1] = i - j;
1569 if(i <= 0) { break; }
1570 }
1571 j = i;
1572 do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0);
1573 ISAb[SA[i]] = j;
1574 }
1575
1576 /* Construct the inverse suffix array of type B* suffixes using trsort. */
1577 trsort(ISAb, SA, m, 1);
1578
1579 /* Set the sorted order of tyoe B* suffixes. */
1580 for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) {
1581 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { }
1582 if(0 <= i) {
1583 t = i;
1584 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { }
1585 SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
1586 }
1587 }
1588
1589 /* Calculate the index of start/end point of each bucket. */
1590 BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */
1591 for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) {
1592 i = BUCKET_A(c0 + 1) - 1;
1593 for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) {
1594 t = i - BUCKET_B(c0, c1);
1595 BUCKET_B(c0, c1) = i; /* end point */
1596
1597 /* Move all type B* suffixes to the correct position. */
1598 for(i = t, j = BUCKET_BSTAR(c0, c1);
1599 j <= k;
1600 --i, --k) { SA[i] = SA[k]; }
1601 }
1602 BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */
1603 BUCKET_B(c0, c0) = i; /* end point */
1604 }
1605 }
1606
1607 return m;
1608}
1609
1610/* Constructs the suffix array by using the sorted order of type B* suffixes. */
1611static
1612void
1613construct_SA(const unsigned char *T, int *SA,
1614 int *bucket_A, int *bucket_B,
1615 int n, int m) {
1616 int *i, *j, *k;
1617 int s;
1618 int c0, c1, c2;
1619
1620 if(0 < m) {
1621 /* Construct the sorted order of type B suffixes by using
1622 the sorted order of type B* suffixes. */
1623 for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1624 /* Scan the suffix array from right to left. */
1625 for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1626 j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1627 i <= j;
1628 --j) {
1629 if(0 < (s = *j)) {
1630 assert(T[s] == c1);
1631 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1632 assert(T[s - 1] <= T[s]);
1633 *j = ~s;
1634 c0 = T[--s];
1635 if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1636 if(c0 != c2) {
1637 if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1638 k = SA + BUCKET_B(c2 = c0, c1);
1639 }
1640 assert(k < j);
1641 *k-- = s;
1642 } else {
1643 assert(((s == 0) && (T[s] == c1)) || (s < 0));
1644 *j = ~s;
1645 }
1646 }
1647 }
1648 }
1649
1650 /* Construct the suffix array by using
1651 the sorted order of type B suffixes. */
1652 k = SA + BUCKET_A(c2 = T[n - 1]);
1653 *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
1654 /* Scan the suffix array from left to right. */
1655 for(i = SA, j = SA + n; i < j; ++i) {
1656 if(0 < (s = *i)) {
1657 assert(T[s - 1] >= T[s]);
1658 c0 = T[--s];
1659 if((s == 0) || (T[s - 1] < c0)) { s = ~s; }
1660 if(c0 != c2) {
1661 BUCKET_A(c2) = k - SA;
1662 k = SA + BUCKET_A(c2 = c0);
1663 }
1664 assert(i < k);
1665 *k++ = s;
1666 } else {
1667 assert(s < 0);
1668 *i = ~s;
1669 }
1670 }
1671}
1672
1673/* Constructs the burrows-wheeler transformed string directly
1674 by using the sorted order of type B* suffixes. */
1675static
1676int
1677construct_BWT(const unsigned char *T, int *SA,
1678 int *bucket_A, int *bucket_B,
1679 int n, int m) {
1680 int *i, *j, *k, *orig;
1681 int s;
1682 int c0, c1, c2;
1683
1684 if(0 < m) {
1685 /* Construct the sorted order of type B suffixes by using
1686 the sorted order of type B* suffixes. */
1687 for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1688 /* Scan the suffix array from right to left. */
1689 for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1690 j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1691 i <= j;
1692 --j) {
1693 if(0 < (s = *j)) {
1694 assert(T[s] == c1);
1695 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1696 assert(T[s - 1] <= T[s]);
1697 c0 = T[--s];
1698 *j = ~((int)c0);
1699 if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1700 if(c0 != c2) {
1701 if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1702 k = SA + BUCKET_B(c2 = c0, c1);
1703 }
1704 assert(k < j);
1705 *k-- = s;
1706 } else if(s != 0) {
1707 *j = ~s;
1708#ifndef NDEBUG
1709 } else {
1710 assert(T[s] == c1);
1711#endif
1712 }
1713 }
1714 }
1715 }
1716
1717 /* Construct the BWTed string by using
1718 the sorted order of type B suffixes. */
1719 k = SA + BUCKET_A(c2 = T[n - 1]);
1720 *k++ = (T[n - 2] < c2) ? ~((int)T[n - 2]) : (n - 1);
1721 /* Scan the suffix array from left to right. */
1722 for(i = SA, j = SA + n, orig = SA; i < j; ++i) {
1723 if(0 < (s = *i)) {
1724 assert(T[s - 1] >= T[s]);
1725 c0 = T[--s];
1726 *i = c0;
1727 if((0 < s) && (T[s - 1] < c0)) { s = ~((int)T[s - 1]); }
1728 if(c0 != c2) {
1729 BUCKET_A(c2) = k - SA;
1730 k = SA + BUCKET_A(c2 = c0);
1731 }
1732 assert(i < k);
1733 *k++ = s;
1734 } else if(s != 0) {
1735 *i = ~s;
1736 } else {
1737 orig = i;
1738 }
1739 }
1740
1741 return orig - SA;
1742}
1743
1744/* Constructs the burrows-wheeler transformed string directly
1745 by using the sorted order of type B* suffixes. */
1746static
1747int
1748construct_BWT_indexes(const unsigned char *T, int *SA,
1749 int *bucket_A, int *bucket_B,
1750 int n, int m,
1751 unsigned char * num_indexes, int * indexes) {
1752 int *i, *j, *k, *orig;
1753 int s;
1754 int c0, c1, c2;
1755
1756 int mod = n / 8;
1757 {
1758 mod |= mod >> 1; mod |= mod >> 2;
1759 mod |= mod >> 4; mod |= mod >> 8;
1760 mod |= mod >> 16; mod >>= 1;
1761
1762 *num_indexes = (unsigned char)((n - 1) / (mod + 1));
1763 }
1764
1765 if(0 < m) {
1766 /* Construct the sorted order of type B suffixes by using
1767 the sorted order of type B* suffixes. */
1768 for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1769 /* Scan the suffix array from right to left. */
1770 for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1771 j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1772 i <= j;
1773 --j) {
1774 if(0 < (s = *j)) {
1775 assert(T[s] == c1);
1776 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1777 assert(T[s - 1] <= T[s]);
1778
1779 if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = j - SA;
1780
1781 c0 = T[--s];
1782 *j = ~((int)c0);
1783 if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1784 if(c0 != c2) {
1785 if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1786 k = SA + BUCKET_B(c2 = c0, c1);
1787 }
1788 assert(k < j);
1789 *k-- = s;
1790 } else if(s != 0) {
1791 *j = ~s;
1792#ifndef NDEBUG
1793 } else {
1794 assert(T[s] == c1);
1795#endif
1796 }
1797 }
1798 }
1799 }
1800
1801 /* Construct the BWTed string by using
1802 the sorted order of type B suffixes. */
1803 k = SA + BUCKET_A(c2 = T[n - 1]);
1804 if (T[n - 2] < c2) {
1805 if (((n - 1) & mod) == 0) indexes[(n - 1) / (mod + 1) - 1] = k - SA;
1806 *k++ = ~((int)T[n - 2]);
1807 }
1808 else {
1809 *k++ = n - 1;
1810 }
1811
1812 /* Scan the suffix array from left to right. */
1813 for(i = SA, j = SA + n, orig = SA; i < j; ++i) {
1814 if(0 < (s = *i)) {
1815 assert(T[s - 1] >= T[s]);
1816
1817 if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = i - SA;
1818
1819 c0 = T[--s];
1820 *i = c0;
1821 if(c0 != c2) {
1822 BUCKET_A(c2) = k - SA;
1823 k = SA + BUCKET_A(c2 = c0);
1824 }
1825 assert(i < k);
1826 if((0 < s) && (T[s - 1] < c0)) {
1827 if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = k - SA;
1828 *k++ = ~((int)T[s - 1]);
1829 } else
1830 *k++ = s;
1831 } else if(s != 0) {
1832 *i = ~s;
1833 } else {
1834 orig = i;
1835 }
1836 }
1837
1838 return orig - SA;
1839}
1840
1841
1842/*---------------------------------------------------------------------------*/
1843
1844/*- Function -*/
1845
1846int
1847divsufsort(const unsigned char *T, int *SA, int n, int openMP) {
1848 int *bucket_A, *bucket_B;
1849 int m;
1850 int err = 0;
1851
1852 /* Check arguments. */
1853 if((T == NULL) || (SA == NULL) || (n < 0)) { return -1; }
1854 else if(n == 0) { return 0; }
1855 else if(n == 1) { SA[0] = 0; return 0; }
1856 else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; }
1857
1858 bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int));
1859 bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int));
1860
1861 /* Suffixsort. */
1862 if((bucket_A != NULL) && (bucket_B != NULL)) {
1863 m = sort_typeBstar(T, SA, bucket_A, bucket_B, n, openMP);
1864 construct_SA(T, SA, bucket_A, bucket_B, n, m);
1865 } else {
1866 err = -2;
1867 }
1868
1869 free(bucket_B);
1870 free(bucket_A);
1871
1872 return err;
1873}
1874
1875int
1876divbwt(const unsigned char *T, unsigned char *U, int *A, int n, unsigned char * num_indexes, int * indexes, int openMP) {
1877 int *B;
1878 int *bucket_A, *bucket_B;
1879 int m, pidx, i;
1880
1881 /* Check arguments. */
1882 if((T == NULL) || (U == NULL) || (n < 0)) { return -1; }
1883 else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; }
1884
1885 if((B = A) == NULL) { B = (int *)malloc((size_t)(n + 1) * sizeof(int)); }
1886 bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int));
1887 bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int));
1888
1889 /* Burrows-Wheeler Transform. */
1890 if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) {
1891 m = sort_typeBstar(T, B, bucket_A, bucket_B, n, openMP);
1892
1893 if (num_indexes == NULL || indexes == NULL) {
1894 pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m);
1895 } else {
1896 pidx = construct_BWT_indexes(T, B, bucket_A, bucket_B, n, m, num_indexes, indexes);
1897 }
1898
1899 /* Copy to output string. */
1900 U[0] = T[n - 1];
1901 for(i = 0; i < pidx; ++i) { U[i + 1] = (unsigned char)B[i]; }
1902 for(i += 1; i < n; ++i) { U[i] = (unsigned char)B[i]; }
1903 pidx += 1;
1904 } else {
1905 pidx = -2;
1906 }
1907
1908 free(bucket_B);
1909 free(bucket_A);
1910 if(A == NULL) { free(B); }
1911
1912 return pidx;
1913}
1914