incbwt/qsufsort/qsufsort.c

   1 /* qsufsort.c
   2    Copyright 1999, N. Jesper Larsson, all rights reserved.
   3
   4    This file contains an implementation of the algorithm presented in "Faster
   5    Suffix Sorting" by N. Jesper Larsson (jesper@@cs.lth.se) and Kunihiko
   6    Sadakane (sada@@is.s.u-tokyo.ac.jp).
   7
   8    This software may be used freely for any purpose. However, when distributed,
   9    the original source must be clearly stated, and, when the source code is
  10    distributed, the copyright notice must be retained and any alterations in
  11    the code must be clearly marked. No warranty is given regarding the quality
  12    of this software.*/
  13
  14 /*
  15   Moved everything into a struct to support parallel construction.
  16   Created another version using sint instead of int to index more than 2 GB at a time.
  17   Put everything into a namespace.
  18
  19       - Jouni Siren 2009-03-20
  20 */
  21
  22 #include <climits>
  23 #include "qsufsort.h"
  24
  25
  26 namespace CSA
  27 {
  28
  29
  30 #define KEY(p)          (V[*(p)+(h)])
  31 #define SWAP(p, q)      (tmp=*(p), *(p)=*(q), *(q)=tmp)
  32 #define MED3(a, b, c)   (KEY(a)<KEY(b) ? (KEY(b)<KEY(c) ? b : KEY(a)<KEY(c) ? c : (a)) : (KEY(b)>KEY(c) ? b : KEY(a)>KEY(c) ? c : (a)))
  33
  34
  35 //--------------------------------------------------------------------------
  36
  37 template <class T>
  38 struct QSufSort
  39 {
  40
  41    T *I,                        /* group array, ultimately suffix array.*/
  42    *V,                          /* inverse array, ultimately inverse of I.*/
  43    r,                           /* number of symbols aggregated by transform.*/
  44    h;                           /* length of already-sorted prefixes.*/
  45
  46
  47 /* Subroutine for select_sort_split and sort_split. Sets group numbers for a
  48    group whose lowest position in I is pl and highest position is pm.*/
  49
  50 void update_group(T *pl, T *pm)
  51 {
  52    T g;
  53
  54    g=pm-I;                      /* group number.*/
  55    V[*pl]=g;                    /* update group number of first position.*/
  56    if (pl==pm)
  57       *pl=-1;                   /* one element, sorted group.*/
  58    else
  59       do                        /* more than one element, unsorted group.*/
  60          V[*++pl]=g;            /* update group numbers.*/
  61       while (pl<pm);
  62 }
  63
  64 /* Quadratic sorting method to use for small subarrays. To be able to update
  65    group numbers consistently, a variant of selection sorting is used.*/
  66
  67 void select_sort_split(T *p, T n) {
  68    T *pa, *pb, *pi, *pn;
  69    T f, v, tmp;
  70
  71    pa=p;                        /* pa is start of group being picked out.*/
  72    pn=p+n-1;                    /* pn is last position of subarray.*/
  73    while (pa<pn) {
  74       for (pi=pb=pa+1, f=KEY(pa); pi<=pn; ++pi)
  75          if ((v=KEY(pi))<f) {
  76             f=v;                /* f is smallest key found.*/
  77             SWAP(pi, pa);       /* place smallest element at beginning.*/
  78             pb=pa+1;            /* pb is position for elements equal to f.*/
  79          } else if (v==f) {     /* if equal to smallest key.*/
  80             SWAP(pi, pb);       /* place next to other smallest elements.*/
  81             ++pb;
  82          }
  83       update_group(pa, pb-1);   /* update group values for new group.*/
  84       pa=pb;                    /* continue sorting rest of the subarray.*/
  85    }
  86    if (pa==pn) {                /* check if last part is single element.*/
  87       V[*pa]=pa-I;
  88       *pa=-1;                   /* sorted group.*/
  89    }
  90 }
  91
  92 /* Subroutine for sort_split, algorithm by Bentley & McIlroy.*/
  93
  94 T choose_pivot(T *p, T n) {
  95    T *pl, *pm, *pn;
  96    T s;
  97
  98    pm=p+(n>>1);                 /* small arrays, middle element.*/
  99    if (n>7) {
 100       pl=p;
 101       pn=p+n-1;
 102       if (n>40) {               /* big arrays, pseudomedian of 9.*/
 103          s=n>>3;
 104          pl=MED3(pl, pl+s, pl+s+s);
 105          pm=MED3(pm-s, pm, pm+s);
 106          pn=MED3(pn-s-s, pn-s, pn);
 107       }
 108       pm=MED3(pl, pm, pn);      /* midsize arrays, median of 3.*/
 109    }
 110    return KEY(pm);
 111 }
 112
 113 /* Sorting routine called for each unsorted group. Sorts the array of integers
 114    (suffix numbers) of length n starting at p. The algorithm is a ternary-split
 115    quicksort taken from Bentley & McIlroy, "Engineering a Sort Function",
 116    Software -- Practice and Experience 23(11), 1249-1265 (November 1993). This
 117    function is based on Program 7.*/
 118
 119 void sort_split(T *p, T n)
 120 {
 121    T *pa, *pb, *pc, *pd, *pl, *pm, *pn;
 122    T f, v, s, t, tmp;
 123
 124    if (n<7) {                   /* multi-selection sort smallest arrays.*/
 125       select_sort_split(p, n);
 126       return;
 127    }
 128
 129    v=choose_pivot(p, n);
 130    pa=pb=p;
 131    pc=pd=p+n-1;
 132    while (1) {                  /* split-end partition.*/
 133       while (pb<=pc && (f=KEY(pb))<=v) {
 134          if (f==v) {
 135             SWAP(pa, pb);
 136             ++pa;
 137          }
 138          ++pb;
 139       }
 140       while (pc>=pb && (f=KEY(pc))>=v) {
 141          if (f==v) {
 142             SWAP(pc, pd);
 143             --pd;
 144          }
 145          --pc;
 146       }
 147       if (pb>pc)
 148          break;
 149       SWAP(pb, pc);
 150       ++pb;
 151       --pc;
 152    }
 153    pn=p+n;
 154    if ((s=pa-p)>(t=pb-pa))
 155       s=t;
 156    for (pl=p, pm=pb-s; s; --s, ++pl, ++pm)
 157       SWAP(pl, pm);
 158    if ((s=pd-pc)>(t=pn-pd-1))
 159       s=t;
 160    for (pl=pb, pm=pn-s; s; --s, ++pl, ++pm)
 161       SWAP(pl, pm);
 162
 163    s=pb-pa;
 164    t=pd-pc;
 165    if (s>0)
 166       sort_split(p, s);
 167    update_group(p+s, p+n-t-1);
 168    if (t>0)
 169       sort_split(p+n-t, t);
 170 }
 171
 172 /* Bucketsort for first iteration.
 173
 174    Input: x[0...n-1] holds integers in the range 1...k-1, all of which appear
 175    at least once. x[n] is 0. (This is the corresponding output of transform.) k
 176    must be at most n+1. p is array of size n+1 whose contents are disregarded.
 177
 178    Output: x is V and p is I after the initial sorting stage of the refined
 179    suffix sorting algorithm.*/
 180
 181 void bucketsort(T *x, T *p, T n, T k)
 182 {
 183    T *pi, i, c, d, g;
 184
 185    for (pi=p; pi<p+k; ++pi)
 186       *pi=-1;                   /* mark linked lists empty.*/
 187    for (i=0; i<=n; ++i) {
 188       x[i]=p[c=x[i]];           /* insert in linked list.*/
 189       p[c]=i;
 190    }
 191    for (pi=p+k-1, i=n; pi>=p; --pi) {
 192       d=x[c=*pi];               /* c is position, d is next in list.*/
 193       x[c]=g=i;                 /* last position equals group number.*/
 194       if (d>=0) {               /* if more than one element in group.*/
 195          p[i--]=c;              /* p is permutation for the sorted x.*/
 196          do {
 197             d=x[c=d];           /* next in linked list.*/
 198             x[c]=g;             /* group number in x.*/
 199             p[i--]=c;           /* permutation in p.*/
 200          } while (d>=0);
 201       } else
 202          p[i--]=-1;             /* one element, sorted group.*/
 203    }
 204 }
 205
 206 /* Transforms the alphabet of x by attempting to aggregate several symbols into
 207    one, while preserving the suffix order of x. The alphabet may also be
 208    compacted, so that x on output comprises all integers of the new alphabet
 209    with no skipped numbers.
 210
 211    Input: x is an array of size n+1 whose first n elements are positive
 212    integers in the range l...k-1. p is array of size n+1, used for temporary
 213    storage. q controls aggregation and compaction by defining the maximum value
 214    for any symbol during transformation: q must be at least k-l; if q<=n,
 215    compaction is guaranteed; if k-l>n, compaction is never done; if q is
 216    INT_MAX, the maximum number of symbols are aggregated into one.
 217
 218    Output: Returns an integer j in the range 1...q representing the size of the
 219    new alphabet. If j<=n+1, the alphabet is compacted. The global variable r is
 220    set to the number of old symbols grouped into one. Only x[n] is 0.*/
 221
 222 T transform(T *x, T *p, T n, T k, T l, T q)
 223 {
 224    T b, c, d, e, i, j, m, s;
 225    T *pi, *pj;
 226
 227    for (s=0, i=k-l; i; i>>=1)
 228       ++s;                      /* s is number of bits in old symbol.*/
 229    e=INT_MAX>>s;                /* e is for overflow checking.*/
 230    for (b=d=r=0; r<n && d<=e && (c=d<<s|(k-l))<=q; ++r) {
 231       b=b<<s|(x[r]-l+1);        /* b is start of x in chunk alphabet.*/
 232       d=c;                      /* d is max symbol in chunk alphabet.*/
 233    }
 234    m=(1<<(r-1)*s)-1;            /* m masks off top old symbol from chunk.*/
 235    x[n]=l-1;                    /* emulate zero terminator.*/
 236    if (d<=n) {                  /* if bucketing possible, compact alphabet.*/
 237       for (pi=p; pi<=p+d; ++pi)
 238          *pi=0;                 /* zero transformation table.*/
 239       for (pi=x+r, c=b; pi<=x+n; ++pi) {
 240          p[c]=1;                /* mark used chunk symbol.*/
 241          c=(c&m)<<s|(*pi-l+1);  /* shift in next old symbol in chunk.*/
 242       }
 243       for (i=1; i<r; ++i) {     /* handle last r-1 positions.*/
 244          p[c]=1;                /* mark used chunk symbol.*/
 245          c=(c&m)<<s;            /* shift in next old symbol in chunk.*/
 246       }
 247       for (pi=p, j=1; pi<=p+d; ++pi)
 248          if (*pi)
 249             *pi=j++;            /* j is new alphabet size.*/
 250       for (pi=x, pj=x+r, c=b; pj<=x+n; ++pi, ++pj) {
 251          *pi=p[c];              /* transform to new alphabet.*/
 252          c=(c&m)<<s|(*pj-l+1);  /* shift in next old symbol in chunk.*/
 253       }
 254       while (pi<x+n) {          /* handle last r-1 positions.*/
 255          *pi++=p[c];            /* transform to new alphabet.*/
 256          c=(c&m)<<s;            /* shift right-end zero in chunk.*/
 257       }
 258    } else {                     /* bucketing not possible, don't compact.*/
 259       for (pi=x, pj=x+r, c=b; pj<=x+n; ++pi, ++pj) {
 260          *pi=c;                 /* transform to new alphabet.*/
 261          c=(c&m)<<s|(*pj-l+1);  /* shift in next old symbol in chunk.*/
 262       }
 263       while (pi<x+n) {          /* handle last r-1 positions.*/
 264          *pi++=c;               /* transform to new alphabet.*/
 265          c=(c&m)<<s;            /* shift right-end zero in chunk.*/
 266       }
 267       j=d+1;                    /* new alphabet size.*/
 268    }
 269    x[n]=0;                      /* end-of-string symbol is zero.*/
 270    return j;                    /* return new alphabet size.*/
 271 }
 272
 273 /* Makes suffix array p of x. x becomes inverse of p. p and x are both of size
 274    n+1. Contents of x[0...n-1] are integers in the range l...k-1. Original
 275    contents of x[n] is disregarded, the n-th symbol being regarded as
 276    end-of-string smaller than all other symbols.*/
 277
 278 void suffixsort(T *x, T *p, T n, T k, T l)
 279 {
 280    T *pi, *pk;
 281    T i, j, s, sl;
 282
 283    V=x;                         /* set global values.*/
 284    I=p;
 285
 286    if (n>=k-l) {                /* if bucketing possible,*/
 287       j=transform(V, I, n, k, l, n);
 288       bucketsort(V, I, n, j);   /* bucketsort on first r positions.*/
 289    } else {
 290       transform(V, I, n, k, l, INT_MAX);
 291       for (i=0; i<=n; ++i)
 292          I[i]=i;                /* initialize I with suffix numbers.*/
 293       h=0;
 294       sort_split(I, n+1);       /* quicksort on first r positions.*/
 295    }
 296    h=r;                         /* number of symbols aggregated by transform.*/
 297
 298    while (*I>=-n) {
 299       pi=I;                     /* pi is first position of group.*/
 300       sl=0;                     /* sl is negated length of sorted groups.*/
 301       do {
 302          if ((s=*pi)<0) {
 303             pi-=s;              /* skip over sorted group.*/
 304             sl+=s;              /* add negated length to sl.*/
 305          } else {
 306             if (sl) {
 307                *(pi+sl)=sl;     /* combine sorted groups before pi.*/
 308                sl=0;
 309             }
 310             pk=I+V[s]+1;        /* pk-1 is last position of unsorted group.*/
 311             sort_split(pi, pk-pi);
 312             pi=pk;              /* next group.*/
 313          }
 314       } while (pi<=I+n);
 315       if (sl)                   /* if the array ends with a sorted group.*/
 316          *(pi+sl)=sl;           /* combine sorted groups at end of I.*/
 317       h=2*h;                    /* double sorted-depth.*/
 318    }
 319
 320    for (i=0; i<=n; ++i)         /* reconstruct suffix array from inverse.*/
 321       I[V[i]]=i;
 322 }
 323
 324 };  // QSufSort
 325
 326 //--------------------------------------------------------------------------
 327
 328
 329 void suffixsort(int *x, int *p, int n, int k, int l)
 330 {
 331   QSufSort<int> sorter;
 332   sorter.suffixsort(x, p, n, k, l);
 333 }
 334
 335 void massive_suffixsort(sint *x, sint *p, sint n, sint k, sint l)
 336 {
 337   QSufSort<sint> sorter;
 338   sorter.suffixsort(x, p, n, k, l);
 339 }
 340
 341
 342 } // namespace CSA