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