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 Replaced int -> sint for massive data support and moved into namespace CSA.
17 - Jouni Siren 2009-03-16
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.*/
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)))
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.*/
40 static void update_group(sint *pl, sint *pm)
44 g=pm-I; /* group number.*/
45 V[*pl]=g; /* update group number of first position.*/
47 *pl=-1; /* one element, sorted group.*/
49 do /* more than one element, unsorted group.*/
50 V[*++pl]=g; /* update group numbers.*/
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.*/
57 static void select_sort_split(sint *p, sint n) {
58 sint *pa, *pb, *pi, *pn;
61 pa=p; /* pa is start of group being picked out.*/
62 pn=p+n-1; /* pn is last position of subarray.*/
64 for (pi=pb=pa+1, f=KEY(pa); pi<=pn; ++pi)
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.*/
73 update_group(pa, pb-1); /* update group values for new group.*/
74 pa=pb; /* continue sorting rest of the subarray.*/
76 if (pa==pn) { /* check if last part is single element.*/
78 *pa=-1; /* sorted group.*/
82 /* Subroutine for sort_split, algorithm by Bentley & McIlroy.*/
84 static sint choose_pivot(sint *p, sint n) {
88 pm=p+(n>>1); /* small arrays, middle element.*/
92 if (n>40) { /* big arrays, pseudomedian of 9.*/
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);
98 pm=MED3(pl, pm, pn); /* midsize arrays, median of 3.*/
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.*/
109 static void sort_split(sint *p, sint n)
111 sint *pa, *pb, *pc, *pd, *pl, *pm, *pn;
112 sint f, v, s, t, tmp;
114 if (n<7) { /* multi-selection sort smallest arrays.*/
115 select_sort_split(p, n);
119 v=choose_pivot(p, n);
122 while (1) { /* split-end partition.*/
123 while (pb<=pc && (f=KEY(pb))<=v) {
130 while (pc>=pb && (f=KEY(pc))>=v) {
144 if ((s=pa-p)>(t=pb-pa))
146 for (pl=p, pm=pb-s; s; --s, ++pl, ++pm)
148 if ((s=pd-pc)>(t=pn-pd-1))
150 for (pl=pb, pm=pn-s; s; --s, ++pl, ++pm)
157 update_group(p+s, p+n-t-1);
159 sort_split(p+n-t, t);
162 /* Bucketsort for first iteration.
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.
168 Output: x is V and p is I after the initial sorting stage of the refined
169 suffix sorting algorithm.*/
171 static void bucketsort(sint *x, sint *p, sint n, sint k)
173 sint *pi, i, c, d, g;
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.*/
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.*/
187 d=x[c=d]; /* next in linked list.*/
188 x[c]=g; /* group number in x.*/
189 p[i--]=c; /* permutation in p.*/
192 p[i--]=-1; /* one element, sorted group.*/
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.
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.
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.*/
212 static sint transform(sint *x, sint *p, sint n, sint k, sint l, sint q)
214 sint b, c, d, e, i, j, m, s;
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.*/
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.*/
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.*/
237 for (pi=p, j=1; pi<=p+d; ++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.*/
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.*/
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.*/
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.*/
257 j=d+1; /* new alphabet size.*/
259 x[n]=0; /* end-of-string symbol is zero.*/
260 return j; /* return new alphabet size.*/
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.*/
268 void suffixsort(sint *x, sint *p, sint n, sint k, sint l)
273 V=x; /* set global values.*/
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.*/
280 transform(V, I, n, k, l, INT_MAX);
282 I[i]=i; /* initialize I with suffix numbers.*/
284 sort_split(I, n+1); /* quicksort on first r positions.*/
286 h=r; /* number of symbols aggregated by transform.*/
289 pi=I; /* pi is first position of group.*/
290 sl=0; /* sl is negated length of sorted groups.*/
293 pi-=s; /* skip over sorted group.*/
294 sl+=s; /* add negated length to sl.*/
297 *(pi+sl)=sl; /* combine sorted groups before pi.*/
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.*/
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.*/
310 for (i=0; i<=n; ++i) /* reconstruct suffix array from inverse.*/