|
0
|
1 /* QSufSort.c
|
|
|
2
|
|
|
3 Original source from qsufsort.c
|
|
|
4
|
|
|
5 Copyright 1999, N. Jesper Larsson, all rights reserved.
|
|
|
6
|
|
|
7 This file contains an implementation of the algorithm presented in "Faster
|
|
|
8 Suffix Sorting" by N. Jesper Larsson (jesper@cs.lth.se) and Kunihiko
|
|
|
9 Sadakane (sada@is.s.u-tokyo.ac.jp).
|
|
|
10
|
|
|
11 This software may be used freely for any purpose. However, when distributed,
|
|
|
12 the original source must be clearly stated, and, when the source code is
|
|
|
13 distributed, the copyright notice must be retained and any alterations in
|
|
|
14 the code must be clearly marked. No warranty is given regarding the quality
|
|
|
15 of this software.
|
|
|
16
|
|
|
17 Modified by Wong Chi-Kwong, 2004
|
|
|
18
|
|
|
19 Changes summary: - Used long variable and function names
|
|
|
20 - Removed global variables
|
|
|
21 - Replace pointer references with array references
|
|
|
22 - Used insertion sort in place of selection sort and increased insertion sort threshold
|
|
|
23 - Reconstructing suffix array from inverse becomes an option
|
|
|
24 - Add handling where end-of-text symbol is not necessary < all characters
|
|
|
25 - Removed codes for supporting alphabet size > number of characters
|
|
|
26
|
|
|
27 No warrenty is given regarding the quality of the modifications.
|
|
|
28
|
|
|
29 */
|
|
|
30
|
|
|
31
|
|
|
32 #include <stdio.h>
|
|
|
33 #include <stdlib.h>
|
|
|
34 #include <limits.h>
|
|
|
35 #include "QSufSort.h"
|
|
|
36
|
|
|
37 #define min(value1, value2) ( ((value1) < (value2)) ? (value1) : (value2) )
|
|
|
38 #define med3(a, b, c) ( a<b ? (b<c ? b : a<c ? c : a) : (b>c ? b : a>c ? c : a))
|
|
|
39 #define swap(a, b, t); t = a; a = b; b = t;
|
|
|
40
|
|
|
41 // Static functions
|
|
|
42 static void QSufSortSortSplit(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
43 const qsint_t highestPos, const qsint_t numSortedChar);
|
|
|
44 static qsint_t QSufSortChoosePivot(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
45 const qsint_t highestPos, const qsint_t numSortedChar);
|
|
|
46 static void QSufSortInsertSortSplit(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
47 const qsint_t highestPos, const qsint_t numSortedChar);
|
|
|
48 static void QSufSortBucketSort(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t alphabetSize);
|
|
|
49 static qsint_t QSufSortTransform(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t largestInputSymbol,
|
|
|
50 const qsint_t smallestInputSymbol, const qsint_t maxNewAlphabetSize, qsint_t *numSymbolAggregated);
|
|
|
51
|
|
|
52 /* Makes suffix array p of x. x becomes inverse of p. p and x are both of size
|
|
|
53 n+1. Contents of x[0...n-1] are integers in the range l...k-1. Original
|
|
|
54 contents of x[n] is disregarded, the n-th symbol being regarded as
|
|
|
55 end-of-string smaller than all other symbols.*/
|
|
|
56 void QSufSortSuffixSort(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t largestInputSymbol,
|
|
|
57 const qsint_t smallestInputSymbol, const int skipTransform)
|
|
|
58 {
|
|
|
59 qsint_t i, j;
|
|
|
60 qsint_t s, negatedSortedGroupLength;
|
|
|
61 qsint_t numSymbolAggregated;
|
|
|
62 qsint_t numSortedPos = 1;
|
|
|
63 qsint_t newAlphabetSize;
|
|
|
64
|
|
|
65 if (!skipTransform) {
|
|
|
66 /* bucketing possible*/
|
|
|
67 newAlphabetSize = QSufSortTransform(V, I, numChar, largestInputSymbol, smallestInputSymbol,
|
|
|
68 numChar, &numSymbolAggregated);
|
|
|
69 QSufSortBucketSort(V, I, numChar, newAlphabetSize);
|
|
|
70 I[0] = -1;
|
|
|
71 V[numChar] = 0;
|
|
|
72 numSortedPos = numSymbolAggregated;
|
|
|
73 }
|
|
|
74
|
|
|
75 while ((qsint_t)(I[0]) >= -(qsint_t)numChar) {
|
|
|
76 i = 0;
|
|
|
77 negatedSortedGroupLength = 0;
|
|
|
78 do {
|
|
|
79 s = I[i];
|
|
|
80 if (s < 0) {
|
|
|
81 i -= s; /* skip over sorted group.*/
|
|
|
82 negatedSortedGroupLength += s;
|
|
|
83 } else {
|
|
|
84 if (negatedSortedGroupLength) {
|
|
|
85 I[i+negatedSortedGroupLength] = negatedSortedGroupLength; /* combine preceding sorted groups */
|
|
|
86 negatedSortedGroupLength = 0;
|
|
|
87 }
|
|
|
88 j = V[s] + 1;
|
|
|
89 QSufSortSortSplit(V, I, i, j - 1, numSortedPos);
|
|
|
90 i = j;
|
|
|
91 }
|
|
|
92 } while (i <= numChar);
|
|
|
93 if (negatedSortedGroupLength) {
|
|
|
94 /* array ends with a sorted group.*/
|
|
|
95 I[i+negatedSortedGroupLength] = negatedSortedGroupLength; /* combine sorted groups at end of I.*/
|
|
|
96 }
|
|
|
97 numSortedPos *= 2; /* double sorted-depth.*/
|
|
|
98 }
|
|
|
99 }
|
|
|
100
|
|
|
101 void QSufSortGenerateSaFromInverse(const qsint_t* V, qsint_t* __restrict I, const qsint_t numChar)
|
|
|
102 {
|
|
|
103 qsint_t i;
|
|
|
104 for (i=0; i<=numChar; i++)
|
|
|
105 I[V[i]] = i + 1;
|
|
|
106 }
|
|
|
107
|
|
|
108 /* Sorting routine called for each unsorted group. Sorts the array of integers
|
|
|
109 (suffix numbers) of length n starting at p. The algorithm is a ternary-split
|
|
|
110 quicksort taken from Bentley & McIlroy, "Engineering a Sort Function",
|
|
|
111 Software -- Practice and Experience 23(11), 1249-1265 (November 1993). This
|
|
|
112 function is based on Program 7.*/
|
|
|
113 static void QSufSortSortSplit(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
114 const qsint_t highestPos, const qsint_t numSortedChar) {
|
|
|
115
|
|
|
116 qsint_t a, b, c, d;
|
|
|
117 qsint_t l, m;
|
|
|
118 qsint_t f, v, s, t;
|
|
|
119 qsint_t tmp;
|
|
|
120 qsint_t numItem;
|
|
|
121
|
|
|
122 numItem = highestPos - lowestPos + 1;
|
|
|
123
|
|
|
124 if (numItem <= INSERT_SORT_NUM_ITEM) {
|
|
|
125 QSufSortInsertSortSplit(V, I, lowestPos, highestPos, numSortedChar);
|
|
|
126 return;
|
|
|
127 }
|
|
|
128
|
|
|
129 v = QSufSortChoosePivot(V, I, lowestPos, highestPos, numSortedChar);
|
|
|
130
|
|
|
131 a = b = lowestPos;
|
|
|
132 c = d = highestPos;
|
|
|
133
|
|
|
134 while (1) {
|
|
|
135 while (c >= b && (f = KEY(V, I, b, numSortedChar)) <= v) {
|
|
|
136 if (f == v) {
|
|
|
137 swap(I[a], I[b], tmp);
|
|
|
138 a++;
|
|
|
139 }
|
|
|
140 b++;
|
|
|
141 }
|
|
|
142 while (c >= b && (f = KEY(V, I, c, numSortedChar)) >= v) {
|
|
|
143 if (f == v) {
|
|
|
144 swap(I[c], I[d], tmp);
|
|
|
145 d--;
|
|
|
146 }
|
|
|
147 c--;
|
|
|
148 }
|
|
|
149 if (b > c)
|
|
|
150 break;
|
|
|
151 swap(I[b], I[c], tmp);
|
|
|
152 b++;
|
|
|
153 c--;
|
|
|
154 }
|
|
|
155
|
|
|
156 s = a - lowestPos;
|
|
|
157 t = b - a;
|
|
|
158 s = min(s, t);
|
|
|
159 for (l = lowestPos, m = b - s; m < b; l++, m++) {
|
|
|
160 swap(I[l], I[m], tmp);
|
|
|
161 }
|
|
|
162
|
|
|
163 s = d - c;
|
|
|
164 t = highestPos - d;
|
|
|
165 s = min(s, t);
|
|
|
166 for (l = b, m = highestPos - s + 1; m <= highestPos; l++, m++) {
|
|
|
167 swap(I[l], I[m], tmp);
|
|
|
168 }
|
|
|
169
|
|
|
170 s = b - a;
|
|
|
171 t = d - c;
|
|
|
172 if (s > 0)
|
|
|
173 QSufSortSortSplit(V, I, lowestPos, lowestPos + s - 1, numSortedChar);
|
|
|
174
|
|
|
175 // Update group number for equal portion
|
|
|
176 a = lowestPos + s;
|
|
|
177 b = highestPos - t;
|
|
|
178 if (a == b) {
|
|
|
179 // Sorted group
|
|
|
180 V[I[a]] = a;
|
|
|
181 I[a] = -1;
|
|
|
182 } else {
|
|
|
183 // Unsorted group
|
|
|
184 for (c=a; c<=b; c++)
|
|
|
185 V[I[c]] = b;
|
|
|
186 }
|
|
|
187
|
|
|
188 if (t > 0)
|
|
|
189 QSufSortSortSplit(V, I, highestPos - t + 1, highestPos, numSortedChar);
|
|
|
190
|
|
|
191 }
|
|
|
192
|
|
|
193 /* Algorithm by Bentley & McIlroy.*/
|
|
|
194 static qsint_t QSufSortChoosePivot(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
195 const qsint_t highestPos, const qsint_t numSortedChar) {
|
|
|
196
|
|
|
197 qsint_t m;
|
|
|
198 qsint_t keyl, keym, keyn;
|
|
|
199 qsint_t key1, key2, key3;
|
|
|
200 qsint_t s;
|
|
|
201 qsint_t numItem;
|
|
|
202
|
|
|
203 numItem = highestPos - lowestPos + 1;
|
|
|
204
|
|
|
205 m = lowestPos + numItem / 2;
|
|
|
206
|
|
|
207 s = numItem / 8;
|
|
|
208 key1 = KEY(V, I, lowestPos, numSortedChar);
|
|
|
209 key2 = KEY(V, I, lowestPos+s, numSortedChar);
|
|
|
210 key3 = KEY(V, I, lowestPos+2*s, numSortedChar);
|
|
|
211 keyl = med3(key1, key2, key3);
|
|
|
212 key1 = KEY(V, I, m-s, numSortedChar);
|
|
|
213 key2 = KEY(V, I, m, numSortedChar);
|
|
|
214 key3 = KEY(V, I, m+s, numSortedChar);
|
|
|
215 keym = med3(key1, key2, key3);
|
|
|
216 key1 = KEY(V, I, highestPos-2*s, numSortedChar);
|
|
|
217 key2 = KEY(V, I, highestPos-s, numSortedChar);
|
|
|
218 key3 = KEY(V, I, highestPos, numSortedChar);
|
|
|
219 keyn = med3(key1, key2, key3);
|
|
|
220
|
|
|
221 return med3(keyl, keym, keyn);
|
|
|
222
|
|
|
223
|
|
|
224 }
|
|
|
225
|
|
|
226 /* Quadratic sorting method to use for small subarrays. */
|
|
|
227 static void QSufSortInsertSortSplit(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t lowestPos,
|
|
|
228 const qsint_t highestPos, const qsint_t numSortedChar)
|
|
|
229 {
|
|
|
230 qsint_t i, j;
|
|
|
231 qsint_t tmpKey, tmpPos;
|
|
|
232 qsint_t numItem;
|
|
|
233 qsint_t key[INSERT_SORT_NUM_ITEM], pos[INSERT_SORT_NUM_ITEM];
|
|
|
234 qsint_t negativeSortedLength;
|
|
|
235 qsint_t groupNum;
|
|
|
236
|
|
|
237 numItem = highestPos - lowestPos + 1;
|
|
|
238
|
|
|
239 for (i=0; i<numItem; i++) {
|
|
|
240 pos[i] = I[lowestPos + i];
|
|
|
241 key[i] = V[pos[i] + numSortedChar];
|
|
|
242 }
|
|
|
243
|
|
|
244 for (i=1; i<numItem; i++) {
|
|
|
245 tmpKey = key[i];
|
|
|
246 tmpPos = pos[i];
|
|
|
247 for (j=i; j>0 && key[j-1] > tmpKey; j--) {
|
|
|
248 key[j] = key[j-1];
|
|
|
249 pos[j] = pos[j-1];
|
|
|
250 }
|
|
|
251 key[j] = tmpKey;
|
|
|
252 pos[j] = tmpPos;
|
|
|
253 }
|
|
|
254
|
|
|
255 negativeSortedLength = -1;
|
|
|
256
|
|
|
257 i = numItem - 1;
|
|
|
258 groupNum = highestPos;
|
|
|
259 while (i > 0) {
|
|
|
260 I[i+lowestPos] = pos[i];
|
|
|
261 V[I[i+lowestPos]] = groupNum;
|
|
|
262 if (key[i-1] == key[i]) {
|
|
|
263 negativeSortedLength = 0;
|
|
|
264 } else {
|
|
|
265 if (negativeSortedLength < 0)
|
|
|
266 I[i+lowestPos] = negativeSortedLength;
|
|
|
267 groupNum = i + lowestPos - 1;
|
|
|
268 negativeSortedLength--;
|
|
|
269 }
|
|
|
270 i--;
|
|
|
271 }
|
|
|
272
|
|
|
273 I[lowestPos] = pos[0];
|
|
|
274 V[I[lowestPos]] = groupNum;
|
|
|
275 if (negativeSortedLength < 0)
|
|
|
276 I[lowestPos] = negativeSortedLength;
|
|
|
277 }
|
|
|
278
|
|
|
279 /* Bucketsort for first iteration.
|
|
|
280
|
|
|
281 Input: x[0...n-1] holds integers in the range 1...k-1, all of which appear
|
|
|
282 at least once. x[n] is 0. (This is the corresponding output of transform.) k
|
|
|
283 must be at most n+1. p is array of size n+1 whose contents are disregarded.
|
|
|
284
|
|
|
285 Output: x is V and p is I after the initial sorting stage of the refined
|
|
|
286 suffix sorting algorithm.*/
|
|
|
287
|
|
|
288 static void QSufSortBucketSort(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t alphabetSize)
|
|
|
289 {
|
|
|
290 qsint_t i, c;
|
|
|
291 qsint_t d;
|
|
|
292 qsint_t groupNum;
|
|
|
293 qsint_t currentIndex;
|
|
|
294
|
|
|
295 // mark linked list empty
|
|
|
296 for (i=0; i<alphabetSize; i++)
|
|
|
297 I[i] = -1;
|
|
|
298
|
|
|
299 // insert to linked list
|
|
|
300 for (i=0; i<=numChar; i++) {
|
|
|
301 c = V[i];
|
|
|
302 V[i] = (qsint_t)(I[c]);
|
|
|
303 I[c] = i;
|
|
|
304 }
|
|
|
305
|
|
|
306 currentIndex = numChar;
|
|
|
307 for (i=alphabetSize; i>0; i--) {
|
|
|
308 c = I[i-1];
|
|
|
309 d = (qsint_t)(V[c]);
|
|
|
310 groupNum = currentIndex;
|
|
|
311 V[c] = groupNum;
|
|
|
312 if (d >= 0) {
|
|
|
313 I[currentIndex] = c;
|
|
|
314 while (d >= 0) {
|
|
|
315 c = d;
|
|
|
316 d = V[c];
|
|
|
317 V[c] = groupNum;
|
|
|
318 currentIndex--;
|
|
|
319 I[currentIndex] = c;
|
|
|
320 }
|
|
|
321 } else {
|
|
|
322 // sorted group
|
|
|
323 I[currentIndex] = -1;
|
|
|
324 }
|
|
|
325 currentIndex--;
|
|
|
326 }
|
|
|
327 }
|
|
|
328
|
|
|
329 /* Transforms the alphabet of x by attempting to aggregate several symbols into
|
|
|
330 one, while preserving the suffix order of x. The alphabet may also be
|
|
|
331 compacted, so that x on output comprises all integers of the new alphabet
|
|
|
332 with no skipped numbers.
|
|
|
333
|
|
|
334 Input: x is an array of size n+1 whose first n elements are positive
|
|
|
335 integers in the range l...k-1. p is array of size n+1, used for temporary
|
|
|
336 storage. q controls aggregation and compaction by defining the maximum intue
|
|
|
337 for any symbol during transformation: q must be at least k-l; if q<=n,
|
|
|
338 compaction is guaranteed; if k-l>n, compaction is never done; if q is
|
|
|
339 INT_MAX, the maximum number of symbols are aggregated into one.
|
|
|
340
|
|
|
341 Output: Returns an integer j in the range 1...q representing the size of the
|
|
|
342 new alphabet. If j<=n+1, the alphabet is compacted. The global variable r is
|
|
|
343 set to the number of old symbols grouped into one. Only x[n] is 0.*/
|
|
|
344 static qsint_t QSufSortTransform(qsint_t* __restrict V, qsint_t* __restrict I, const qsint_t numChar, const qsint_t largestInputSymbol,
|
|
|
345 const qsint_t smallestInputSymbol, const qsint_t maxNewAlphabetSize, qsint_t *numSymbolAggregated)
|
|
|
346 {
|
|
|
347 qsint_t c, i, j;
|
|
|
348 qsint_t a; // numSymbolAggregated
|
|
|
349 qsint_t mask;
|
|
|
350 qsint_t minSymbolInChunk = 0, maxSymbolInChunk = 0;
|
|
|
351 qsint_t newAlphabetSize;
|
|
|
352 qsint_t maxNumInputSymbol, maxNumBit, maxSymbol;
|
|
|
353
|
|
|
354 maxNumInputSymbol = largestInputSymbol - smallestInputSymbol + 1;
|
|
|
355
|
|
|
356 for (maxNumBit = 0, i = maxNumInputSymbol; i; i >>= 1) ++maxNumBit;
|
|
|
357 maxSymbol = QSINT_MAX >> maxNumBit;
|
|
|
358
|
|
|
359 c = maxNumInputSymbol;
|
|
|
360 for (a = 0; a < numChar && maxSymbolInChunk <= maxSymbol && c <= maxNewAlphabetSize; a++) {
|
|
|
361 minSymbolInChunk = (minSymbolInChunk << maxNumBit) | (V[a] - smallestInputSymbol + 1);
|
|
|
362 maxSymbolInChunk = c;
|
|
|
363 c = (maxSymbolInChunk << maxNumBit) | maxNumInputSymbol;
|
|
|
364 }
|
|
|
365
|
|
|
366 mask = (1 << (a-1) * maxNumBit) - 1; /* mask masks off top old symbol from chunk.*/
|
|
|
367 V[numChar] = smallestInputSymbol - 1; /* emulate zero terminator.*/
|
|
|
368
|
|
|
369 /* bucketing possible, compact alphabet.*/
|
|
|
370 for (i=0; i<=maxSymbolInChunk; i++)
|
|
|
371 I[i] = 0; /* zero transformation table.*/
|
|
|
372 c = minSymbolInChunk;
|
|
|
373 for (i=a; i<=numChar; i++) {
|
|
|
374 I[c] = 1; /* mark used chunk symbol.*/
|
|
|
375 c = ((c & mask) << maxNumBit) | (V[i] - smallestInputSymbol + 1); /* shift in next old symbol in chunk.*/
|
|
|
376 }
|
|
|
377 for (i=1; i<a; i++) { /* handle last r-1 positions.*/
|
|
|
378 I[c] = 1; /* mark used chunk symbol.*/
|
|
|
379 c = (c & mask) << maxNumBit; /* shift in next old symbol in chunk.*/
|
|
|
380 }
|
|
|
381 newAlphabetSize = 1;
|
|
|
382 for (i=0; i<=maxSymbolInChunk; i++) {
|
|
|
383 if (I[i]) {
|
|
|
384 I[i] = newAlphabetSize;
|
|
|
385 newAlphabetSize++;
|
|
|
386 }
|
|
|
387 }
|
|
|
388 c = minSymbolInChunk;
|
|
|
389 for (i=0, j=a; j<=numChar; i++, j++) {
|
|
|
390 V[i] = I[c]; /* transform to new alphabet.*/
|
|
|
391 c = ((c & mask) << maxNumBit) | (V[j] - smallestInputSymbol + 1); /* shift in next old symbol in chunk.*/
|
|
|
392 }
|
|
|
393 for (; i<numChar; i++) { /* handle last a-1 positions.*/
|
|
|
394 V[i] = I[c]; /* transform to new alphabet.*/
|
|
|
395 c = (c & mask) << maxNumBit; /* shift right-end zero in chunk.*/
|
|
|
396 }
|
|
|
397
|
|
|
398 V[numChar] = 0; /* end-of-string symbol is zero.*/
|
|
|
399
|
|
|
400 *numSymbolAggregated = a;
|
|
|
401 return newAlphabetSize;
|
|
|
402 }
|
