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