2 Copyright 1999, N. Jesper Larsson, all rights reserved.
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).
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
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.
19 - Jouni Siren 2009-03-20
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)))
35 //--------------------------------------------------------------------------
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.*/
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.*/
50 void update_group(T *pl, T *pm)
54 g=pm-I; /* group number.*/
55 V[*pl]=g; /* update group number of first position.*/
57 *pl=-1; /* one element, sorted group.*/
59 do /* more than one element, unsorted group.*/
60 V[*++pl]=g; /* update group numbers.*/
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.*/
67 void select_sort_split(T *p, T n) {
71 pa=p; /* pa is start of group being picked out.*/
72 pn=p+n-1; /* pn is last position of subarray.*/
74 for (pi=pb=pa+1, f=KEY(pa); pi<=pn; ++pi)
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.*/
83 update_group(pa, pb-1); /* update group values for new group.*/
84 pa=pb; /* continue sorting rest of the subarray.*/
86 if (pa==pn) { /* check if last part is single element.*/
88 *pa=-1; /* sorted group.*/
92 /* Subroutine for sort_split, algorithm by Bentley & McIlroy.*/
94 T choose_pivot(T *p, T n) {
98 pm=p+(n>>1); /* small arrays, middle element.*/
102 if (n>40) { /* big arrays, pseudomedian of 9.*/
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);
108 pm=MED3(pl, pm, pn); /* midsize arrays, median of 3.*/
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.*/
119 void sort_split(T *p, T n)
121 T *pa, *pb, *pc, *pd, *pl, *pm, *pn;
124 if (n<7) { /* multi-selection sort smallest arrays.*/
125 select_sort_split(p, n);
129 v=choose_pivot(p, n);
132 while (1) { /* split-end partition.*/
133 while (pb<=pc && (f=KEY(pb))<=v) {
140 while (pc>=pb && (f=KEY(pc))>=v) {
154 if ((s=pa-p)>(t=pb-pa))
156 for (pl=p, pm=pb-s; s; --s, ++pl, ++pm)
158 if ((s=pd-pc)>(t=pn-pd-1))
160 for (pl=pb, pm=pn-s; s; --s, ++pl, ++pm)
167 update_group(p+s, p+n-t-1);
169 sort_split(p+n-t, t);
172 /* Bucketsort for first iteration.
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.
178 Output: x is V and p is I after the initial sorting stage of the refined
179 suffix sorting algorithm.*/
181 void bucketsort(T *x, T *p, T n, T k)
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.*/
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.*/
197 d=x[c=d]; /* next in linked list.*/
198 x[c]=g; /* group number in x.*/
199 p[i--]=c; /* permutation in p.*/
202 p[i--]=-1; /* one element, sorted group.*/
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.
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.
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.*/
222 T transform(T *x, T *p, T n, T k, T l, T q)
224 T b, c, d, e, i, j, m, s;
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.*/
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.*/
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.*/
247 for (pi=p, j=1; pi<=p+d; ++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.*/
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.*/
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.*/
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.*/
267 j=d+1; /* new alphabet size.*/
269 x[n]=0; /* end-of-string symbol is zero.*/
270 return j; /* return new alphabet size.*/
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.*/
278 void suffixsort(T *x, T *p, T n, T k, T l)
283 V=x; /* set global values.*/
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.*/
290 transform(V, I, n, k, l, INT_MAX);
292 I[i]=i; /* initialize I with suffix numbers.*/
294 sort_split(I, n+1); /* quicksort on first r positions.*/
296 h=r; /* number of symbols aggregated by transform.*/
299 pi=I; /* pi is first position of group.*/
300 sl=0; /* sl is negated length of sorted groups.*/
303 pi-=s; /* skip over sorted group.*/
304 sl+=s; /* add negated length to sl.*/
307 *(pi+sl)=sl; /* combine sorted groups before pi.*/
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.*/
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.*/
320 for (i=0; i<=n; ++i) /* reconstruct suffix array from inverse.*/
326 //--------------------------------------------------------------------------
329 void suffixsort(int *x, int *p, int n, int k, int l)
331 QSufSort<int> sorter;
332 sorter.suffixsort(x, p, n, k, l);
335 void massive_suffixsort(sint *x, sint *p, sint n, sint k, sint l)
337 QSufSort<sint> sorter;
338 sorter.suffixsort(x, p, n, k, l);