annotate bwa-0.6.2/QSufSort.c @ 2:a294fbfcb1db draft default tip

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