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