wimlib.net Git - wimlib/blob - 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
 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 inline
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 inline
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 inline
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 }