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