source: thirdparty/blosc/internal-complibs/zstd-0.7.4/dictBuilder/divsufsort.c @ 8ebc79b

Revision 8ebc79b, 53.3 KB checked in by Hal Finkel <hfinkel@…>, 8 years ago (diff)

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