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