CoCalc -- divsufsort.c

GitHub Repository: Kitware/CMake
Path: blob/master/Utilities/cmzstd/lib/dictBuilder/divsufsort.c
³¹⁵⁶ views
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
/*- Compiler specifics -*/
28
#ifdef __clang__
29
#pragma clang diagnostic ignored "-Wshorten-64-to-32"
30
#endif
31

32
#if defined(_MSC_VER)
33
#  pragma warning(disable : 4244)
34
#  pragma warning(disable : 4127)    /* C4127 : Condition expression is constant */
35
#endif
36

37

38
/*- Dependencies -*/
39
#include <assert.h>
40
#include <stdio.h>
41
#include <stdlib.h>
42

43
#ifdef __clang_analyzer__
44
#include <string.h>
45
#endif
46

47
#include "divsufsort.h"
48

49
/*- Constants -*/
50
#if defined(INLINE)
51
# undef INLINE
52
#endif
53
#if !defined(INLINE)
54
# define INLINE __inline
55
#endif
56
#if defined(ALPHABET_SIZE) && (ALPHABET_SIZE < 1)
57
# undef ALPHABET_SIZE
58
#endif
59
#if !defined(ALPHABET_SIZE)
60
# define ALPHABET_SIZE (256)
61
#endif
62
#define BUCKET_A_SIZE (ALPHABET_SIZE)
63
#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
64
#if defined(SS_INSERTIONSORT_THRESHOLD)
65
# if SS_INSERTIONSORT_THRESHOLD < 1
66
#  undef SS_INSERTIONSORT_THRESHOLD
67
#  define SS_INSERTIONSORT_THRESHOLD (1)
68
# endif
69
#else
70
# define SS_INSERTIONSORT_THRESHOLD (8)
71
#endif
72
#if defined(SS_BLOCKSIZE)
73
# if SS_BLOCKSIZE < 0
74
#  undef SS_BLOCKSIZE
75
#  define SS_BLOCKSIZE (0)
76
# elif 32768 <= SS_BLOCKSIZE
77
#  undef SS_BLOCKSIZE
78
#  define SS_BLOCKSIZE (32767)
79
# endif
80
#else
81
# define SS_BLOCKSIZE (1024)
82
#endif
83
/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */
84
#if SS_BLOCKSIZE == 0
85
# define SS_MISORT_STACKSIZE (96)
86
#elif SS_BLOCKSIZE <= 4096
87
# define SS_MISORT_STACKSIZE (16)
88
#else
89
# define SS_MISORT_STACKSIZE (24)
90
#endif
91
#define SS_SMERGE_STACKSIZE (32)
92
#define TR_INSERTIONSORT_THRESHOLD (8)
93
#define TR_STACKSIZE (64)
94

95

96
/*- Macros -*/
97
#ifndef SWAP
98
# define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0)
99
#endif /* SWAP */
100
#ifndef MIN
101
# define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b))
102
#endif /* MIN */
103
#ifndef MAX
104
# define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b))
105
#endif /* MAX */
106
#define STACK_PUSH(_a, _b, _c, _d)\
107
  do {\
108
    assert(ssize < STACK_SIZE);\
109
    stack[ssize].a = (_a), stack[ssize].b = (_b),\
110
    stack[ssize].c = (_c), stack[ssize++].d = (_d);\
111
  } while(0)
112
#define STACK_PUSH5(_a, _b, _c, _d, _e)\
113
  do {\
114
    assert(ssize < STACK_SIZE);\
115
    stack[ssize].a = (_a), stack[ssize].b = (_b),\
116
    stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\
117
  } while(0)
118
#define STACK_POP(_a, _b, _c, _d)\
119
  do {\
120
    assert(0 <= ssize);\
121
    if(ssize == 0) { return; }\
122
    (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
123
    (_c) = stack[ssize].c, (_d) = stack[ssize].d;\
124
  } while(0)
125
#define STACK_POP5(_a, _b, _c, _d, _e)\
126
  do {\
127
    assert(0 <= ssize);\
128
    if(ssize == 0) { return; }\
129
    (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
130
    (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\
131
  } while(0)
132
#define BUCKET_A(_c0) bucket_A[(_c0)]
133
#if ALPHABET_SIZE == 256
134
#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
135
#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
136
#else
137
#define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)])
138
#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)])
139
#endif
140

141

142
/*- Private Functions -*/
143

144
static const int lg_table[256]= {
145
 -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,
146
  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,
147
  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,
148
  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,
149
  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,
150
  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,
151
  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,
152
  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
153
};
154

155
#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
156

157
static INLINE
158
int
159
ss_ilg(int n) {
160
#if SS_BLOCKSIZE == 0
161
  return (n & 0xffff0000) ?
162
          ((n & 0xff000000) ?
163
            24 + lg_table[(n >> 24) & 0xff] :
164
            16 + lg_table[(n >> 16) & 0xff]) :
165
          ((n & 0x0000ff00) ?
166
             8 + lg_table[(n >>  8) & 0xff] :
167
             0 + lg_table[(n >>  0) & 0xff]);
168
#elif SS_BLOCKSIZE < 256
169
  return lg_table[n];
170
#else
171
  return (n & 0xff00) ?
172
          8 + lg_table[(n >> 8) & 0xff] :
173
          0 + lg_table[(n >> 0) & 0xff];
174
#endif
175
}
176

177
#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
178

179
#if SS_BLOCKSIZE != 0
180

181
static const int sqq_table[256] = {
182
  0,  16,  22,  27,  32,  35,  39,  42,  45,  48,  50,  53,  55,  57,  59,  61,
183
 64,  65,  67,  69,  71,  73,  75,  76,  78,  80,  81,  83,  84,  86,  87,  89,
184
 90,  91,  93,  94,  96,  97,  98,  99, 101, 102, 103, 104, 106, 107, 108, 109,
185
110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
186
128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
187
143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155,
188
156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
189
169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180,
190
181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191,
191
192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201,
192
202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211,
193
212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221,
194
221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230,
195
230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
196
239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247,
197
247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255
198
};
199

200
static INLINE
201
int
202
ss_isqrt(int x) {
203
  int y, e;
204

205
  if(x >= (SS_BLOCKSIZE * SS_BLOCKSIZE)) { return SS_BLOCKSIZE; }
206
  e = (x & 0xffff0000) ?
207
        ((x & 0xff000000) ?
208
          24 + lg_table[(x >> 24) & 0xff] :
209
          16 + lg_table[(x >> 16) & 0xff]) :
210
        ((x & 0x0000ff00) ?
211
           8 + lg_table[(x >>  8) & 0xff] :
212
           0 + lg_table[(x >>  0) & 0xff]);
213

214
  if(e >= 16) {
215
    y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
216
    if(e >= 24) { y = (y + 1 + x / y) >> 1; }
217
    y = (y + 1 + x / y) >> 1;
218
  } else if(e >= 8) {
219
    y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
220
  } else {
221
    return sqq_table[x] >> 4;
222
  }
223

224
  return (x < (y * y)) ? y - 1 : y;
225
}
226

227
#endif /* SS_BLOCKSIZE != 0 */
228

229

230
/*---------------------------------------------------------------------------*/
231

232
/* Compares two suffixes. */
233
static INLINE
234
int
235
ss_compare(const unsigned char *T,
236
           const int *p1, const int *p2,
237
           int depth) {
238
  const unsigned char *U1, *U2, *U1n, *U2n;
239

240
  for(U1 = T + depth + *p1,
241
      U2 = T + depth + *p2,
242
      U1n = T + *(p1 + 1) + 2,
243
      U2n = T + *(p2 + 1) + 2;
244
      (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
245
      ++U1, ++U2) {
246
  }
247

248
  return U1 < U1n ?
249
        (U2 < U2n ? *U1 - *U2 : 1) :
250
        (U2 < U2n ? -1 : 0);
251
}
252

253

254
/*---------------------------------------------------------------------------*/
255

256
#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
257

258
/* Insertionsort for small size groups */
259
static
260
void
261
ss_insertionsort(const unsigned char *T, const int *PA,
262
                 int *first, int *last, int depth) {
263
  int *i, *j;
264
  int t;
265
  int r;
266

267
  for(i = last - 2; first <= i; --i) {
268
    for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) {
269
      do { *(j - 1) = *j; } while((++j < last) && (*j < 0));
270
      if(last <= j) { break; }
271
    }
272
    if(r == 0) { *j = ~*j; }
273
    *(j - 1) = t;
274
  }
275
}
276

277
#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */
278

279

280
/*---------------------------------------------------------------------------*/
281

282
#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
283

284
static INLINE
285
void
286
ss_fixdown(const unsigned char *Td, const int *PA,
287
           int *SA, int i, int size) {
288
  int j, k;
289
  int v;
290
  int c, d, e;
291

292
  for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
293
    d = Td[PA[SA[k = j++]]];
294
    if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; }
295
    if(d <= c) { break; }
296
  }
297
  SA[i] = v;
298
}
299

300
/* Simple top-down heapsort. */
301
static
302
void
303
ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) {
304
  int i, m;
305
  int t;
306

307
  m = size;
308
  if((size % 2) == 0) {
309
    m--;
310
    if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); }
311
  }
312

313
  for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); }
314
  if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); }
315
  for(i = m - 1; 0 < i; --i) {
316
    t = SA[0], SA[0] = SA[i];
317
    ss_fixdown(Td, PA, SA, 0, i);
318
    SA[i] = t;
319
  }
320
}
321

322

323
/*---------------------------------------------------------------------------*/
324

325
/* Returns the median of three elements. */
326
static INLINE
327
int *
328
ss_median3(const unsigned char *Td, const int *PA,
329
           int *v1, int *v2, int *v3) {
330
  int *t;
331
  if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); }
332
  if(Td[PA[*v2]] > Td[PA[*v3]]) {
333
    if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; }
334
    else { return v3; }
335
  }
336
  return v2;
337
}
338

339
/* Returns the median of five elements. */
340
static INLINE
341
int *
342
ss_median5(const unsigned char *Td, const int *PA,
343
           int *v1, int *v2, int *v3, int *v4, int *v5) {
344
  int *t;
345
  if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); }
346
  if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); }
347
  if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); }
348
  if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); }
349
  if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); }
350
  if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; }
351
  return v3;
352
}
353

354
/* Returns the pivot element. */
355
static INLINE
356
int *
357
ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) {
358
  int *middle;
359
  int t;
360

361
  t = last - first;
362
  middle = first + t / 2;
363

364
  if(t <= 512) {
365
    if(t <= 32) {
366
      return ss_median3(Td, PA, first, middle, last - 1);
367
    } else {
368
      t >>= 2;
369
      return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
370
    }
371
  }
372
  t >>= 3;
373
  first  = ss_median3(Td, PA, first, first + t, first + (t << 1));
374
  middle = ss_median3(Td, PA, middle - t, middle, middle + t);
375
  last   = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
376
  return ss_median3(Td, PA, first, middle, last);
377
}
378

379

380
/*---------------------------------------------------------------------------*/
381

382
/* Binary partition for substrings. */
383
static INLINE
384
int *
385
ss_partition(const int *PA,
386
                    int *first, int *last, int depth) {
387
  int *a, *b;
388
  int t;
389
  for(a = first - 1, b = last;;) {
390
    for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
391
    for(; (a < --b) && ((PA[*b] + depth) <  (PA[*b + 1] + 1));) { }
392
    if(b <= a) { break; }
393
    t = ~*b;
394
    *b = *a;
395
    *a = t;
396
  }
397
  if(first < a) { *first = ~*first; }
398
  return a;
399
}
400

401
/* Multikey introsort for medium size groups. */
402
static
403
void
404
ss_mintrosort(const unsigned char *T, const int *PA,
405
              int *first, int *last,
406
              int depth) {
407
#define STACK_SIZE SS_MISORT_STACKSIZE
408
  struct { int *a, *b, c; int d; } stack[STACK_SIZE];
409
  const unsigned char *Td;
410
  int *a, *b, *c, *d, *e, *f;
411
  int s, t;
412
  int ssize;
413
  int limit;
414
  int v, x = 0;
415

416
  for(ssize = 0, limit = ss_ilg(last - first);;) {
417

418
    if((last - first) <= SS_INSERTIONSORT_THRESHOLD) {
419
#if 1 < SS_INSERTIONSORT_THRESHOLD
420
      if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); }
421
#endif
422
      STACK_POP(first, last, depth, limit);
423
      continue;
424
    }
425

426
    Td = T + depth;
427
    if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); }
428
    if(limit < 0) {
429
      for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) {
430
        if((x = Td[PA[*a]]) != v) {
431
          if(1 < (a - first)) { break; }
432
          v = x;
433
          first = a;
434
        }
435
      }
436
      if(Td[PA[*first] - 1] < v) {
437
        first = ss_partition(PA, first, a, depth);
438
      }
439
      if((a - first) <= (last - a)) {
440
        if(1 < (a - first)) {
441
          STACK_PUSH(a, last, depth, -1);
442
          last = a, depth += 1, limit = ss_ilg(a - first);
443
        } else {
444
          first = a, limit = -1;
445
        }
446
      } else {
447
        if(1 < (last - a)) {
448
          STACK_PUSH(first, a, depth + 1, ss_ilg(a - first));
449
          first = a, limit = -1;
450
        } else {
451
          last = a, depth += 1, limit = ss_ilg(a - first);
452
        }
453
      }
454
      continue;
455
    }
456

457
    /* choose pivot */
458
    a = ss_pivot(Td, PA, first, last);
459
    v = Td[PA[*a]];
460
    SWAP(*first, *a);
461

462
    /* partition */
463
    for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { }
464
    if(((a = b) < last) && (x < v)) {
465
      for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) {
466
        if(x == v) { SWAP(*b, *a); ++a; }
467
      }
468
    }
469
    for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { }
470
    if((b < (d = c)) && (x > v)) {
471
      for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
472
        if(x == v) { SWAP(*c, *d); --d; }
473
      }
474
    }
475
    for(; b < c;) {
476
      SWAP(*b, *c);
477
      for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) {
478
        if(x == v) { SWAP(*b, *a); ++a; }
479
      }
480
      for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
481
        if(x == v) { SWAP(*c, *d); --d; }
482
      }
483
    }
484

485
    if(a <= d) {
486
      c = b - 1;
487

488
      if((s = a - first) > (t = b - a)) { s = t; }
489
      for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
490
      if((s = d - c) > (t = last - d - 1)) { s = t; }
491
      for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
492

493
      a = first + (b - a), c = last - (d - c);
494
      b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth);
495

496
      if((a - first) <= (last - c)) {
497
        if((last - c) <= (c - b)) {
498
          STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
499
          STACK_PUSH(c, last, depth, limit);
500
          last = a;
501
        } else if((a - first) <= (c - b)) {
502
          STACK_PUSH(c, last, depth, limit);
503
          STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
504
          last = a;
505
        } else {
506
          STACK_PUSH(c, last, depth, limit);
507
          STACK_PUSH(first, a, depth, limit);
508
          first = b, last = c, depth += 1, limit = ss_ilg(c - b);
509
        }
510
      } else {
511
        if((a - first) <= (c - b)) {
512
          STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
513
          STACK_PUSH(first, a, depth, limit);
514
          first = c;
515
        } else if((last - c) <= (c - b)) {
516
          STACK_PUSH(first, a, depth, limit);
517
          STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
518
          first = c;
519
        } else {
520
          STACK_PUSH(first, a, depth, limit);
521
          STACK_PUSH(c, last, depth, limit);
522
          first = b, last = c, depth += 1, limit = ss_ilg(c - b);
523
        }
524
      }
525
    } else {
526
      limit += 1;
527
      if(Td[PA[*first] - 1] < v) {
528
        first = ss_partition(PA, first, last, depth);
529
        limit = ss_ilg(last - first);
530
      }
531
      depth += 1;
532
    }
533
  }
534
#undef STACK_SIZE
535
}
536

537
#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
538

539

540
/*---------------------------------------------------------------------------*/
541

542
#if SS_BLOCKSIZE != 0
543

544
static INLINE
545
void
546
ss_blockswap(int *a, int *b, int n) {
547
  int t;
548
  for(; 0 < n; --n, ++a, ++b) {
549
    t = *a, *a = *b, *b = t;
550
  }
551
}
552

553
static INLINE
554
void
555
ss_rotate(int *first, int *middle, int *last) {
556
  int *a, *b, t;
557
  int l, r;
558
  l = middle - first, r = last - middle;
559
  for(; (0 < l) && (0 < r);) {
560
    if(l == r) { ss_blockswap(first, middle, l); break; }
561
    if(l < r) {
562
      a = last - 1, b = middle - 1;
563
      t = *a;
564
      do {
565
        *a-- = *b, *b-- = *a;
566
        if(b < first) {
567
          *a = t;
568
          last = a;
569
          if((r -= l + 1) <= l) { break; }
570
          a -= 1, b = middle - 1;
571
          t = *a;
572
        }
573
      } while(1);
574
    } else {
575
      a = first, b = middle;
576
      t = *a;
577
      do {
578
        *a++ = *b, *b++ = *a;
579
        if(last <= b) {
580
          *a = t;
581
          first = a + 1;
582
          if((l -= r + 1) <= r) { break; }
583
          a += 1, b = middle;
584
          t = *a;
585
        }
586
      } while(1);
587
    }
588
  }
589
}
590

591

592
/*---------------------------------------------------------------------------*/
593

594
static
595
void
596
ss_inplacemerge(const unsigned char *T, const int *PA,
597
                int *first, int *middle, int *last,
598
                int depth) {
599
  const int *p;
600
  int *a, *b;
601
  int len, half;
602
  int q, r;
603
  int x;
604

605
  for(;;) {
606
    if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); }
607
    else                { x = 0; p = PA +  *(last - 1); }
608
    for(a = first, len = middle - first, half = len >> 1, r = -1;
609
        0 < len;
610
        len = half, half >>= 1) {
611
      b = a + half;
612
      q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
613
      if(q < 0) {
614
        a = b + 1;
615
        half -= (len & 1) ^ 1;
616
      } else {
617
        r = q;
618
      }
619
    }
620
    if(a < middle) {
621
      if(r == 0) { *a = ~*a; }
622
      ss_rotate(a, middle, last);
623
      last -= middle - a;
624
      middle = a;
625
      if(first == middle) { break; }
626
    }
627
    --last;
628
    if(x != 0) { while(*--last < 0) { } }
629
    if(middle == last) { break; }
630
  }
631
}
632

633

634
/*---------------------------------------------------------------------------*/
635

636
/* Merge-forward with internal buffer. */
637
static
638
void
639
ss_mergeforward(const unsigned char *T, const int *PA,
640
                int *first, int *middle, int *last,
641
                int *buf, int depth) {
642
  int *a, *b, *c, *bufend;
643
  int t;
644
  int r;
645

646
  bufend = buf + (middle - first) - 1;
647
  ss_blockswap(buf, first, middle - first);
648

649
  for(t = *(a = first), b = buf, c = middle;;) {
650
    r = ss_compare(T, PA + *b, PA + *c, depth);
651
    if(r < 0) {
652
      do {
653
        *a++ = *b;
654
        if(bufend <= b) { *bufend = t; return; }
655
        *b++ = *a;
656
      } while(*b < 0);
657
    } else if(r > 0) {
658
      do {
659
        *a++ = *c, *c++ = *a;
660
        if(last <= c) {
661
          while(b < bufend) { *a++ = *b, *b++ = *a; }
662
          *a = *b, *b = t;
663
          return;
664
        }
665
      } while(*c < 0);
666
    } else {
667
      *c = ~*c;
668
      do {
669
        *a++ = *b;
670
        if(bufend <= b) { *bufend = t; return; }
671
        *b++ = *a;
672
      } while(*b < 0);
673

674
      do {
675
        *a++ = *c, *c++ = *a;
676
        if(last <= c) {
677
          while(b < bufend) { *a++ = *b, *b++ = *a; }
678
          *a = *b, *b = t;
679
          return;
680
        }
681
      } while(*c < 0);
682
    }
683
  }
684
}
685

686
/* Merge-backward with internal buffer. */
687
static
688
void
689
ss_mergebackward(const unsigned char *T, const int *PA,
690
                 int *first, int *middle, int *last,
691
                 int *buf, int depth) {
692
  const int *p1, *p2;
693
  int *a, *b, *c, *bufend;
694
  int t;
695
  int r;
696
  int x;
697

698
  bufend = buf + (last - middle) - 1;
699
  ss_blockswap(buf, middle, last - middle);
700

701
  x = 0;
702
  if(*bufend < 0)       { p1 = PA + ~*bufend; x |= 1; }
703
  else                  { p1 = PA +  *bufend; }
704
  if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; }
705
  else                  { p2 = PA +  *(middle - 1); }
706
  for(t = *(a = last - 1), b = bufend, c = middle - 1;;) {
707
    r = ss_compare(T, p1, p2, depth);
708
    if(0 < r) {
709
      if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
710
      *a-- = *b;
711
      if(b <= buf) { *buf = t; break; }
712
      *b-- = *a;
713
      if(*b < 0) { p1 = PA + ~*b; x |= 1; }
714
      else       { p1 = PA +  *b; }
715
    } else if(r < 0) {
716
      if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
717
      *a-- = *c, *c-- = *a;
718
      if(c < first) {
719
        while(buf < b) { *a-- = *b, *b-- = *a; }
720
        *a = *b, *b = t;
721
        break;
722
      }
723
      if(*c < 0) { p2 = PA + ~*c; x |= 2; }
724
      else       { p2 = PA +  *c; }
725
    } else {
726
      if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
727
      *a-- = ~*b;
728
      if(b <= buf) { *buf = t; break; }
729
      *b-- = *a;
730
      if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
731
      *a-- = *c, *c-- = *a;
732
      if(c < first) {
733
        while(buf < b) { *a-- = *b, *b-- = *a; }
734
        *a = *b, *b = t;
735
        break;
736
      }
737
      if(*b < 0) { p1 = PA + ~*b; x |= 1; }
738
      else       { p1 = PA +  *b; }
739
      if(*c < 0) { p2 = PA + ~*c; x |= 2; }
740
      else       { p2 = PA +  *c; }
741
    }
742
  }
743
}
744

745
/* D&C based merge. */
746
static
747
void
748
ss_swapmerge(const unsigned char *T, const int *PA,
749
             int *first, int *middle, int *last,
750
             int *buf, int bufsize, int depth) {
751
#define STACK_SIZE SS_SMERGE_STACKSIZE
752
#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
753
#define MERGE_CHECK(a, b, c)\
754
  do {\
755
    if(((c) & 1) ||\
756
       (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\
757
      *(a) = ~*(a);\
758
    }\
759
    if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\
760
      *(b) = ~*(b);\
761
    }\
762
  } while(0)
763
  struct { int *a, *b, *c; int d; } stack[STACK_SIZE];
764
  int *l, *r, *lm, *rm;
765
  int m, len, half;
766
  int ssize;
767
  int check, next;
768

769
  for(check = 0, ssize = 0;;) {
770
    if((last - middle) <= bufsize) {
771
      if((first < middle) && (middle < last)) {
772
        ss_mergebackward(T, PA, first, middle, last, buf, depth);
773
      }
774
      MERGE_CHECK(first, last, check);
775
      STACK_POP(first, middle, last, check);
776
      continue;
777
    }
778

779
    if((middle - first) <= bufsize) {
780
      if(first < middle) {
781
        ss_mergeforward(T, PA, first, middle, last, buf, depth);
782
      }
783
      MERGE_CHECK(first, last, check);
784
      STACK_POP(first, middle, last, check);
785
      continue;
786
    }
787

788
    for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1;
789
        0 < len;
790
        len = half, half >>= 1) {
791
      if(ss_compare(T, PA + GETIDX(*(middle + m + half)),
792
                       PA + GETIDX(*(middle - m - half - 1)), depth) < 0) {
793
        m += half + 1;
794
        half -= (len & 1) ^ 1;
795
      }
796
    }
797

798
    if(0 < m) {
799
      lm = middle - m, rm = middle + m;
800
      ss_blockswap(lm, middle, m);
801
      l = r = middle, next = 0;
802
      if(rm < last) {
803
        if(*rm < 0) {
804
          *rm = ~*rm;
805
          if(first < lm) { for(; *--l < 0;) { } next |= 4; }
806
          next |= 1;
807
        } else if(first < lm) {
808
          for(; *r < 0; ++r) { }
809
          next |= 2;
810
        }
811
      }
812

813
      if((l - first) <= (last - r)) {
814
        STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
815
        middle = lm, last = l, check = (check & 3) | (next & 4);
816
      } else {
817
        if((next & 2) && (r == middle)) { next ^= 6; }
818
        STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
819
        first = r, middle = rm, check = (next & 3) | (check & 4);
820
      }
821
    } else {
822
      if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) {
823
        *middle = ~*middle;
824
      }
825
      MERGE_CHECK(first, last, check);
826
      STACK_POP(first, middle, last, check);
827
    }
828
  }
829
#undef STACK_SIZE
830
}
831

832
#endif /* SS_BLOCKSIZE != 0 */
833

834

835
/*---------------------------------------------------------------------------*/
836

837
/* Substring sort */
838
static
839
void
840
sssort(const unsigned char *T, const int *PA,
841
       int *first, int *last,
842
       int *buf, int bufsize,
843
       int depth, int n, int lastsuffix) {
844
  int *a;
845
#if SS_BLOCKSIZE != 0
846
  int *b, *middle, *curbuf;
847
  int j, k, curbufsize, limit;
848
#endif
849
  int i;
850

851
  if(lastsuffix != 0) { ++first; }
852

853
#if SS_BLOCKSIZE == 0
854
  ss_mintrosort(T, PA, first, last, depth);
855
#else
856
  if((bufsize < SS_BLOCKSIZE) &&
857
      (bufsize < (last - first)) &&
858
      (bufsize < (limit = ss_isqrt(last - first)))) {
859
    if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; }
860
    buf = middle = last - limit, bufsize = limit;
861
  } else {
862
    middle = last, limit = 0;
863
  }
864
  for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) {
865
#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
866
    ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth);
867
#elif 1 < SS_BLOCKSIZE
868
    ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth);
869
#endif
870
    curbufsize = last - (a + SS_BLOCKSIZE);
871
    curbuf = a + SS_BLOCKSIZE;
872
    if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
873
    for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) {
874
      ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
875
    }
876
  }
877
#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
878
  ss_mintrosort(T, PA, a, middle, depth);
879
#elif 1 < SS_BLOCKSIZE
880
  ss_insertionsort(T, PA, a, middle, depth);
881
#endif
882
  for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) {
883
    if(i & 1) {
884
      ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
885
      a -= k;
886
    }
887
  }
888
  if(limit != 0) {
889
#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
890
    ss_mintrosort(T, PA, middle, last, depth);
891
#elif 1 < SS_BLOCKSIZE
892
    ss_insertionsort(T, PA, middle, last, depth);
893
#endif
894
    ss_inplacemerge(T, PA, first, middle, last, depth);
895
  }
896
#endif
897

898
  if(lastsuffix != 0) {
899
    /* Insert last type B* suffix. */
900
    int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
901
    for(a = first, i = *(first - 1);
902
        (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth)));
903
        ++a) {
904
      *(a - 1) = *a;
905
    }
906
    *(a - 1) = i;
907
  }
908
}
909

910

911
/*---------------------------------------------------------------------------*/
912

913
static INLINE
914
int
915
tr_ilg(int n) {
916
  return (n & 0xffff0000) ?
917
          ((n & 0xff000000) ?
918
            24 + lg_table[(n >> 24) & 0xff] :
919
            16 + lg_table[(n >> 16) & 0xff]) :
920
          ((n & 0x0000ff00) ?
921
             8 + lg_table[(n >>  8) & 0xff] :
922
             0 + lg_table[(n >>  0) & 0xff]);
923
}
924

925

926
/*---------------------------------------------------------------------------*/
927

928
/* Simple insertionsort for small size groups. */
929
static
930
void
931
tr_insertionsort(const int *ISAd, int *first, int *last) {
932
  int *a, *b;
933
  int t, r;
934

935
  for(a = first + 1; a < last; ++a) {
936
    for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) {
937
      do { *(b + 1) = *b; } while((first <= --b) && (*b < 0));
938
      if(b < first) { break; }
939
    }
940
    if(r == 0) { *b = ~*b; }
941
    *(b + 1) = t;
942
  }
943
}
944

945

946
/*---------------------------------------------------------------------------*/
947

948
static INLINE
949
void
950
tr_fixdown(const int *ISAd, int *SA, int i, int size) {
951
  int j, k;
952
  int v;
953
  int c, d, e;
954

955
  for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
956
    d = ISAd[SA[k = j++]];
957
    if(d < (e = ISAd[SA[j]])) { k = j; d = e; }
958
    if(d <= c) { break; }
959
  }
960
  SA[i] = v;
961
}
962

963
/* Simple top-down heapsort. */
964
static
965
void
966
tr_heapsort(const int *ISAd, int *SA, int size) {
967
  int i, m;
968
  int t;
969

970
  m = size;
971
  if((size % 2) == 0) {
972
    m--;
973
    if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); }
974
  }
975

976
  for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); }
977
  if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); }
978
  for(i = m - 1; 0 < i; --i) {
979
    t = SA[0], SA[0] = SA[i];
980
    tr_fixdown(ISAd, SA, 0, i);
981
    SA[i] = t;
982
  }
983
}
984

985

986
/*---------------------------------------------------------------------------*/
987

988
/* Returns the median of three elements. */
989
static INLINE
990
int *
991
tr_median3(const int *ISAd, int *v1, int *v2, int *v3) {
992
  int *t;
993
  if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); }
994
  if(ISAd[*v2] > ISAd[*v3]) {
995
    if(ISAd[*v1] > ISAd[*v3]) { return v1; }
996
    else { return v3; }
997
  }
998
  return v2;
999
}
1000

1001
/* Returns the median of five elements. */
1002
static INLINE
1003
int *
1004
tr_median5(const int *ISAd,
1005
           int *v1, int *v2, int *v3, int *v4, int *v5) {
1006
  int *t;
1007
  if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); }
1008
  if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); }
1009
  if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4); SWAP(v3, v5); }
1010
  if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3); }
1011
  if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4); SWAP(v3, v5); }
1012
  if(ISAd[*v3] > ISAd[*v4]) { return v4; }
1013
  return v3;
1014
}
1015

1016
/* Returns the pivot element. */
1017
static INLINE
1018
int *
1019
tr_pivot(const int *ISAd, int *first, int *last) {
1020
  int *middle;
1021
  int t;
1022

1023
  t = last - first;
1024
  middle = first + t / 2;
1025

1026
  if(t <= 512) {
1027
    if(t <= 32) {
1028
      return tr_median3(ISAd, first, middle, last - 1);
1029
    } else {
1030
      t >>= 2;
1031
      return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
1032
    }
1033
  }
1034
  t >>= 3;
1035
  first  = tr_median3(ISAd, first, first + t, first + (t << 1));
1036
  middle = tr_median3(ISAd, middle - t, middle, middle + t);
1037
  last   = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
1038
  return tr_median3(ISAd, first, middle, last);
1039
}
1040

1041

1042
/*---------------------------------------------------------------------------*/
1043

1044
typedef struct _trbudget_t trbudget_t;
1045
struct _trbudget_t {
1046
  int chance;
1047
  int remain;
1048
  int incval;
1049
  int count;
1050
};
1051

1052
static INLINE
1053
void
1054
trbudget_init(trbudget_t *budget, int chance, int incval) {
1055
  budget->chance = chance;
1056
  budget->remain = budget->incval = incval;
1057
}
1058

1059
static INLINE
1060
int
1061
trbudget_check(trbudget_t *budget, int size) {
1062
  if(size <= budget->remain) { budget->remain -= size; return 1; }
1063
  if(budget->chance == 0) { budget->count += size; return 0; }
1064
  budget->remain += budget->incval - size;
1065
  budget->chance -= 1;
1066
  return 1;
1067
}
1068

1069

1070
/*---------------------------------------------------------------------------*/
1071

1072
static INLINE
1073
void
1074
tr_partition(const int *ISAd,
1075
             int *first, int *middle, int *last,
1076
             int **pa, int **pb, int v) {
1077
  int *a, *b, *c, *d, *e, *f;
1078
  int t, s;
1079
  int x = 0;
1080

1081
  for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { }
1082
  if(((a = b) < last) && (x < v)) {
1083
    for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
1084
      if(x == v) { SWAP(*b, *a); ++a; }
1085
    }
1086
  }
1087
  for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { }
1088
  if((b < (d = c)) && (x > v)) {
1089
    for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1090
      if(x == v) { SWAP(*c, *d); --d; }
1091
    }
1092
  }
1093
  for(; b < c;) {
1094
    SWAP(*b, *c);
1095
    for(; (++b < c) && ((x = ISAd[*b]) <= v);) {
1096
      if(x == v) { SWAP(*b, *a); ++a; }
1097
    }
1098
    for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
1099
      if(x == v) { SWAP(*c, *d); --d; }
1100
    }
1101
  }
1102

1103
  if(a <= d) {
1104
    c = b - 1;
1105
    if((s = a - first) > (t = b - a)) { s = t; }
1106
    for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1107
    if((s = d - c) > (t = last - d - 1)) { s = t; }
1108
    for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
1109
    first += (b - a), last -= (d - c);
1110
  }
1111
  *pa = first, *pb = last;
1112
}
1113

1114
static
1115
void
1116
tr_copy(int *ISA, const int *SA,
1117
        int *first, int *a, int *b, int *last,
1118
        int depth) {
1119
  /* sort suffixes of middle partition
1120
     by using sorted order of suffixes of left and right partition. */
1121
  int *c, *d, *e;
1122
  int s, v;
1123

1124
  v = b - SA - 1;
1125
  for(c = first, d = a - 1; c <= d; ++c) {
1126
    #ifdef __clang_analyzer__
1127
    assert(c);
1128
    #endif
1129
    if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1130
      *++d = s;
1131
      ISA[s] = d - SA;
1132
    }
1133
  }
1134
  for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1135
    if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1136
      *--d = s;
1137
      ISA[s] = d - SA;
1138
    }
1139
  }
1140
}
1141

1142
static
1143
void
1144
tr_partialcopy(int *ISA, const int *SA,
1145
               int *first, int *a, int *b, int *last,
1146
               int depth) {
1147
  int *c, *d, *e;
1148
  int s, v;
1149
  int rank, lastrank, newrank = -1;
1150

1151
  v = b - SA - 1;
1152
  lastrank = -1;
1153
  for(c = first, d = a - 1; c <= d; ++c) {
1154
    if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1155
      *++d = s;
1156
      rank = ISA[s + depth];
1157
      if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1158
      ISA[s] = newrank;
1159
    }
1160
  }
1161

1162
  lastrank = -1;
1163
  for(e = d; first <= e; --e) {
1164
    rank = ISA[*e];
1165
    if(lastrank != rank) { lastrank = rank; newrank = e - SA; }
1166
    if(newrank != rank) { ISA[*e] = newrank; }
1167
  }
1168

1169
  lastrank = -1;
1170
  for(c = last - 1, e = d + 1, d = b; e < d; --c) {
1171
    if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
1172
      *--d = s;
1173
      rank = ISA[s + depth];
1174
      if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
1175
      ISA[s] = newrank;
1176
    }
1177
  }
1178
}
1179

1180
static
1181
void
1182
tr_introsort(int *ISA, const int *ISAd,
1183
             int *SA, int *first, int *last,
1184
             trbudget_t *budget) {
1185
#define STACK_SIZE TR_STACKSIZE
1186
  struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE];
1187
  int *a, *b, *c;
1188
  int t;
1189
  int v, x = 0;
1190
  int incr = ISAd - ISA;
1191
  int limit, next;
1192
  int ssize, trlink = -1;
1193

1194
  #ifdef __clang_analyzer__
1195
  memset(stack, 0, sizeof(stack));
1196
  #endif
1197

1198
  for(ssize = 0, limit = tr_ilg(last - first);;) {
1199

1200
    if(limit < 0) {
1201
      if(limit == -1) {
1202
        /* tandem repeat partition */
1203
        tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1);
1204

1205
        /* update ranks */
1206
        if(a < last) {
1207
          for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1208
        }
1209
        if(b < last) {
1210
          for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; }
1211
        }
1212

1213
        /* push */
1214
        if(1 < (b - a)) {
1215
          STACK_PUSH5(NULL, a, b, 0, 0);
1216
          STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
1217
          trlink = ssize - 2;
1218
        }
1219
        if((a - first) <= (last - b)) {
1220
          if(1 < (a - first)) {
1221
            STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink);
1222
            last = a, limit = tr_ilg(a - first);
1223
          } else if(1 < (last - b)) {
1224
            first = b, limit = tr_ilg(last - b);
1225
          } else {
1226
            STACK_POP5(ISAd, first, last, limit, trlink);
1227
          }
1228
        } else {
1229
          if(1 < (last - b)) {
1230
            STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink);
1231
            first = b, limit = tr_ilg(last - b);
1232
          } else if(1 < (a - first)) {
1233
            last = a, limit = tr_ilg(a - first);
1234
          } else {
1235
            STACK_POP5(ISAd, first, last, limit, trlink);
1236
          }
1237
        }
1238
      } else if(limit == -2) {
1239
        /* tandem repeat copy */
1240
        a = stack[--ssize].b, b = stack[ssize].c;
1241
        if(stack[ssize].d == 0) {
1242
          tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
1243
        } else {
1244
          if(0 <= trlink) { stack[trlink].d = -1; }
1245
          tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
1246
        }
1247
        STACK_POP5(ISAd, first, last, limit, trlink);
1248
      } else {
1249
        /* sorted partition */
1250
        if(0 <= *first) {
1251
          a = first;
1252
          do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a));
1253
          first = a;
1254
        }
1255
        if(first < last) {
1256
          a = first; do { *a = ~*a; } while(*++a < 0);
1257
          next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1;
1258
          if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } }
1259

1260
          /* push */
1261
          if(trbudget_check(budget, a - first)) {
1262
            if((a - first) <= (last - a)) {
1263
              STACK_PUSH5(ISAd, a, last, -3, trlink);
1264
              ISAd += incr, last = a, limit = next;
1265
            } else {
1266
              if(1 < (last - a)) {
1267
                STACK_PUSH5(ISAd + incr, first, a, next, trlink);
1268
                first = a, limit = -3;
1269
              } else {
1270
                ISAd += incr, last = a, limit = next;
1271
              }
1272
            }
1273
          } else {
1274
            if(0 <= trlink) { stack[trlink].d = -1; }
1275
            if(1 < (last - a)) {
1276
              first = a, limit = -3;
1277
            } else {
1278
              STACK_POP5(ISAd, first, last, limit, trlink);
1279
            }
1280
          }
1281
        } else {
1282
          STACK_POP5(ISAd, first, last, limit, trlink);
1283
        }
1284
      }
1285
      continue;
1286
    }
1287

1288
    if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
1289
      tr_insertionsort(ISAd, first, last);
1290
      limit = -3;
1291
      continue;
1292
    }
1293

1294
    if(limit-- == 0) {
1295
      tr_heapsort(ISAd, first, last - first);
1296
      for(a = last - 1; first < a; a = b) {
1297
        for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
1298
      }
1299
      limit = -3;
1300
      continue;
1301
    }
1302

1303
    /* choose pivot */
1304
    a = tr_pivot(ISAd, first, last);
1305
    SWAP(*first, *a);
1306
    v = ISAd[*first];
1307

1308
    /* partition */
1309
    tr_partition(ISAd, first, first + 1, last, &a, &b, v);
1310
    if((last - first) != (b - a)) {
1311
      next = (ISA[*a] != v) ? tr_ilg(b - a) : -1;
1312

1313
      /* update ranks */
1314
      for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
1315
      if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } }
1316

1317
      /* push */
1318
      if((1 < (b - a)) && (trbudget_check(budget, b - a))) {
1319
        if((a - first) <= (last - b)) {
1320
          if((last - b) <= (b - a)) {
1321
            if(1 < (a - first)) {
1322
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1323
              STACK_PUSH5(ISAd, b, last, limit, trlink);
1324
              last = a;
1325
            } else if(1 < (last - b)) {
1326
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1327
              first = b;
1328
            } else {
1329
              ISAd += incr, first = a, last = b, limit = next;
1330
            }
1331
          } else if((a - first) <= (b - a)) {
1332
            if(1 < (a - first)) {
1333
              STACK_PUSH5(ISAd, b, last, limit, trlink);
1334
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1335
              last = a;
1336
            } else {
1337
              STACK_PUSH5(ISAd, b, last, limit, trlink);
1338
              ISAd += incr, first = a, last = b, limit = next;
1339
            }
1340
          } else {
1341
            STACK_PUSH5(ISAd, b, last, limit, trlink);
1342
            STACK_PUSH5(ISAd, first, a, limit, trlink);
1343
            ISAd += incr, first = a, last = b, limit = next;
1344
          }
1345
        } else {
1346
          if((a - first) <= (b - a)) {
1347
            if(1 < (last - b)) {
1348
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1349
              STACK_PUSH5(ISAd, first, a, limit, trlink);
1350
              first = b;
1351
            } else if(1 < (a - first)) {
1352
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1353
              last = a;
1354
            } else {
1355
              ISAd += incr, first = a, last = b, limit = next;
1356
            }
1357
          } else if((last - b) <= (b - a)) {
1358
            if(1 < (last - b)) {
1359
              STACK_PUSH5(ISAd, first, a, limit, trlink);
1360
              STACK_PUSH5(ISAd + incr, a, b, next, trlink);
1361
              first = b;
1362
            } else {
1363
              STACK_PUSH5(ISAd, first, a, limit, trlink);
1364
              ISAd += incr, first = a, last = b, limit = next;
1365
            }
1366
          } else {
1367
            STACK_PUSH5(ISAd, first, a, limit, trlink);
1368
            STACK_PUSH5(ISAd, b, last, limit, trlink);
1369
            ISAd += incr, first = a, last = b, limit = next;
1370
          }
1371
        }
1372
      } else {
1373
        if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
1374
        if((a - first) <= (last - b)) {
1375
          if(1 < (a - first)) {
1376
            STACK_PUSH5(ISAd, b, last, limit, trlink);
1377
            last = a;
1378
          } else if(1 < (last - b)) {
1379
            first = b;
1380
          } else {
1381
            STACK_POP5(ISAd, first, last, limit, trlink);
1382
          }
1383
        } else {
1384
          if(1 < (last - b)) {
1385
            STACK_PUSH5(ISAd, first, a, limit, trlink);
1386
            first = b;
1387
          } else if(1 < (a - first)) {
1388
            last = a;
1389
          } else {
1390
            STACK_POP5(ISAd, first, last, limit, trlink);
1391
          }
1392
        }
1393
      }
1394
    } else {
1395
      if(trbudget_check(budget, last - first)) {
1396
        limit = tr_ilg(last - first), ISAd += incr;
1397
      } else {
1398
        if(0 <= trlink) { stack[trlink].d = -1; }
1399
        STACK_POP5(ISAd, first, last, limit, trlink);
1400
      }
1401
    }
1402
  }
1403
#undef STACK_SIZE
1404
}
1405

1406

1407

1408
/*---------------------------------------------------------------------------*/
1409

1410
/* Tandem repeat sort */
1411
static
1412
void
1413
trsort(int *ISA, int *SA, int n, int depth) {
1414
  int *ISAd;
1415
  int *first, *last;
1416
  trbudget_t budget;
1417
  int t, skip, unsorted;
1418

1419
  trbudget_init(&budget, tr_ilg(n) * 2 / 3, n);
1420
/*  trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */
1421
  for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) {
1422
    first = SA;
1423
    skip = 0;
1424
    unsorted = 0;
1425
    do {
1426
      if((t = *first) < 0) { first -= t; skip += t; }
1427
      else {
1428
        if(skip != 0) { *(first + skip) = skip; skip = 0; }
1429
        last = SA + ISA[t] + 1;
1430
        if(1 < (last - first)) {
1431
          budget.count = 0;
1432
          tr_introsort(ISA, ISAd, SA, first, last, &budget);
1433
          if(budget.count != 0) { unsorted += budget.count; }
1434
          else { skip = first - last; }
1435
        } else if((last - first) == 1) {
1436
          skip = -1;
1437
        }
1438
        first = last;
1439
      }
1440
    } while(first < (SA + n));
1441
    if(skip != 0) { *(first + skip) = skip; }
1442
    if(unsorted == 0) { break; }
1443
  }
1444
}
1445

1446

1447
/*---------------------------------------------------------------------------*/
1448

1449
/* Sorts suffixes of type B*. */
1450
static
1451
int
1452
sort_typeBstar(const unsigned char *T, int *SA,
1453
               int *bucket_A, int *bucket_B,
1454
               int n, int openMP) {
1455
  int *PAb, *ISAb, *buf;
1456
#ifdef LIBBSC_OPENMP
1457
  int *curbuf;
1458
  int l;
1459
#endif
1460
  int i, j, k, t, m, bufsize;
1461
  int c0, c1;
1462
#ifdef LIBBSC_OPENMP
1463
  int d0, d1;
1464
#endif
1465
  (void)openMP;
1466

1467
  /* Initialize bucket arrays. */
1468
  for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; }
1469
  for(i = 0; i < BUCKET_B_SIZE; ++i) { bucket_B[i] = 0; }
1470

1471
  /* Count the number of occurrences of the first one or two characters of each
1472
     type A, B and B* suffix. Moreover, store the beginning position of all
1473
     type B* suffixes into the array SA. */
1474
  for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) {
1475
    /* type A suffix. */
1476
    do { ++BUCKET_A(c1 = c0); } while((0 <= --i) && ((c0 = T[i]) >= c1));
1477
    if(0 <= i) {
1478
      /* type B* suffix. */
1479
      ++BUCKET_BSTAR(c0, c1);
1480
      SA[--m] = i;
1481
      /* type B suffix. */
1482
      for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) {
1483
        ++BUCKET_B(c0, c1);
1484
      }
1485
    }
1486
  }
1487
  m = n - m;
1488
/*
1489
note:
1490
  A type B* suffix is lexicographically smaller than a type B suffix that
1491
  begins with the same first two characters.
1492
*/
1493

1494
  /* Calculate the index of start/end point of each bucket. */
1495
  for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0) {
1496
    t = i + BUCKET_A(c0);
1497
    BUCKET_A(c0) = i + j; /* start point */
1498
    i = t + BUCKET_B(c0, c0);
1499
    for(c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1) {
1500
      j += BUCKET_BSTAR(c0, c1);
1501
      BUCKET_BSTAR(c0, c1) = j; /* end point */
1502
      i += BUCKET_B(c0, c1);
1503
    }
1504
  }
1505

1506
  if(0 < m) {
1507
    /* Sort the type B* suffixes by their first two characters. */
1508
    PAb = SA + n - m; ISAb = SA + m;
1509
    for(i = m - 2; 0 <= i; --i) {
1510
      t = PAb[i], c0 = T[t], c1 = T[t + 1];
1511
      SA[--BUCKET_BSTAR(c0, c1)] = i;
1512
    }
1513
    t = PAb[m - 1], c0 = T[t], c1 = T[t + 1];
1514
    SA[--BUCKET_BSTAR(c0, c1)] = m - 1;
1515

1516
    /* Sort the type B* substrings using sssort. */
1517
#ifdef LIBBSC_OPENMP
1518
    if (openMP)
1519
    {
1520
        buf = SA + m;
1521
        c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m;
1522
#pragma omp parallel default(shared) private(bufsize, curbuf, k, l, d0, d1)
1523
        {
1524
          bufsize = (n - (2 * m)) / omp_get_num_threads();
1525
          curbuf = buf + omp_get_thread_num() * bufsize;
1526
          k = 0;
1527
          for(;;) {
1528
            #pragma omp critical(sssort_lock)
1529
            {
1530
              if(0 < (l = j)) {
1531
                d0 = c0, d1 = c1;
1532
                do {
1533
                  k = BUCKET_BSTAR(d0, d1);
1534
                  if(--d1 <= d0) {
1535
                    d1 = ALPHABET_SIZE - 1;
1536
                    if(--d0 < 0) { break; }
1537
                  }
1538
                } while(((l - k) <= 1) && (0 < (l = k)));
1539
                c0 = d0, c1 = d1, j = k;
1540
              }
1541
            }
1542
            if(l == 0) { break; }
1543
            sssort(T, PAb, SA + k, SA + l,
1544
                   curbuf, bufsize, 2, n, *(SA + k) == (m - 1));
1545
          }
1546
        }
1547
    }
1548
    else
1549
    {
1550
        buf = SA + m, bufsize = n - (2 * m);
1551
        for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
1552
          for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
1553
            i = BUCKET_BSTAR(c0, c1);
1554
            if(1 < (j - i)) {
1555
              sssort(T, PAb, SA + i, SA + j,
1556
                     buf, bufsize, 2, n, *(SA + i) == (m - 1));
1557
            }
1558
          }
1559
        }
1560
    }
1561
#else
1562
    buf = SA + m, bufsize = n - (2 * m);
1563
    for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
1564
      for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
1565
        i = BUCKET_BSTAR(c0, c1);
1566
        if(1 < (j - i)) {
1567
          sssort(T, PAb, SA + i, SA + j,
1568
                 buf, bufsize, 2, n, *(SA + i) == (m - 1));
1569
        }
1570
      }
1571
    }
1572
#endif
1573

1574
    /* Compute ranks of type B* substrings. */
1575
    for(i = m - 1; 0 <= i; --i) {
1576
      if(0 <= SA[i]) {
1577
        j = i;
1578
        do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i]));
1579
        SA[i + 1] = i - j;
1580
        if(i <= 0) { break; }
1581
      }
1582
      j = i;
1583
      do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0);
1584
      ISAb[SA[i]] = j;
1585
    }
1586

1587
    /* Construct the inverse suffix array of type B* suffixes using trsort. */
1588
    trsort(ISAb, SA, m, 1);
1589

1590
    /* Set the sorted order of type B* suffixes. */
1591
    for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) {
1592
      for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { }
1593
      if(0 <= i) {
1594
        t = i;
1595
        for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { }
1596
        SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
1597
      }
1598
    }
1599

1600
    /* Calculate the index of start/end point of each bucket. */
1601
    BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */
1602
    for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) {
1603
      i = BUCKET_A(c0 + 1) - 1;
1604
      for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) {
1605
        t = i - BUCKET_B(c0, c1);
1606
        BUCKET_B(c0, c1) = i; /* end point */
1607

1608
        /* Move all type B* suffixes to the correct position. */
1609
        for(i = t, j = BUCKET_BSTAR(c0, c1);
1610
            j <= k;
1611
            --i, --k) { SA[i] = SA[k]; }
1612
      }
1613
      BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */
1614
      BUCKET_B(c0, c0) = i; /* end point */
1615
    }
1616
  }
1617

1618
  return m;
1619
}
1620

1621
/* Constructs the suffix array by using the sorted order of type B* suffixes. */
1622
static
1623
void
1624
construct_SA(const unsigned char *T, int *SA,
1625
             int *bucket_A, int *bucket_B,
1626
             int n, int m) {
1627
  int *i, *j, *k;
1628
  int s;
1629
  int c0, c1, c2;
1630

1631
  if(0 < m) {
1632
    /* Construct the sorted order of type B suffixes by using
1633
       the sorted order of type B* suffixes. */
1634
    for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1635
      /* Scan the suffix array from right to left. */
1636
      for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1637
          j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1638
          i <= j;
1639
          --j) {
1640
        if(0 < (s = *j)) {
1641
          assert(T[s] == c1);
1642
          assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1643
          assert(T[s - 1] <= T[s]);
1644
          *j = ~s;
1645
          c0 = T[--s];
1646
          if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1647
          if(c0 != c2) {
1648
            if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1649
            k = SA + BUCKET_B(c2 = c0, c1);
1650
          }
1651
          assert(k < j); assert(k != NULL);
1652
          *k-- = s;
1653
        } else {
1654
          assert(((s == 0) && (T[s] == c1)) || (s < 0));
1655
          *j = ~s;
1656
        }
1657
      }
1658
    }
1659
  }
1660

1661
  /* Construct the suffix array by using
1662
     the sorted order of type B suffixes. */
1663
  k = SA + BUCKET_A(c2 = T[n - 1]);
1664
  *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
1665
  /* Scan the suffix array from left to right. */
1666
  for(i = SA, j = SA + n; i < j; ++i) {
1667
    if(0 < (s = *i)) {
1668
      assert(T[s - 1] >= T[s]);
1669
      c0 = T[--s];
1670
      if((s == 0) || (T[s - 1] < c0)) { s = ~s; }
1671
      if(c0 != c2) {
1672
        BUCKET_A(c2) = k - SA;
1673
        k = SA + BUCKET_A(c2 = c0);
1674
      }
1675
      assert(i < k);
1676
      *k++ = s;
1677
    } else {
1678
      assert(s < 0);
1679
      *i = ~s;
1680
    }
1681
  }
1682
}
1683

1684
/* Constructs the burrows-wheeler transformed string directly
1685
   by using the sorted order of type B* suffixes. */
1686
static
1687
int
1688
construct_BWT(const unsigned char *T, int *SA,
1689
              int *bucket_A, int *bucket_B,
1690
              int n, int m) {
1691
  int *i, *j, *k, *orig;
1692
  int s;
1693
  int c0, c1, c2;
1694

1695
  if(0 < m) {
1696
    /* Construct the sorted order of type B suffixes by using
1697
       the sorted order of type B* suffixes. */
1698
    for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1699
      /* Scan the suffix array from right to left. */
1700
      for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1701
          j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1702
          i <= j;
1703
          --j) {
1704
        if(0 < (s = *j)) {
1705
          assert(T[s] == c1);
1706
          assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1707
          assert(T[s - 1] <= T[s]);
1708
          c0 = T[--s];
1709
          *j = ~((int)c0);
1710
          if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1711
          if(c0 != c2) {
1712
            if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1713
            k = SA + BUCKET_B(c2 = c0, c1);
1714
          }
1715
          assert(k < j); assert(k != NULL);
1716
          *k-- = s;
1717
        } else if(s != 0) {
1718
          *j = ~s;
1719
#ifndef NDEBUG
1720
        } else {
1721
          assert(T[s] == c1);
1722
#endif
1723
        }
1724
      }
1725
    }
1726
  }
1727

1728
  /* Construct the BWTed string by using
1729
     the sorted order of type B suffixes. */
1730
  k = SA + BUCKET_A(c2 = T[n - 1]);
1731
  *k++ = (T[n - 2] < c2) ? ~((int)T[n - 2]) : (n - 1);
1732
  /* Scan the suffix array from left to right. */
1733
  for(i = SA, j = SA + n, orig = SA; i < j; ++i) {
1734
    if(0 < (s = *i)) {
1735
      assert(T[s - 1] >= T[s]);
1736
      c0 = T[--s];
1737
      *i = c0;
1738
      if((0 < s) && (T[s - 1] < c0)) { s = ~((int)T[s - 1]); }
1739
      if(c0 != c2) {
1740
        BUCKET_A(c2) = k - SA;
1741
        k = SA + BUCKET_A(c2 = c0);
1742
      }
1743
      assert(i < k);
1744
      *k++ = s;
1745
    } else if(s != 0) {
1746
      *i = ~s;
1747
    } else {
1748
      orig = i;
1749
    }
1750
  }
1751

1752
  return orig - SA;
1753
}
1754

1755
/* Constructs the burrows-wheeler transformed string directly
1756
   by using the sorted order of type B* suffixes. */
1757
static
1758
int
1759
construct_BWT_indexes(const unsigned char *T, int *SA,
1760
                      int *bucket_A, int *bucket_B,
1761
                      int n, int m,
1762
                      unsigned char * num_indexes, int * indexes) {
1763
  int *i, *j, *k, *orig;
1764
  int s;
1765
  int c0, c1, c2;
1766

1767
  int mod = n / 8;
1768
  {
1769
      mod |= mod >> 1;  mod |= mod >> 2;
1770
      mod |= mod >> 4;  mod |= mod >> 8;
1771
      mod |= mod >> 16; mod >>= 1;
1772

1773
      *num_indexes = (unsigned char)((n - 1) / (mod + 1));
1774
  }
1775

1776
  if(0 < m) {
1777
    /* Construct the sorted order of type B suffixes by using
1778
       the sorted order of type B* suffixes. */
1779
    for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
1780
      /* Scan the suffix array from right to left. */
1781
      for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
1782
          j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
1783
          i <= j;
1784
          --j) {
1785
        if(0 < (s = *j)) {
1786
          assert(T[s] == c1);
1787
          assert(((s + 1) < n) && (T[s] <= T[s + 1]));
1788
          assert(T[s - 1] <= T[s]);
1789

1790
          if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = j - SA;
1791

1792
          c0 = T[--s];
1793
          *j = ~((int)c0);
1794
          if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
1795
          if(c0 != c2) {
1796
            if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
1797
            k = SA + BUCKET_B(c2 = c0, c1);
1798
          }
1799
          assert(k < j); assert(k != NULL);
1800
          *k-- = s;
1801
        } else if(s != 0) {
1802
          *j = ~s;
1803
#ifndef NDEBUG
1804
        } else {
1805
          assert(T[s] == c1);
1806
#endif
1807
        }
1808
      }
1809
    }
1810
  }
1811

1812
  /* Construct the BWTed string by using
1813
     the sorted order of type B suffixes. */
1814
  k = SA + BUCKET_A(c2 = T[n - 1]);
1815
  if (T[n - 2] < c2) {
1816
    if (((n - 1) & mod) == 0) indexes[(n - 1) / (mod + 1) - 1] = k - SA;
1817
    *k++ = ~((int)T[n - 2]);
1818
  }
1819
  else {
1820
    *k++ = n - 1;
1821
  }
1822

1823
  /* Scan the suffix array from left to right. */
1824
  for(i = SA, j = SA + n, orig = SA; i < j; ++i) {
1825
    if(0 < (s = *i)) {
1826
      assert(T[s - 1] >= T[s]);
1827

1828
      if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = i - SA;
1829

1830
      c0 = T[--s];
1831
      *i = c0;
1832
      if(c0 != c2) {
1833
        BUCKET_A(c2) = k - SA;
1834
        k = SA + BUCKET_A(c2 = c0);
1835
      }
1836
      assert(i < k);
1837
      if((0 < s) && (T[s - 1] < c0)) {
1838
          if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = k - SA;
1839
          *k++ = ~((int)T[s - 1]);
1840
      } else
1841
        *k++ = s;
1842
    } else if(s != 0) {
1843
      *i = ~s;
1844
    } else {
1845
      orig = i;
1846
    }
1847
  }
1848

1849
  return orig - SA;
1850
}
1851

1852

1853
/*---------------------------------------------------------------------------*/
1854

1855
/*- Function -*/
1856

1857
int
1858
divsufsort(const unsigned char *T, int *SA, int n, int openMP) {
1859
  int *bucket_A, *bucket_B;
1860
  int m;
1861
  int err = 0;
1862

1863
  /* Check arguments. */
1864
  if((T == NULL) || (SA == NULL) || (n < 0)) { return -1; }
1865
  else if(n == 0) { return 0; }
1866
  else if(n == 1) { SA[0] = 0; return 0; }
1867
  else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; }
1868

1869
  bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int));
1870
  bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int));
1871

1872
  /* Suffixsort. */
1873
  if((bucket_A != NULL) && (bucket_B != NULL)) {
1874
    m = sort_typeBstar(T, SA, bucket_A, bucket_B, n, openMP);
1875
    construct_SA(T, SA, bucket_A, bucket_B, n, m);
1876
  } else {
1877
    err = -2;
1878
  }
1879

1880
  free(bucket_B);
1881
  free(bucket_A);
1882

1883
  return err;
1884
}
1885

1886
int
1887
divbwt(const unsigned char *T, unsigned char *U, int *A, int n, unsigned char * num_indexes, int * indexes, int openMP) {
1888
  int *B;
1889
  int *bucket_A, *bucket_B;
1890
  int m, pidx, i;
1891

1892
  /* Check arguments. */
1893
  if((T == NULL) || (U == NULL) || (n < 0)) { return -1; }
1894
  else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; }
1895

1896
  if((B = A) == NULL) { B = (int *)malloc((size_t)(n + 1) * sizeof(int)); }
1897
  bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int));
1898
  bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int));
1899

1900
  /* Burrows-Wheeler Transform. */
1901
  if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) {
1902
    m = sort_typeBstar(T, B, bucket_A, bucket_B, n, openMP);
1903

1904
    if (num_indexes == NULL || indexes == NULL) {
1905
        pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m);
1906
    } else {
1907
        pidx = construct_BWT_indexes(T, B, bucket_A, bucket_B, n, m, num_indexes, indexes);
1908
    }
1909

1910
    /* Copy to output string. */
1911
    U[0] = T[n - 1];
1912
    for(i = 0; i < pidx; ++i) { U[i + 1] = (unsigned char)B[i]; }
1913
    for(i += 1; i < n; ++i) { U[i] = (unsigned char)B[i]; }
1914
    pidx += 1;
1915
  } else {
1916
    pidx = -2;
1917
  }
1918

1919
  free(bucket_B);
1920
  free(bucket_A);
1921
  if(A == NULL) { free(B); }
1922

1923
  return pidx;
1924
}
1925

1926
Product

Resources

Company
