2 * divsufsort.c for libdivsufsort-lite
3 * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
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
14 * The above copyright notice and this permission notice shall be
15 * included in all copies or substantial portions of the Software.
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.
29 #include "wimlib/divsufsort.h"
33 #if defined(ALPHABET_SIZE) && (ALPHABET_SIZE < 1)
36 #if !defined(ALPHABET_SIZE)
37 # define ALPHABET_SIZE (256)
39 #define BUCKET_A_SIZE (ALPHABET_SIZE)
40 #define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
41 #if defined(SS_INSERTIONSORT_THRESHOLD)
42 # if SS_INSERTIONSORT_THRESHOLD < 1
43 # undef SS_INSERTIONSORT_THRESHOLD
44 # define SS_INSERTIONSORT_THRESHOLD (1)
47 # define SS_INSERTIONSORT_THRESHOLD (8)
49 #if defined(SS_BLOCKSIZE)
52 # define SS_BLOCKSIZE (0)
53 # elif 32768 <= SS_BLOCKSIZE
55 # define SS_BLOCKSIZE (32767)
58 # define SS_BLOCKSIZE (1024)
60 /* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */
62 # define SS_MISORT_STACKSIZE (96)
63 #elif SS_BLOCKSIZE <= 4096
64 # define SS_MISORT_STACKSIZE (16)
66 # define SS_MISORT_STACKSIZE (24)
68 #define SS_SMERGE_STACKSIZE (32)
69 #define TR_INSERTIONSORT_THRESHOLD (8)
70 #define TR_STACKSIZE (64)
75 # define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0)
78 # define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b))
81 # define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b))
83 #define STACK_PUSH(_a, _b, _c, _d)\
85 assert(ssize < STACK_SIZE);\
86 stack[ssize].a = (_a), stack[ssize].b = (_b),\
87 stack[ssize].c = (_c), stack[ssize++].d = (_d);\
89 #define STACK_PUSH5(_a, _b, _c, _d, _e)\
91 assert(ssize < STACK_SIZE);\
92 stack[ssize].a = (_a), stack[ssize].b = (_b),\
93 stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\
95 #define STACK_POP(_a, _b, _c, _d)\
98 if(ssize == 0) { return; }\
99 (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
100 (_c) = stack[ssize].c, (_d) = stack[ssize].d;\
102 #define STACK_POP5(_a, _b, _c, _d, _e)\
105 if(ssize == 0) { return; }\
106 (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
107 (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\
109 #define BUCKET_A(_c0) bucket_A[(_c0)]
110 #if ALPHABET_SIZE == 256
111 #define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
112 #define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
114 #define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)])
115 #define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)])
119 /*- Private Functions -*/
121 static const int lg_table[256]= {
122 -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,
123 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,
124 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,
125 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,
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,
127 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
128 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
129 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
132 #if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
137 #if SS_BLOCKSIZE == 0
138 return (n & 0xffff0000) ?
140 24 + lg_table[(n >> 24) & 0xff] :
141 16 + lg_table[(n >> 16) & 0xff]) :
143 8 + lg_table[(n >> 8) & 0xff] :
144 0 + lg_table[(n >> 0) & 0xff]);
145 #elif SS_BLOCKSIZE < 256
148 return (n & 0xff00) ?
149 8 + lg_table[(n >> 8) & 0xff] :
150 0 + lg_table[(n >> 0) & 0xff];
154 #endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
156 #if SS_BLOCKSIZE != 0
158 static const int sqq_table[256] = {
159 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61,
160 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89,
161 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109,
162 110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
163 128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
164 143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155,
165 156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
166 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180,
167 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191,
168 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201,
169 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,
171 221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230,
172 230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
173 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247,
174 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255
182 if(x >= (SS_BLOCKSIZE * SS_BLOCKSIZE)) { return SS_BLOCKSIZE; }
183 e = (x & 0xffff0000) ?
185 24 + lg_table[(x >> 24) & 0xff] :
186 16 + lg_table[(x >> 16) & 0xff]) :
188 8 + lg_table[(x >> 8) & 0xff] :
189 0 + lg_table[(x >> 0) & 0xff]);
192 y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
193 if(e >= 24) { y = (y + 1 + x / y) >> 1; }
194 y = (y + 1 + x / y) >> 1;
196 y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
198 return sqq_table[x] >> 4;
201 return (x < (y * y)) ? y - 1 : y;
204 #endif /* SS_BLOCKSIZE != 0 */
207 /*---------------------------------------------------------------------------*/
209 /* Compares two suffixes. */
212 ss_compare(const unsigned char *T,
213 const int *p1, const int *p2,
215 const unsigned char *U1, *U2, *U1n, *U2n;
217 for(U1 = T + depth + *p1,
218 U2 = T + depth + *p2,
219 U1n = T + *(p1 + 1) + 2,
220 U2n = T + *(p2 + 1) + 2;
221 (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
226 (U2 < U2n ? *U1 - *U2 : 1) :
231 /*---------------------------------------------------------------------------*/
233 #if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
235 /* Insertionsort for small size groups */
238 ss_insertionsort(const unsigned char *T, const int *PA,
239 int *first, int *last, int depth) {
244 for(i = last - 2; first <= i; --i) {
245 for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) {
246 do { *(j - 1) = *j; } while((++j < last) && (*j < 0));
247 if(last <= j) { break; }
249 if(r == 0) { *j = ~*j; }
254 #endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */
257 /*---------------------------------------------------------------------------*/
259 #if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
263 ss_fixdown(const unsigned char *Td, const int *PA,
264 int *SA, int i, int size) {
269 for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
270 d = Td[PA[SA[k = j++]]];
271 if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; }
272 if(d <= c) { break; }
277 /* Simple top-down heapsort. */
280 ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) {
285 if((size % 2) == 0) {
287 if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); }
290 for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); }
291 if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); }
292 for(i = m - 1; 0 < i; --i) {
293 t = SA[0], SA[0] = SA[i];
294 ss_fixdown(Td, PA, SA, 0, i);
300 /*---------------------------------------------------------------------------*/
302 /* Returns the median of three elements. */
305 ss_median3(const unsigned char *Td, const int *PA,
306 int *v1, int *v2, int *v3) {
308 if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); }
309 if(Td[PA[*v2]] > Td[PA[*v3]]) {
310 if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; }
316 /* Returns the median of five elements. */
319 ss_median5(const unsigned char *Td, const int *PA,
320 int *v1, int *v2, int *v3, int *v4, int *v5) {
322 if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); }
323 if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); }
324 if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); }
325 if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); }
326 if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); }
327 if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; }
331 /* Returns the pivot element. */
334 ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) {
339 middle = first + t / 2;
343 return ss_median3(Td, PA, first, middle, last - 1);
346 return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
350 first = ss_median3(Td, PA, first, first + t, first + (t << 1));
351 middle = ss_median3(Td, PA, middle - t, middle, middle + t);
352 last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
353 return ss_median3(Td, PA, first, middle, last);
357 /*---------------------------------------------------------------------------*/
359 /* Binary partition for substrings. */
362 ss_partition(const int *PA,
363 int *first, int *last, int depth) {
366 for(a = first - 1, b = last;;) {
367 for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
368 for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { }
369 if(b <= a) { break; }
374 if(first < a) { *first = ~*first; }
378 /* Multikey introsort for medium size groups. */
381 ss_mintrosort(const unsigned char *T, const int *PA,
382 int *first, int *last,
384 #define STACK_SIZE SS_MISORT_STACKSIZE
385 struct { int *a, *b, c; int d; } stack[STACK_SIZE];
386 const unsigned char *Td;
387 int *a, *b, *c, *d, *e, *f;
393 for(ssize = 0, limit = ss_ilg(last - first);;) {
395 if((last - first) <= SS_INSERTIONSORT_THRESHOLD) {
396 #if 1 < SS_INSERTIONSORT_THRESHOLD
397 if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); }
399 STACK_POP(first, last, depth, limit);
404 if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); }
406 for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) {
407 if((x = Td[PA[*a]]) != v) {
408 if(1 < (a - first)) { break; }
413 if(Td[PA[*first] - 1] < v) {
414 first = ss_partition(PA, first, a, depth);
416 if((a - first) <= (last - a)) {
417 if(1 < (a - first)) {
418 STACK_PUSH(a, last, depth, -1);
419 last = a, depth += 1, limit = ss_ilg(a - first);
421 first = a, limit = -1;
425 STACK_PUSH(first, a, depth + 1, ss_ilg(a - first));
426 first = a, limit = -1;
428 last = a, depth += 1, limit = ss_ilg(a - first);
435 a = ss_pivot(Td, PA, first, last);
440 for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { }
441 if(((a = b) < last) && (x < v)) {
442 for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) {
443 if(x == v) { SWAP(*b, *a); ++a; }
446 for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { }
447 if((b < (d = c)) && (x > v)) {
448 for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
449 if(x == v) { SWAP(*c, *d); --d; }
454 for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) {
455 if(x == v) { SWAP(*b, *a); ++a; }
457 for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
458 if(x == v) { SWAP(*c, *d); --d; }
465 if((s = a - first) > (t = b - a)) { s = t; }
466 for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
467 if((s = d - c) > (t = last - d - 1)) { s = t; }
468 for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
470 a = first + (b - a), c = last - (d - c);
471 b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth);
473 if((a - first) <= (last - c)) {
474 if((last - c) <= (c - b)) {
475 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
476 STACK_PUSH(c, last, depth, limit);
478 } else if((a - first) <= (c - b)) {
479 STACK_PUSH(c, last, depth, limit);
480 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
483 STACK_PUSH(c, last, depth, limit);
484 STACK_PUSH(first, a, depth, limit);
485 first = b, last = c, depth += 1, limit = ss_ilg(c - b);
488 if((a - first) <= (c - b)) {
489 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
490 STACK_PUSH(first, a, depth, limit);
492 } else if((last - c) <= (c - b)) {
493 STACK_PUSH(first, a, depth, limit);
494 STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
497 STACK_PUSH(first, a, depth, limit);
498 STACK_PUSH(c, last, depth, limit);
499 first = b, last = c, depth += 1, limit = ss_ilg(c - b);
504 if(Td[PA[*first] - 1] < v) {
505 first = ss_partition(PA, first, last, depth);
506 limit = ss_ilg(last - first);
514 #endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
517 /*---------------------------------------------------------------------------*/
519 #if SS_BLOCKSIZE != 0
523 ss_blockswap(int *a, int *b, int n) {
525 for(; 0 < n; --n, ++a, ++b) {
526 t = *a, *a = *b, *b = t;
532 ss_rotate(int *first, int *middle, int *last) {
535 l = middle - first, r = last - middle;
536 for(; (0 < l) && (0 < r);) {
537 if(l == r) { ss_blockswap(first, middle, l); break; }
539 a = last - 1, b = middle - 1;
542 *a-- = *b, *b-- = *a;
546 if((r -= l + 1) <= l) { break; }
547 a -= 1, b = middle - 1;
552 a = first, b = middle;
555 *a++ = *b, *b++ = *a;
559 if((l -= r + 1) <= r) { break; }
569 /*---------------------------------------------------------------------------*/
573 ss_inplacemerge(const unsigned char *T, const int *PA,
574 int *first, int *middle, int *last,
583 if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); }
584 else { x = 0; p = PA + *(last - 1); }
585 for(a = first, len = middle - first, half = len >> 1, r = -1;
587 len = half, half >>= 1) {
589 q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
592 half -= (len & 1) ^ 1;
598 if(r == 0) { *a = ~*a; }
599 ss_rotate(a, middle, last);
602 if(first == middle) { break; }
605 if(x != 0) { while(*--last < 0) { } }
606 if(middle == last) { break; }
611 /*---------------------------------------------------------------------------*/
613 /* Merge-forward with internal buffer. */
616 ss_mergeforward(const unsigned char *T, const int *PA,
617 int *first, int *middle, int *last,
618 int *buf, int depth) {
619 int *a, *b, *c, *bufend;
623 bufend = buf + (middle - first) - 1;
624 ss_blockswap(buf, first, middle - first);
626 for(t = *(a = first), b = buf, c = middle;;) {
627 r = ss_compare(T, PA + *b, PA + *c, depth);
631 if(bufend <= b) { *bufend = t; return; }
636 *a++ = *c, *c++ = *a;
638 while(b < bufend) { *a++ = *b, *b++ = *a; }
647 if(bufend <= b) { *bufend = t; return; }
652 *a++ = *c, *c++ = *a;
654 while(b < bufend) { *a++ = *b, *b++ = *a; }
663 /* Merge-backward with internal buffer. */
666 ss_mergebackward(const unsigned char *T, const int *PA,
667 int *first, int *middle, int *last,
668 int *buf, int depth) {
670 int *a, *b, *c, *bufend;
675 bufend = buf + (last - middle) - 1;
676 ss_blockswap(buf, middle, last - middle);
679 if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; }
680 else { p1 = PA + *bufend; }
681 if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; }
682 else { p2 = PA + *(middle - 1); }
683 for(t = *(a = last - 1), b = bufend, c = middle - 1;;) {
684 r = ss_compare(T, p1, p2, depth);
686 if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
688 if(b <= buf) { *buf = t; break; }
690 if(*b < 0) { p1 = PA + ~*b; x |= 1; }
691 else { p1 = PA + *b; }
693 if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
694 *a-- = *c, *c-- = *a;
696 while(buf < b) { *a-- = *b, *b-- = *a; }
700 if(*c < 0) { p2 = PA + ~*c; x |= 2; }
701 else { p2 = PA + *c; }
703 if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
705 if(b <= buf) { *buf = t; break; }
707 if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
708 *a-- = *c, *c-- = *a;
710 while(buf < b) { *a-- = *b, *b-- = *a; }
714 if(*b < 0) { p1 = PA + ~*b; x |= 1; }
715 else { p1 = PA + *b; }
716 if(*c < 0) { p2 = PA + ~*c; x |= 2; }
717 else { p2 = PA + *c; }
722 /* D&C based merge. */
725 ss_swapmerge(const unsigned char *T, const int *PA,
726 int *first, int *middle, int *last,
727 int *buf, int bufsize, int depth) {
728 #define STACK_SIZE SS_SMERGE_STACKSIZE
729 #define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
730 #define MERGE_CHECK(a, b, c)\
733 (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\
736 if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\
740 struct { int *a, *b, *c; int d; } stack[STACK_SIZE];
741 int *l, *r, *lm, *rm;
746 for(check = 0, ssize = 0;;) {
747 if((last - middle) <= bufsize) {
748 if((first < middle) && (middle < last)) {
749 ss_mergebackward(T, PA, first, middle, last, buf, depth);
751 MERGE_CHECK(first, last, check);
752 STACK_POP(first, middle, last, check);
756 if((middle - first) <= bufsize) {
758 ss_mergeforward(T, PA, first, middle, last, buf, depth);
760 MERGE_CHECK(first, last, check);
761 STACK_POP(first, middle, last, check);
765 for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1;
767 len = half, half >>= 1) {
768 if(ss_compare(T, PA + GETIDX(*(middle + m + half)),
769 PA + GETIDX(*(middle - m - half - 1)), depth) < 0) {
771 half -= (len & 1) ^ 1;
776 lm = middle - m, rm = middle + m;
777 ss_blockswap(lm, middle, m);
778 l = r = middle, next = 0;
782 if(first < lm) { for(; *--l < 0;) { } next |= 4; }
784 } else if(first < lm) {
785 for(; *r < 0; ++r) { }
790 if((l - first) <= (last - r)) {
791 STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
792 middle = lm, last = l, check = (check & 3) | (next & 4);
794 if((next & 2) && (r == middle)) { next ^= 6; }
795 STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
796 first = r, middle = rm, check = (next & 3) | (check & 4);
799 if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) {
802 MERGE_CHECK(first, last, check);
803 STACK_POP(first, middle, last, check);
809 #endif /* SS_BLOCKSIZE != 0 */
812 /*---------------------------------------------------------------------------*/
817 sssort(const unsigned char *T, const int *PA,
818 int *first, int *last,
819 int *buf, int bufsize,
820 int depth, int n, int lastsuffix) {
822 #if SS_BLOCKSIZE != 0
823 int *b, *middle, *curbuf;
824 int j, k, curbufsize, limit;
828 if(lastsuffix != 0) { ++first; }
830 #if SS_BLOCKSIZE == 0
831 ss_mintrosort(T, PA, first, last, depth);
833 if((bufsize < SS_BLOCKSIZE) &&
834 (bufsize < (last - first)) &&
835 (bufsize < (limit = ss_isqrt(last - first)))) {
836 if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; }
837 buf = middle = last - limit, bufsize = limit;
839 middle = last, limit = 0;
841 for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) {
842 #if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
843 ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth);
844 #elif 1 < SS_BLOCKSIZE
845 ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth);
847 curbufsize = last - (a + SS_BLOCKSIZE);
848 curbuf = a + SS_BLOCKSIZE;
849 if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
850 for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) {
851 ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
854 #if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
855 ss_mintrosort(T, PA, a, middle, depth);
856 #elif 1 < SS_BLOCKSIZE
857 ss_insertionsort(T, PA, a, middle, depth);
859 for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) {
861 ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
866 #if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
867 ss_mintrosort(T, PA, middle, last, depth);
868 #elif 1 < SS_BLOCKSIZE
869 ss_insertionsort(T, PA, middle, last, depth);
871 ss_inplacemerge(T, PA, first, middle, last, depth);
875 if(lastsuffix != 0) {
876 /* Insert last type B* suffix. */
877 int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
878 for(a = first, i = *(first - 1);
879 (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth)));
888 /*---------------------------------------------------------------------------*/
893 return (n & 0xffff0000) ?
895 24 + lg_table[(n >> 24) & 0xff] :
896 16 + lg_table[(n >> 16) & 0xff]) :
898 8 + lg_table[(n >> 8) & 0xff] :
899 0 + lg_table[(n >> 0) & 0xff]);
903 /*---------------------------------------------------------------------------*/
905 /* Simple insertionsort for small size groups. */
908 tr_insertionsort(const int *ISAd, int *first, int *last) {
912 for(a = first + 1; a < last; ++a) {
913 for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) {
914 do { *(b + 1) = *b; } while((first <= --b) && (*b < 0));
915 if(b < first) { break; }
917 if(r == 0) { *b = ~*b; }
923 /*---------------------------------------------------------------------------*/
927 tr_fixdown(const int *ISAd, int *SA, int i, int size) {
932 for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
933 d = ISAd[SA[k = j++]];
934 if(d < (e = ISAd[SA[j]])) { k = j; d = e; }
935 if(d <= c) { break; }
940 /* Simple top-down heapsort. */
943 tr_heapsort(const int *ISAd, int *SA, int size) {
948 if((size % 2) == 0) {
950 if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); }
953 for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); }
954 if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); }
955 for(i = m - 1; 0 < i; --i) {
956 t = SA[0], SA[0] = SA[i];
957 tr_fixdown(ISAd, SA, 0, i);
963 /*---------------------------------------------------------------------------*/
965 /* Returns the median of three elements. */
968 tr_median3(const int *ISAd, int *v1, int *v2, int *v3) {
970 if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); }
971 if(ISAd[*v2] > ISAd[*v3]) {
972 if(ISAd[*v1] > ISAd[*v3]) { return v1; }
978 /* Returns the median of five elements. */
981 tr_median5(const int *ISAd,
982 int *v1, int *v2, int *v3, int *v4, int *v5) {
984 if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); }
985 if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); }
986 if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4); SWAP(v3, v5); }
987 if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3); }
988 if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4); SWAP(v3, v5); }
989 if(ISAd[*v3] > ISAd[*v4]) { return v4; }
993 /* Returns the pivot element. */
996 tr_pivot(const int *ISAd, int *first, int *last) {
1001 middle = first + t / 2;
1005 return tr_median3(ISAd, first, middle, last - 1);
1008 return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1012 first = tr_median3(ISAd, first, first + t, first + (t << 1));
1013 middle = tr_median3(ISAd, middle - t, middle, middle + t);
1014 last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
1015 return tr_median3(ISAd, first, middle, last);
1019 /*---------------------------------------------------------------------------*/
1021 typedef struct _trbudget_t trbudget_t;
1022 struct _trbudget_t {
1031 trbudget_init(trbudget_t *budget, int chance, int incval) {
1032 budget->chance = chance;
1033 budget->remain = budget->incval = incval;
1038 trbudget_check(trbudget_t *budget, int size) {
1039 if(size <= budget->remain) { budget->remain -= size; return 1; }
1040 if(budget->chance == 0) { budget->count += size; return 0; }
1041 budget->remain += budget->incval - size;
1042 budget->chance -= 1;
1047 /*---------------------------------------------------------------------------*/
1051 tr_partition(const int *ISAd,
1052 int *first, int *middle, int *last,
1053 int **pa, int **pb, int v) {
1054 int *a, *b, *c, *d, *e, *f;
1058 for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { }
1059 if(((a = b) < last) && (x < v)) {
1060 for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
1061 if(x == v) { SWAP(*b, *a); ++a; }
1064 for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { }
1065 if((b < (d = c)) && (x > v)) {
1066 for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1067 if(x == v) { SWAP(*c, *d); --d; }
1072 for(; (++b < c) && ((x = ISAd[*b]) <= v);) {
1073 if(x == v) { SWAP(*b, *a); ++a; }
1075 for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1076 if(x == v) { SWAP(*c, *d); --d; }
1082 if((s = a - first) > (t = b - a)) { s = t; }
1083 for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1084 if((s = d - c) > (t = last - d - 1)) { s = t; }
1085 for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1086 first += (b - a), last -= (d - c);
1088 *pa = first, *pb = last;
1093 tr_copy(int *ISA, const int *SA,
1094 int *first, int *a, int *b, int *last,
1096 /* sort suffixes of middle partition
1097 by using sorted order of suffixes of left and right partition. */
1102 for(c = first, d = a - 1; c <= d; ++c) {
1103 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1108 for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1109 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1118 tr_partialcopy(int *ISA, const int *SA,
1119 int *first, int *a, int *b, int *last,
1123 int rank, lastrank, newrank = -1;
1127 for(c = first, d = a - 1; c <= d; ++c) {
1128 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1130 rank = ISA[s + depth];
1131 if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1137 for(e = d; first <= e; --e) {
1139 if(lastrank != rank) { lastrank = rank; newrank = e - SA; }
1140 if(newrank != rank) { ISA[*e] = newrank; }
1144 for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1145 if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1147 rank = ISA[s + depth];
1148 if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1156 tr_introsort(int *ISA, const int *ISAd,
1157 int *SA, int *first, int *last,
1158 trbudget_t *budget) {
1159 #define STACK_SIZE TR_STACKSIZE
1160 struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE];
1164 int incr = ISAd - ISA;
1166 int ssize, trlink = -1;
1168 for(ssize = 0, limit = tr_ilg(last - first);;) {
1172 /* tandem repeat partition */
1173 tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1);
1177 for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1180 for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
1185 STACK_PUSH5(NULL, a, b, 0, 0);
1186 STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1189 if((a - first) <= (last - b)) {
1190 if(1 < (a - first)) {
1191 STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink);
1192 last = a, limit = tr_ilg(a - first);
1193 } else if(1 < (last - b)) {
1194 first = b, limit = tr_ilg(last - b);
1196 STACK_POP5(ISAd, first, last, limit, trlink);
1199 if(1 < (last - b)) {
1200 STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink);
1201 first = b, limit = tr_ilg(last - b);
1202 } else if(1 < (a - first)) {
1203 last = a, limit = tr_ilg(a - first);
1205 STACK_POP5(ISAd, first, last, limit, trlink);
1208 } else if(limit == -2) {
1209 /* tandem repeat copy */
1210 a = stack[--ssize].b, b = stack[ssize].c;
1211 if(stack[ssize].d == 0) {
1212 tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
1214 if(0 <= trlink) { stack[trlink].d = -1; }
1215 tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
1217 STACK_POP5(ISAd, first, last, limit, trlink);
1219 /* sorted partition */
1222 do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a));
1226 a = first; do { *a = ~*a; } while(*++a < 0);
1227 next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1;
1228 if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } }
1231 if(trbudget_check(budget, a - first)) {
1232 if((a - first) <= (last - a)) {
1233 STACK_PUSH5(ISAd, a, last, -3, trlink);
1234 ISAd += incr, last = a, limit = next;
1236 if(1 < (last - a)) {
1237 STACK_PUSH5(ISAd + incr, first, a, next, trlink);
1238 first = a, limit = -3;
1240 ISAd += incr, last = a, limit = next;
1244 if(0 <= trlink) { stack[trlink].d = -1; }
1245 if(1 < (last - a)) {
1246 first = a, limit = -3;
1248 STACK_POP5(ISAd, first, last, limit, trlink);
1252 STACK_POP5(ISAd, first, last, limit, trlink);
1258 if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
1259 tr_insertionsort(ISAd, first, last);
1265 tr_heapsort(ISAd, first, last - first);
1266 for(a = last - 1; first < a; a = b) {
1267 for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
1274 a = tr_pivot(ISAd, first, last);
1279 tr_partition(ISAd, first, first + 1, last, &a, &b, v);
1280 if((last - first) != (b - a)) {
1281 next = (ISA[*a] != v) ? tr_ilg(b - a) : -1;
1284 for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1285 if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } }
1288 if((1 < (b - a)) && (trbudget_check(budget, b - a))) {
1289 if((a - first) <= (last - b)) {
1290 if((last - b) <= (b - a)) {
1291 if(1 < (a - first)) {
1292 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1293 STACK_PUSH5(ISAd, b, last, limit, trlink);
1295 } else if(1 < (last - b)) {
1296 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1299 ISAd += incr, first = a, last = b, limit = next;
1301 } else if((a - first) <= (b - a)) {
1302 if(1 < (a - first)) {
1303 STACK_PUSH5(ISAd, b, last, limit, trlink);
1304 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1307 STACK_PUSH5(ISAd, b, last, limit, trlink);
1308 ISAd += incr, first = a, last = b, limit = next;
1311 STACK_PUSH5(ISAd, b, last, limit, trlink);
1312 STACK_PUSH5(ISAd, first, a, limit, trlink);
1313 ISAd += incr, first = a, last = b, limit = next;
1316 if((a - first) <= (b - a)) {
1317 if(1 < (last - b)) {
1318 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1319 STACK_PUSH5(ISAd, first, a, limit, trlink);
1321 } else if(1 < (a - first)) {
1322 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1325 ISAd += incr, first = a, last = b, limit = next;
1327 } else if((last - b) <= (b - a)) {
1328 if(1 < (last - b)) {
1329 STACK_PUSH5(ISAd, first, a, limit, trlink);
1330 STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1333 STACK_PUSH5(ISAd, first, a, limit, trlink);
1334 ISAd += incr, first = a, last = b, limit = next;
1337 STACK_PUSH5(ISAd, first, a, limit, trlink);
1338 STACK_PUSH5(ISAd, b, last, limit, trlink);
1339 ISAd += incr, first = a, last = b, limit = next;
1343 if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
1344 if((a - first) <= (last - b)) {
1345 if(1 < (a - first)) {
1346 STACK_PUSH5(ISAd, b, last, limit, trlink);
1348 } else if(1 < (last - b)) {
1351 STACK_POP5(ISAd, first, last, limit, trlink);
1354 if(1 < (last - b)) {
1355 STACK_PUSH5(ISAd, first, a, limit, trlink);
1357 } else if(1 < (a - first)) {
1360 STACK_POP5(ISAd, first, last, limit, trlink);
1365 if(trbudget_check(budget, last - first)) {
1366 limit = tr_ilg(last - first), ISAd += incr;
1368 if(0 <= trlink) { stack[trlink].d = -1; }
1369 STACK_POP5(ISAd, first, last, limit, trlink);
1378 /*---------------------------------------------------------------------------*/
1380 /* Tandem repeat sort */
1383 trsort(int *ISA, int *SA, int n, int depth) {
1387 int t, skip, unsorted;
1389 trbudget_init(&budget, tr_ilg(n) * 2 / 3, n);
1390 /* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */
1391 for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) {
1396 if((t = *first) < 0) { first -= t; skip += t; }
1398 if(skip != 0) { *(first + skip) = skip; skip = 0; }
1399 last = SA + ISA[t] + 1;
1400 if(1 < (last - first)) {
1402 tr_introsort(ISA, ISAd, SA, first, last, &budget);
1403 if(budget.count != 0) { unsorted += budget.count; }
1404 else { skip = first - last; }
1405 } else if((last - first) == 1) {
1410 } while(first < (SA + n));
1411 if(skip != 0) { *(first + skip) = skip; }
1412 if(unsorted == 0) { break; }
1417 /*---------------------------------------------------------------------------*/
1419 /* Sorts suffixes of type B*. */
1422 sort_typeBstar(const unsigned char *T, int *SA,
1423 int *bucket_A, int *bucket_B,
1425 int *PAb, *ISAb, *buf;
1426 int i, j, k, t, m, bufsize;
1429 /* Initialize bucket arrays. */
1430 for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; }
1431 for(i = 0; i < BUCKET_B_SIZE; ++i) { bucket_B[i] = 0; }
1433 /* Count the number of occurrences of the first one or two characters of each
1434 type A, B and B* suffix. Moreover, store the beginning position of all
1435 type B* suffixes into the array SA. */
1436 for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) {
1437 /* type A suffix. */
1438 do { ++BUCKET_A(c1 = c0); } while((0 <= --i) && ((c0 = T[i]) >= c1));
1440 /* type B* suffix. */
1441 ++BUCKET_BSTAR(c0, c1);
1443 /* type B suffix. */
1444 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) {
1452 A type B* suffix is lexicographically smaller than a type B suffix that
1453 begins with the same first two characters.
1456 /* Calculate the index of start/end point of each bucket. */
1457 for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0) {
1458 t = i + BUCKET_A(c0);
1459 BUCKET_A(c0) = i + j; /* start point */
1460 i = t + BUCKET_B(c0, c0);
1461 for(c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1) {
1462 j += BUCKET_BSTAR(c0, c1);
1463 BUCKET_BSTAR(c0, c1) = j; /* end point */
1464 i += BUCKET_B(c0, c1);
1469 /* Sort the type B* suffixes by their first two characters. */
1470 PAb = SA + n - m; ISAb = SA + m;
1471 for(i = m - 2; 0 <= i; --i) {
1472 t = PAb[i], c0 = T[t], c1 = T[t + 1];
1473 SA[--BUCKET_BSTAR(c0, c1)] = i;
1475 t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
1476 SA[--BUCKET_BSTAR(c0, c1)] = m - 1;
1478 /* Sort the type B* substrings using sssort. */
1479 buf = SA + m, bufsize = n - (2 * m);
1480 for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
1481 for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
1482 i = BUCKET_BSTAR(c0, c1);
1484 sssort(T, PAb, SA + i, SA + j,
1485 buf, bufsize, 2, n, *(SA + i) == (m - 1));
1490 /* Compute ranks of type B* substrings. */
1491 for(i = m - 1; 0 <= i; --i) {
1494 do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i]));
1496 if(i <= 0) { break; }
1499 do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0);
1503 /* Construct the inverse suffix array of type B* suffixes using trsort. */
1504 trsort(ISAb, SA, m, 1);
1506 /* Set the sorted order of tyoe B* suffixes. */
1507 for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) {
1508 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { }
1511 for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { }
1512 SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
1516 /* Calculate the index of start/end point of each bucket. */
1517 BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */
1518 for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) {
1519 i = BUCKET_A(c0 + 1) - 1;
1520 for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) {
1521 t = i - BUCKET_B(c0, c1);
1522 BUCKET_B(c0, c1) = i; /* end point */
1524 /* Move all type B* suffixes to the correct position. */
1525 for(i = t, j = BUCKET_BSTAR(c0, c1);
1527 --i, --k) { SA[i] = SA[k]; }
1529 BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */
1530 BUCKET_B(c0, c0) = i; /* end point */
1537 /* Constructs the suffix array by using the sorted order of type B* suffixes. */
1540 construct_SA(const unsigned char *T, int *SA,
1541 int *bucket_A, int *bucket_B,
1548 /* Construct the sorted order of type B suffixes by using
1549 the sorted order of type B* suffixes. */
1550 for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1551 /* Scan the suffix array from right to left. */
1552 for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1553 j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1558 assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1559 assert(T[s - 1] <= T[s]);
1562 if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1564 if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1565 k = SA + BUCKET_B(c2 = c0, c1);
1570 assert(((s == 0) && (T[s] == c1)) || (s < 0));
1577 /* Construct the suffix array by using
1578 the sorted order of type B suffixes. */
1579 k = SA + BUCKET_A(c2 = T[n - 1]);
1580 *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
1581 /* Scan the suffix array from left to right. */
1582 for(i = SA, j = SA + n; i < j; ++i) {
1584 assert(T[s - 1] >= T[s]);
1586 if((s == 0) || (T[s - 1] < c0)) { s = ~s; }
1588 BUCKET_A(c2) = k - SA;
1589 k = SA + BUCKET_A(c2 = c0);
1600 /*---------------------------------------------------------------------------*/
1604 /* XXX Modified from original: use provided temporary space instead of
1607 divsufsort(const unsigned char *T, int *SA, int n,
1608 int *bucket_A, int *bucket_B)
1627 m = sort_typeBstar(T, SA, bucket_A, bucket_B, n);
1628 construct_SA(T, SA, bucket_A, bucket_B, n, m);