File: | home/bhubbard/working/src/ceph/src/zstd/lib/dictBuilder/divsufsort.c |
Warning: | line 1789, column 16 Dereference of null pointer |
[?] Use j/k keys for keyboard navigation
1 | /* | |||
2 | * divsufsort.c for libdivsufsort-lite | |||
3 | * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved. | |||
4 | * | |||
5 | * Permission is hereby granted, free of charge, to any person | |||
6 | * obtaining a copy of this software and associated documentation | |||
7 | * files (the "Software"), to deal in the Software without | |||
8 | * restriction, including without limitation the rights to use, | |||
9 | * copy, modify, merge, publish, distribute, sublicense, and/or sell | |||
10 | * copies of the Software, and to permit persons to whom the | |||
11 | * Software is furnished to do so, subject to the following | |||
12 | * conditions: | |||
13 | * | |||
14 | * The above copyright notice and this permission notice shall be | |||
15 | * included in all copies or substantial portions of the Software. | |||
16 | * | |||
17 | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, | |||
18 | * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES | |||
19 | * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND | |||
20 | * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT | |||
21 | * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, | |||
22 | * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |||
23 | * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR | |||
24 | * OTHER DEALINGS IN THE SOFTWARE. | |||
25 | */ | |||
26 | ||||
27 | /*- Compiler specifics -*/ | |||
28 | #ifdef __clang__1 | |||
29 | #pragma clang diagnostic ignored "-Wshorten-64-to-32" | |||
30 | #endif | |||
31 | ||||
32 | #if defined(_MSC_VER) | |||
33 | # pragma warning(disable : 4244) | |||
34 | # pragma warning(disable : 4127) /* C4127 : Condition expression is constant */ | |||
35 | #endif | |||
36 | ||||
37 | ||||
38 | /*- Dependencies -*/ | |||
39 | #include <assert.h> | |||
40 | #include <stdio.h> | |||
41 | #include <stdlib.h> | |||
42 | ||||
43 | #include "divsufsort.h" | |||
44 | ||||
45 | /*- Constants -*/ | |||
46 | #if defined(INLINE__inline) | |||
47 | # undef INLINE__inline | |||
48 | #endif | |||
49 | #if !defined(INLINE__inline) | |||
50 | # define INLINE__inline __inline | |||
51 | #endif | |||
52 | #if defined(ALPHABET_SIZE(256)) && (ALPHABET_SIZE(256) < 1) | |||
53 | # undef ALPHABET_SIZE(256) | |||
54 | #endif | |||
55 | #if !defined(ALPHABET_SIZE(256)) | |||
56 | # define ALPHABET_SIZE(256) (256) | |||
57 | #endif | |||
58 | #define BUCKET_A_SIZE((256)) (ALPHABET_SIZE(256)) | |||
59 | #define BUCKET_B_SIZE((256) * (256)) (ALPHABET_SIZE(256) * ALPHABET_SIZE(256)) | |||
60 | #if defined(SS_INSERTIONSORT_THRESHOLD(8)) | |||
61 | # if SS_INSERTIONSORT_THRESHOLD(8) < 1 | |||
62 | # undef SS_INSERTIONSORT_THRESHOLD(8) | |||
63 | # define SS_INSERTIONSORT_THRESHOLD(8) (1) | |||
64 | # endif | |||
65 | #else | |||
66 | # define SS_INSERTIONSORT_THRESHOLD(8) (8) | |||
67 | #endif | |||
68 | #if defined(SS_BLOCKSIZE(1024)) | |||
69 | # if SS_BLOCKSIZE(1024) < 0 | |||
70 | # undef SS_BLOCKSIZE(1024) | |||
71 | # define SS_BLOCKSIZE(1024) (0) | |||
72 | # elif 32768 <= SS_BLOCKSIZE(1024) | |||
73 | # undef SS_BLOCKSIZE(1024) | |||
74 | # define SS_BLOCKSIZE(1024) (32767) | |||
75 | # endif | |||
76 | #else | |||
77 | # define SS_BLOCKSIZE(1024) (1024) | |||
78 | #endif | |||
79 | /* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */ | |||
80 | #if SS_BLOCKSIZE(1024) == 0 | |||
81 | # define SS_MISORT_STACKSIZE(16) (96) | |||
82 | #elif SS_BLOCKSIZE(1024) <= 4096 | |||
83 | # define SS_MISORT_STACKSIZE(16) (16) | |||
84 | #else | |||
85 | # define SS_MISORT_STACKSIZE(16) (24) | |||
86 | #endif | |||
87 | #define SS_SMERGE_STACKSIZE(32) (32) | |||
88 | #define TR_INSERTIONSORT_THRESHOLD(8) (8) | |||
89 | #define TR_STACKSIZE(64) (64) | |||
90 | ||||
91 | ||||
92 | /*- Macros -*/ | |||
93 | #ifndef SWAP | |||
94 | # define SWAP(_a, _b)do { t = (_a); (_a) = (_b); (_b) = t; } while(0) do { t = (_a); (_a) = (_b); (_b) = t; } while(0) | |||
95 | #endif /* SWAP */ | |||
96 | #ifndef MIN | |||
97 | # define MIN(_a, _b)(((_a) < (_b)) ? (_a) : (_b)) (((_a) < (_b)) ? (_a) : (_b)) | |||
98 | #endif /* MIN */ | |||
99 | #ifndef MAX | |||
100 | # define MAX(_a, _b)(((_a) > (_b)) ? (_a) : (_b)) (((_a) > (_b)) ? (_a) : (_b)) | |||
101 | #endif /* MAX */ | |||
102 | #define STACK_PUSH(_a, _b, _c, _d)do { ((void) (0)); stack[ssize].a = (_a), stack[ssize].b = (_b ), stack[ssize].c = (_c), stack[ssize++].d = (_d); } while(0)\ | |||
103 | do {\ | |||
104 | assert(ssize < STACK_SIZE)((void) (0));\ | |||
105 | stack[ssize].a = (_a), stack[ssize].b = (_b),\ | |||
106 | stack[ssize].c = (_c), stack[ssize++].d = (_d);\ | |||
107 | } while(0) | |||
108 | #define STACK_PUSH5(_a, _b, _c, _d, _e)do { ((void) (0)); stack[ssize].a = (_a), stack[ssize].b = (_b ), stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++ ].e = (_e); } while(0)\ | |||
109 | do {\ | |||
110 | assert(ssize < STACK_SIZE)((void) (0));\ | |||
111 | stack[ssize].a = (_a), stack[ssize].b = (_b),\ | |||
112 | stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\ | |||
113 | } while(0) | |||
114 | #define STACK_POP(_a, _b, _c, _d)do { ((void) (0)); if(ssize == 0) { return; } (_a) = stack[-- ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d; } while(0)\ | |||
115 | do {\ | |||
116 | assert(0 <= ssize)((void) (0));\ | |||
117 | if(ssize == 0) { return; }\ | |||
118 | (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ | |||
119 | (_c) = stack[ssize].c, (_d) = stack[ssize].d;\ | |||
120 | } while(0) | |||
121 | #define STACK_POP5(_a, _b, _c, _d, _e)do { ((void) (0)); if(ssize == 0) { return; } (_a) = stack[-- ssize].a, (_b) = stack[ssize].b, (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e; } while(0)\ | |||
122 | do {\ | |||
123 | assert(0 <= ssize)((void) (0));\ | |||
124 | if(ssize == 0) { return; }\ | |||
125 | (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ | |||
126 | (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\ | |||
127 | } while(0) | |||
128 | #define BUCKET_A(_c0)bucket_A[(_c0)] bucket_A[(_c0)] | |||
129 | #if ALPHABET_SIZE(256) == 256 | |||
130 | #define BUCKET_B(_c0, _c1)(bucket_B[((_c1) << 8) | (_c0)]) (bucket_B[((_c1) << 8) | (_c0)]) | |||
131 | #define BUCKET_BSTAR(_c0, _c1)(bucket_B[((_c0) << 8) | (_c1)]) (bucket_B[((_c0) << 8) | (_c1)]) | |||
132 | #else | |||
133 | #define BUCKET_B(_c0, _c1)(bucket_B[((_c1) << 8) | (_c0)]) (bucket_B[(_c1) * ALPHABET_SIZE(256) + (_c0)]) | |||
134 | #define BUCKET_BSTAR(_c0, _c1)(bucket_B[((_c0) << 8) | (_c1)]) (bucket_B[(_c0) * ALPHABET_SIZE(256) + (_c1)]) | |||
135 | #endif | |||
136 | ||||
137 | ||||
138 | /*- Private Functions -*/ | |||
139 | ||||
140 | static const int lg_table[256]= { | |||
141 | -1,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, | |||
142 | 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, | |||
143 | 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, | |||
144 | 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, | |||
145 | 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, | |||
146 | 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, | |||
147 | 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, | |||
148 | 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7 | |||
149 | }; | |||
150 | ||||
151 | #if (SS_BLOCKSIZE(1024) == 0) || (SS_INSERTIONSORT_THRESHOLD(8) < SS_BLOCKSIZE(1024)) | |||
152 | ||||
153 | static INLINE__inline | |||
154 | int | |||
155 | ss_ilg(int n) { | |||
156 | #if SS_BLOCKSIZE(1024) == 0 | |||
157 | return (n & 0xffff0000) ? | |||
158 | ((n & 0xff000000) ? | |||
159 | 24 + lg_table[(n >> 24) & 0xff] : | |||
160 | 16 + lg_table[(n >> 16) & 0xff]) : | |||
161 | ((n & 0x0000ff00) ? | |||
162 | 8 + lg_table[(n >> 8) & 0xff] : | |||
163 | 0 + lg_table[(n >> 0) & 0xff]); | |||
164 | #elif SS_BLOCKSIZE(1024) < 256 | |||
165 | return lg_table[n]; | |||
166 | #else | |||
167 | return (n & 0xff00) ? | |||
168 | 8 + lg_table[(n >> 8) & 0xff] : | |||
169 | 0 + lg_table[(n >> 0) & 0xff]; | |||
170 | #endif | |||
171 | } | |||
172 | ||||
173 | #endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ | |||
174 | ||||
175 | #if SS_BLOCKSIZE(1024) != 0 | |||
176 | ||||
177 | static const int sqq_table[256] = { | |||
178 | 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, | |||
179 | 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, | |||
180 | 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109, | |||
181 | 110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, | |||
182 | 128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, | |||
183 | 143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, | |||
184 | 156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, | |||
185 | 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, | |||
186 | 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, | |||
187 | 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, | |||
188 | 202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211, | |||
189 | 212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, | |||
190 | 221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, | |||
191 | 230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, | |||
192 | 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, | |||
193 | 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255 | |||
194 | }; | |||
195 | ||||
196 | static INLINE__inline | |||
197 | int | |||
198 | ss_isqrt(int x) { | |||
199 | int y, e; | |||
200 | ||||
201 | if(x >= (SS_BLOCKSIZE(1024) * SS_BLOCKSIZE(1024))) { return SS_BLOCKSIZE(1024); } | |||
202 | e = (x & 0xffff0000) ? | |||
203 | ((x & 0xff000000) ? | |||
204 | 24 + lg_table[(x >> 24) & 0xff] : | |||
205 | 16 + lg_table[(x >> 16) & 0xff]) : | |||
206 | ((x & 0x0000ff00) ? | |||
207 | 8 + lg_table[(x >> 8) & 0xff] : | |||
208 | 0 + lg_table[(x >> 0) & 0xff]); | |||
209 | ||||
210 | if(e >= 16) { | |||
211 | y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7); | |||
212 | if(e >= 24) { y = (y + 1 + x / y) >> 1; } | |||
213 | y = (y + 1 + x / y) >> 1; | |||
214 | } else if(e >= 8) { | |||
215 | y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1; | |||
216 | } else { | |||
217 | return sqq_table[x] >> 4; | |||
218 | } | |||
219 | ||||
220 | return (x < (y * y)) ? y - 1 : y; | |||
221 | } | |||
222 | ||||
223 | #endif /* SS_BLOCKSIZE != 0 */ | |||
224 | ||||
225 | ||||
226 | /*---------------------------------------------------------------------------*/ | |||
227 | ||||
228 | /* Compares two suffixes. */ | |||
229 | static INLINE__inline | |||
230 | int | |||
231 | ss_compare(const unsigned char *T, | |||
232 | const int *p1, const int *p2, | |||
233 | int depth) { | |||
234 | const unsigned char *U1, *U2, *U1n, *U2n; | |||
235 | ||||
236 | for(U1 = T + depth + *p1, | |||
237 | U2 = T + depth + *p2, | |||
238 | U1n = T + *(p1 + 1) + 2, | |||
239 | U2n = T + *(p2 + 1) + 2; | |||
240 | (U1 < U1n) && (U2 < U2n) && (*U1 == *U2); | |||
241 | ++U1, ++U2) { | |||
242 | } | |||
243 | ||||
244 | return U1 < U1n ? | |||
245 | (U2 < U2n ? *U1 - *U2 : 1) : | |||
246 | (U2 < U2n ? -1 : 0); | |||
247 | } | |||
248 | ||||
249 | ||||
250 | /*---------------------------------------------------------------------------*/ | |||
251 | ||||
252 | #if (SS_BLOCKSIZE(1024) != 1) && (SS_INSERTIONSORT_THRESHOLD(8) != 1) | |||
253 | ||||
254 | /* Insertionsort for small size groups */ | |||
255 | static | |||
256 | void | |||
257 | ss_insertionsort(const unsigned char *T, const int *PA, | |||
258 | int *first, int *last, int depth) { | |||
259 | int *i, *j; | |||
260 | int t; | |||
261 | int r; | |||
262 | ||||
263 | for(i = last - 2; first <= i; --i) { | |||
264 | for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) { | |||
265 | do { *(j - 1) = *j; } while((++j < last) && (*j < 0)); | |||
266 | if(last <= j) { break; } | |||
267 | } | |||
268 | if(r == 0) { *j = ~*j; } | |||
269 | *(j - 1) = t; | |||
270 | } | |||
271 | } | |||
272 | ||||
273 | #endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */ | |||
274 | ||||
275 | ||||
276 | /*---------------------------------------------------------------------------*/ | |||
277 | ||||
278 | #if (SS_BLOCKSIZE(1024) == 0) || (SS_INSERTIONSORT_THRESHOLD(8) < SS_BLOCKSIZE(1024)) | |||
279 | ||||
280 | static INLINE__inline | |||
281 | void | |||
282 | ss_fixdown(const unsigned char *Td, const int *PA, | |||
283 | int *SA, int i, int size) { | |||
284 | int j, k; | |||
285 | int v; | |||
286 | int c, d, e; | |||
287 | ||||
288 | for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { | |||
289 | d = Td[PA[SA[k = j++]]]; | |||
290 | if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; } | |||
291 | if(d <= c) { break; } | |||
292 | } | |||
293 | SA[i] = v; | |||
294 | } | |||
295 | ||||
296 | /* Simple top-down heapsort. */ | |||
297 | static | |||
298 | void | |||
299 | ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) { | |||
300 | int i, m; | |||
301 | int t; | |||
302 | ||||
303 | m = size; | |||
304 | if((size % 2) == 0) { | |||
305 | m--; | |||
306 | if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2])do { t = (SA[m]); (SA[m]) = (SA[m / 2]); (SA[m / 2]) = t; } while (0); } | |||
307 | } | |||
308 | ||||
309 | for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); } | |||
310 | if((size % 2) == 0) { SWAP(SA[0], SA[m])do { t = (SA[0]); (SA[0]) = (SA[m]); (SA[m]) = t; } while(0); ss_fixdown(Td, PA, SA, 0, m); } | |||
311 | for(i = m - 1; 0 < i; --i) { | |||
312 | t = SA[0], SA[0] = SA[i]; | |||
313 | ss_fixdown(Td, PA, SA, 0, i); | |||
314 | SA[i] = t; | |||
315 | } | |||
316 | } | |||
317 | ||||
318 | ||||
319 | /*---------------------------------------------------------------------------*/ | |||
320 | ||||
321 | /* Returns the median of three elements. */ | |||
322 | static INLINE__inline | |||
323 | int * | |||
324 | ss_median3(const unsigned char *Td, const int *PA, | |||
325 | int *v1, int *v2, int *v3) { | |||
326 | int *t; | |||
327 | if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2)do { t = (v1); (v1) = (v2); (v2) = t; } while(0); } | |||
328 | if(Td[PA[*v2]] > Td[PA[*v3]]) { | |||
329 | if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; } | |||
330 | else { return v3; } | |||
331 | } | |||
332 | return v2; | |||
333 | } | |||
334 | ||||
335 | /* Returns the median of five elements. */ | |||
336 | static INLINE__inline | |||
337 | int * | |||
338 | ss_median5(const unsigned char *Td, const int *PA, | |||
339 | int *v1, int *v2, int *v3, int *v4, int *v5) { | |||
340 | int *t; | |||
341 | if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3)do { t = (v2); (v2) = (v3); (v3) = t; } while(0); } | |||
342 | if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5)do { t = (v4); (v4) = (v5); (v5) = t; } while(0); } | |||
343 | if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4)do { t = (v2); (v2) = (v4); (v4) = t; } while(0); SWAP(v3, v5)do { t = (v3); (v3) = (v5); (v5) = t; } while(0); } | |||
344 | if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3)do { t = (v1); (v1) = (v3); (v3) = t; } while(0); } | |||
345 | if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4)do { t = (v1); (v1) = (v4); (v4) = t; } while(0); SWAP(v3, v5)do { t = (v3); (v3) = (v5); (v5) = t; } while(0); } | |||
346 | if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; } | |||
347 | return v3; | |||
348 | } | |||
349 | ||||
350 | /* Returns the pivot element. */ | |||
351 | static INLINE__inline | |||
352 | int * | |||
353 | ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) { | |||
354 | int *middle; | |||
355 | int t; | |||
356 | ||||
357 | t = last - first; | |||
358 | middle = first + t / 2; | |||
359 | ||||
360 | if(t <= 512) { | |||
361 | if(t <= 32) { | |||
362 | return ss_median3(Td, PA, first, middle, last - 1); | |||
363 | } else { | |||
364 | t >>= 2; | |||
365 | return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1); | |||
366 | } | |||
367 | } | |||
368 | t >>= 3; | |||
369 | first = ss_median3(Td, PA, first, first + t, first + (t << 1)); | |||
370 | middle = ss_median3(Td, PA, middle - t, middle, middle + t); | |||
371 | last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1); | |||
372 | return ss_median3(Td, PA, first, middle, last); | |||
373 | } | |||
374 | ||||
375 | ||||
376 | /*---------------------------------------------------------------------------*/ | |||
377 | ||||
378 | /* Binary partition for substrings. */ | |||
379 | static INLINE__inline | |||
380 | int * | |||
381 | ss_partition(const int *PA, | |||
382 | int *first, int *last, int depth) { | |||
383 | int *a, *b; | |||
384 | int t; | |||
385 | for(a = first - 1, b = last;;) { | |||
386 | for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; } | |||
387 | for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { } | |||
388 | if(b <= a) { break; } | |||
389 | t = ~*b; | |||
390 | *b = *a; | |||
391 | *a = t; | |||
392 | } | |||
393 | if(first < a) { *first = ~*first; } | |||
394 | return a; | |||
395 | } | |||
396 | ||||
397 | /* Multikey introsort for medium size groups. */ | |||
398 | static | |||
399 | void | |||
400 | ss_mintrosort(const unsigned char *T, const int *PA, | |||
401 | int *first, int *last, | |||
402 | int depth) { | |||
403 | #define STACK_SIZE SS_MISORT_STACKSIZE(16) | |||
404 | struct { int *a, *b, c; int d; } stack[STACK_SIZE]; | |||
405 | const unsigned char *Td; | |||
406 | int *a, *b, *c, *d, *e, *f; | |||
407 | int s, t; | |||
408 | int ssize; | |||
409 | int limit; | |||
410 | int v, x = 0; | |||
411 | ||||
412 | for(ssize = 0, limit = ss_ilg(last - first);;) { | |||
413 | ||||
414 | if((last - first) <= SS_INSERTIONSORT_THRESHOLD(8)) { | |||
415 | #if 1 < SS_INSERTIONSORT_THRESHOLD(8) | |||
416 | if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); } | |||
417 | #endif | |||
418 | STACK_POP(first, last, depth, limit)do { ((void) (0)); if(ssize == 0) { return; } (first) = stack [--ssize].a, (last) = stack[ssize].b, (depth) = stack[ssize]. c, (limit) = stack[ssize].d; } while(0); | |||
419 | continue; | |||
420 | } | |||
421 | ||||
422 | Td = T + depth; | |||
423 | if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); } | |||
424 | if(limit < 0) { | |||
425 | for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) { | |||
426 | if((x = Td[PA[*a]]) != v) { | |||
427 | if(1 < (a - first)) { break; } | |||
428 | v = x; | |||
429 | first = a; | |||
430 | } | |||
431 | } | |||
432 | if(Td[PA[*first] - 1] < v) { | |||
433 | first = ss_partition(PA, first, a, depth); | |||
434 | } | |||
435 | if((a - first) <= (last - a)) { | |||
436 | if(1 < (a - first)) { | |||
437 | STACK_PUSH(a, last, depth, -1)do { ((void) (0)); stack[ssize].a = (a), stack[ssize].b = (last ), stack[ssize].c = (depth), stack[ssize++].d = (-1); } while (0); | |||
438 | last = a, depth += 1, limit = ss_ilg(a - first); | |||
439 | } else { | |||
440 | first = a, limit = -1; | |||
441 | } | |||
442 | } else { | |||
443 | if(1 < (last - a)) { | |||
444 | STACK_PUSH(first, a, depth + 1, ss_ilg(a - first))do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (a), stack[ssize].c = (depth + 1), stack[ssize++].d = (ss_ilg (a - first)); } while(0); | |||
445 | first = a, limit = -1; | |||
446 | } else { | |||
447 | last = a, depth += 1, limit = ss_ilg(a - first); | |||
448 | } | |||
449 | } | |||
450 | continue; | |||
451 | } | |||
452 | ||||
453 | /* choose pivot */ | |||
454 | a = ss_pivot(Td, PA, first, last); | |||
455 | v = Td[PA[*a]]; | |||
456 | SWAP(*first, *a)do { t = (*first); (*first) = (*a); (*a) = t; } while(0); | |||
457 | ||||
458 | /* partition */ | |||
459 | for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { } | |||
460 | if(((a = b) < last) && (x < v)) { | |||
461 | for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) { | |||
462 | if(x == v) { SWAP(*b, *a)do { t = (*b); (*b) = (*a); (*a) = t; } while(0); ++a; } | |||
463 | } | |||
464 | } | |||
465 | for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { } | |||
466 | if((b < (d = c)) && (x > v)) { | |||
467 | for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { | |||
468 | if(x == v) { SWAP(*c, *d)do { t = (*c); (*c) = (*d); (*d) = t; } while(0); --d; } | |||
469 | } | |||
470 | } | |||
471 | for(; b < c;) { | |||
472 | SWAP(*b, *c)do { t = (*b); (*b) = (*c); (*c) = t; } while(0); | |||
473 | for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) { | |||
474 | if(x == v) { SWAP(*b, *a)do { t = (*b); (*b) = (*a); (*a) = t; } while(0); ++a; } | |||
475 | } | |||
476 | for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { | |||
477 | if(x == v) { SWAP(*c, *d)do { t = (*c); (*c) = (*d); (*d) = t; } while(0); --d; } | |||
478 | } | |||
479 | } | |||
480 | ||||
481 | if(a <= d) { | |||
482 | c = b - 1; | |||
483 | ||||
484 | if((s = a - first) > (t = b - a)) { s = t; } | |||
485 | for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f)do { t = (*e); (*e) = (*f); (*f) = t; } while(0); } | |||
486 | if((s = d - c) > (t = last - d - 1)) { s = t; } | |||
487 | for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f)do { t = (*e); (*e) = (*f); (*f) = t; } while(0); } | |||
488 | ||||
489 | a = first + (b - a), c = last - (d - c); | |||
490 | b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth); | |||
491 | ||||
492 | if((a - first) <= (last - c)) { | |||
493 | if((last - c) <= (c - b)) { | |||
494 | STACK_PUSH(b, c, depth + 1, ss_ilg(c - b))do { ((void) (0)); stack[ssize].a = (b), stack[ssize].b = (c) , stack[ssize].c = (depth + 1), stack[ssize++].d = (ss_ilg(c - b)); } while(0); | |||
495 | STACK_PUSH(c, last, depth, limit)do { ((void) (0)); stack[ssize].a = (c), stack[ssize].b = (last ), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while (0); | |||
496 | last = a; | |||
497 | } else if((a - first) <= (c - b)) { | |||
498 | STACK_PUSH(c, last, depth, limit)do { ((void) (0)); stack[ssize].a = (c), stack[ssize].b = (last ), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while (0); | |||
499 | STACK_PUSH(b, c, depth + 1, ss_ilg(c - b))do { ((void) (0)); stack[ssize].a = (b), stack[ssize].b = (c) , stack[ssize].c = (depth + 1), stack[ssize++].d = (ss_ilg(c - b)); } while(0); | |||
500 | last = a; | |||
501 | } else { | |||
502 | STACK_PUSH(c, last, depth, limit)do { ((void) (0)); stack[ssize].a = (c), stack[ssize].b = (last ), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while (0); | |||
503 | STACK_PUSH(first, a, depth, limit)do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (a), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while(0); | |||
504 | first = b, last = c, depth += 1, limit = ss_ilg(c - b); | |||
505 | } | |||
506 | } else { | |||
507 | if((a - first) <= (c - b)) { | |||
508 | STACK_PUSH(b, c, depth + 1, ss_ilg(c - b))do { ((void) (0)); stack[ssize].a = (b), stack[ssize].b = (c) , stack[ssize].c = (depth + 1), stack[ssize++].d = (ss_ilg(c - b)); } while(0); | |||
509 | STACK_PUSH(first, a, depth, limit)do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (a), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while(0); | |||
510 | first = c; | |||
511 | } else if((last - c) <= (c - b)) { | |||
512 | STACK_PUSH(first, a, depth, limit)do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (a), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while(0); | |||
513 | STACK_PUSH(b, c, depth + 1, ss_ilg(c - b))do { ((void) (0)); stack[ssize].a = (b), stack[ssize].b = (c) , stack[ssize].c = (depth + 1), stack[ssize++].d = (ss_ilg(c - b)); } while(0); | |||
514 | first = c; | |||
515 | } else { | |||
516 | STACK_PUSH(first, a, depth, limit)do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (a), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while(0); | |||
517 | STACK_PUSH(c, last, depth, limit)do { ((void) (0)); stack[ssize].a = (c), stack[ssize].b = (last ), stack[ssize].c = (depth), stack[ssize++].d = (limit); } while (0); | |||
518 | first = b, last = c, depth += 1, limit = ss_ilg(c - b); | |||
519 | } | |||
520 | } | |||
521 | } else { | |||
522 | limit += 1; | |||
523 | if(Td[PA[*first] - 1] < v) { | |||
524 | first = ss_partition(PA, first, last, depth); | |||
525 | limit = ss_ilg(last - first); | |||
526 | } | |||
527 | depth += 1; | |||
528 | } | |||
529 | } | |||
530 | #undef STACK_SIZE | |||
531 | } | |||
532 | ||||
533 | #endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ | |||
534 | ||||
535 | ||||
536 | /*---------------------------------------------------------------------------*/ | |||
537 | ||||
538 | #if SS_BLOCKSIZE(1024) != 0 | |||
539 | ||||
540 | static INLINE__inline | |||
541 | void | |||
542 | ss_blockswap(int *a, int *b, int n) { | |||
543 | int t; | |||
544 | for(; 0 < n; --n, ++a, ++b) { | |||
545 | t = *a, *a = *b, *b = t; | |||
546 | } | |||
547 | } | |||
548 | ||||
549 | static INLINE__inline | |||
550 | void | |||
551 | ss_rotate(int *first, int *middle, int *last) { | |||
552 | int *a, *b, t; | |||
553 | int l, r; | |||
554 | l = middle - first, r = last - middle; | |||
555 | for(; (0 < l) && (0 < r);) { | |||
556 | if(l == r) { ss_blockswap(first, middle, l); break; } | |||
557 | if(l < r) { | |||
558 | a = last - 1, b = middle - 1; | |||
559 | t = *a; | |||
560 | do { | |||
561 | *a-- = *b, *b-- = *a; | |||
562 | if(b < first) { | |||
563 | *a = t; | |||
564 | last = a; | |||
565 | if((r -= l + 1) <= l) { break; } | |||
566 | a -= 1, b = middle - 1; | |||
567 | t = *a; | |||
568 | } | |||
569 | } while(1); | |||
570 | } else { | |||
571 | a = first, b = middle; | |||
572 | t = *a; | |||
573 | do { | |||
574 | *a++ = *b, *b++ = *a; | |||
575 | if(last <= b) { | |||
576 | *a = t; | |||
577 | first = a + 1; | |||
578 | if((l -= r + 1) <= r) { break; } | |||
579 | a += 1, b = middle; | |||
580 | t = *a; | |||
581 | } | |||
582 | } while(1); | |||
583 | } | |||
584 | } | |||
585 | } | |||
586 | ||||
587 | ||||
588 | /*---------------------------------------------------------------------------*/ | |||
589 | ||||
590 | static | |||
591 | void | |||
592 | ss_inplacemerge(const unsigned char *T, const int *PA, | |||
593 | int *first, int *middle, int *last, | |||
594 | int depth) { | |||
595 | const int *p; | |||
596 | int *a, *b; | |||
597 | int len, half; | |||
598 | int q, r; | |||
599 | int x; | |||
600 | ||||
601 | for(;;) { | |||
602 | if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); } | |||
603 | else { x = 0; p = PA + *(last - 1); } | |||
604 | for(a = first, len = middle - first, half = len >> 1, r = -1; | |||
605 | 0 < len; | |||
606 | len = half, half >>= 1) { | |||
607 | b = a + half; | |||
608 | q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth); | |||
609 | if(q < 0) { | |||
610 | a = b + 1; | |||
611 | half -= (len & 1) ^ 1; | |||
612 | } else { | |||
613 | r = q; | |||
614 | } | |||
615 | } | |||
616 | if(a < middle) { | |||
617 | if(r == 0) { *a = ~*a; } | |||
618 | ss_rotate(a, middle, last); | |||
619 | last -= middle - a; | |||
620 | middle = a; | |||
621 | if(first == middle) { break; } | |||
622 | } | |||
623 | --last; | |||
624 | if(x != 0) { while(*--last < 0) { } } | |||
625 | if(middle == last) { break; } | |||
626 | } | |||
627 | } | |||
628 | ||||
629 | ||||
630 | /*---------------------------------------------------------------------------*/ | |||
631 | ||||
632 | /* Merge-forward with internal buffer. */ | |||
633 | static | |||
634 | void | |||
635 | ss_mergeforward(const unsigned char *T, const int *PA, | |||
636 | int *first, int *middle, int *last, | |||
637 | int *buf, int depth) { | |||
638 | int *a, *b, *c, *bufend; | |||
639 | int t; | |||
640 | int r; | |||
641 | ||||
642 | bufend = buf + (middle - first) - 1; | |||
643 | ss_blockswap(buf, first, middle - first); | |||
644 | ||||
645 | for(t = *(a = first), b = buf, c = middle;;) { | |||
646 | r = ss_compare(T, PA + *b, PA + *c, depth); | |||
647 | if(r < 0) { | |||
648 | do { | |||
649 | *a++ = *b; | |||
650 | if(bufend <= b) { *bufend = t; return; } | |||
651 | *b++ = *a; | |||
652 | } while(*b < 0); | |||
653 | } else if(r > 0) { | |||
654 | do { | |||
655 | *a++ = *c, *c++ = *a; | |||
656 | if(last <= c) { | |||
657 | while(b < bufend) { *a++ = *b, *b++ = *a; } | |||
658 | *a = *b, *b = t; | |||
659 | return; | |||
660 | } | |||
661 | } while(*c < 0); | |||
662 | } else { | |||
663 | *c = ~*c; | |||
664 | do { | |||
665 | *a++ = *b; | |||
666 | if(bufend <= b) { *bufend = t; return; } | |||
667 | *b++ = *a; | |||
668 | } while(*b < 0); | |||
669 | ||||
670 | do { | |||
671 | *a++ = *c, *c++ = *a; | |||
672 | if(last <= c) { | |||
673 | while(b < bufend) { *a++ = *b, *b++ = *a; } | |||
674 | *a = *b, *b = t; | |||
675 | return; | |||
676 | } | |||
677 | } while(*c < 0); | |||
678 | } | |||
679 | } | |||
680 | } | |||
681 | ||||
682 | /* Merge-backward with internal buffer. */ | |||
683 | static | |||
684 | void | |||
685 | ss_mergebackward(const unsigned char *T, const int *PA, | |||
686 | int *first, int *middle, int *last, | |||
687 | int *buf, int depth) { | |||
688 | const int *p1, *p2; | |||
689 | int *a, *b, *c, *bufend; | |||
690 | int t; | |||
691 | int r; | |||
692 | int x; | |||
693 | ||||
694 | bufend = buf + (last - middle) - 1; | |||
695 | ss_blockswap(buf, middle, last - middle); | |||
696 | ||||
697 | x = 0; | |||
698 | if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; } | |||
699 | else { p1 = PA + *bufend; } | |||
700 | if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; } | |||
701 | else { p2 = PA + *(middle - 1); } | |||
702 | for(t = *(a = last - 1), b = bufend, c = middle - 1;;) { | |||
703 | r = ss_compare(T, p1, p2, depth); | |||
704 | if(0 < r) { | |||
705 | if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } | |||
706 | *a-- = *b; | |||
707 | if(b <= buf) { *buf = t; break; } | |||
708 | *b-- = *a; | |||
709 | if(*b < 0) { p1 = PA + ~*b; x |= 1; } | |||
710 | else { p1 = PA + *b; } | |||
711 | } else if(r < 0) { | |||
712 | if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } | |||
713 | *a-- = *c, *c-- = *a; | |||
714 | if(c < first) { | |||
715 | while(buf < b) { *a-- = *b, *b-- = *a; } | |||
716 | *a = *b, *b = t; | |||
717 | break; | |||
718 | } | |||
719 | if(*c < 0) { p2 = PA + ~*c; x |= 2; } | |||
720 | else { p2 = PA + *c; } | |||
721 | } else { | |||
722 | if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } | |||
723 | *a-- = ~*b; | |||
724 | if(b <= buf) { *buf = t; break; } | |||
725 | *b-- = *a; | |||
726 | if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } | |||
727 | *a-- = *c, *c-- = *a; | |||
728 | if(c < first) { | |||
729 | while(buf < b) { *a-- = *b, *b-- = *a; } | |||
730 | *a = *b, *b = t; | |||
731 | break; | |||
732 | } | |||
733 | if(*b < 0) { p1 = PA + ~*b; x |= 1; } | |||
734 | else { p1 = PA + *b; } | |||
735 | if(*c < 0) { p2 = PA + ~*c; x |= 2; } | |||
736 | else { p2 = PA + *c; } | |||
737 | } | |||
738 | } | |||
739 | } | |||
740 | ||||
741 | /* D&C based merge. */ | |||
742 | static | |||
743 | void | |||
744 | ss_swapmerge(const unsigned char *T, const int *PA, | |||
745 | int *first, int *middle, int *last, | |||
746 | int *buf, int bufsize, int depth) { | |||
747 | #define STACK_SIZE SS_SMERGE_STACKSIZE(32) | |||
748 | #define GETIDX(a)((0 <= (a)) ? (a) : (~(a))) ((0 <= (a)) ? (a) : (~(a))) | |||
749 | #define MERGE_CHECK(a, b, c)do { if(((c) & 1) || (((c) & 2) && (ss_compare (T, PA + ((0 <= (*((a) - 1))) ? (*((a) - 1)) : (~(*((a) - 1 )))), PA + *(a), depth) == 0))) { *(a) = ~*(a); } if(((c) & 4) && ((ss_compare(T, PA + ((0 <= (*((b) - 1))) ? (*((b) - 1)) : (~(*((b) - 1)))), PA + *(b), depth) == 0))) { *(b) = ~*(b); } } while(0)\ | |||
750 | do {\ | |||
751 | if(((c) & 1) ||\ | |||
752 | (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1))((0 <= (*((a) - 1))) ? (*((a) - 1)) : (~(*((a) - 1)))), PA + *(a), depth) == 0))) {\ | |||
753 | *(a) = ~*(a);\ | |||
754 | }\ | |||
755 | if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1))((0 <= (*((b) - 1))) ? (*((b) - 1)) : (~(*((b) - 1)))), PA + *(b), depth) == 0))) {\ | |||
756 | *(b) = ~*(b);\ | |||
757 | }\ | |||
758 | } while(0) | |||
759 | struct { int *a, *b, *c; int d; } stack[STACK_SIZE]; | |||
760 | int *l, *r, *lm, *rm; | |||
761 | int m, len, half; | |||
762 | int ssize; | |||
763 | int check, next; | |||
764 | ||||
765 | for(check = 0, ssize = 0;;) { | |||
766 | if((last - middle) <= bufsize) { | |||
767 | if((first < middle) && (middle < last)) { | |||
768 | ss_mergebackward(T, PA, first, middle, last, buf, depth); | |||
769 | } | |||
770 | MERGE_CHECK(first, last, check)do { if(((check) & 1) || (((check) & 2) && (ss_compare (T, PA + ((0 <= (*((first) - 1))) ? (*((first) - 1)) : (~( *((first) - 1)))), PA + *(first), depth) == 0))) { *(first) = ~*(first); } if(((check) & 4) && ((ss_compare(T, PA + ((0 <= (*((last) - 1))) ? (*((last) - 1)) : (~(*((last ) - 1)))), PA + *(last), depth) == 0))) { *(last) = ~*(last); } } while(0); | |||
771 | STACK_POP(first, middle, last, check)do { ((void) (0)); if(ssize == 0) { return; } (first) = stack [--ssize].a, (middle) = stack[ssize].b, (last) = stack[ssize] .c, (check) = stack[ssize].d; } while(0); | |||
772 | continue; | |||
773 | } | |||
774 | ||||
775 | if((middle - first) <= bufsize) { | |||
776 | if(first < middle) { | |||
777 | ss_mergeforward(T, PA, first, middle, last, buf, depth); | |||
778 | } | |||
779 | MERGE_CHECK(first, last, check)do { if(((check) & 1) || (((check) & 2) && (ss_compare (T, PA + ((0 <= (*((first) - 1))) ? (*((first) - 1)) : (~( *((first) - 1)))), PA + *(first), depth) == 0))) { *(first) = ~*(first); } if(((check) & 4) && ((ss_compare(T, PA + ((0 <= (*((last) - 1))) ? (*((last) - 1)) : (~(*((last ) - 1)))), PA + *(last), depth) == 0))) { *(last) = ~*(last); } } while(0); | |||
780 | STACK_POP(first, middle, last, check)do { ((void) (0)); if(ssize == 0) { return; } (first) = stack [--ssize].a, (middle) = stack[ssize].b, (last) = stack[ssize] .c, (check) = stack[ssize].d; } while(0); | |||
781 | continue; | |||
782 | } | |||
783 | ||||
784 | for(m = 0, len = MIN(middle - first, last - middle)(((middle - first) < (last - middle)) ? (middle - first) : (last - middle)), half = len >> 1; | |||
785 | 0 < len; | |||
786 | len = half, half >>= 1) { | |||
787 | if(ss_compare(T, PA + GETIDX(*(middle + m + half))((0 <= (*(middle + m + half))) ? (*(middle + m + half)) : ( ~(*(middle + m + half)))), | |||
788 | PA + GETIDX(*(middle - m - half - 1))((0 <= (*(middle - m - half - 1))) ? (*(middle - m - half - 1)) : (~(*(middle - m - half - 1)))), depth) < 0) { | |||
789 | m += half + 1; | |||
790 | half -= (len & 1) ^ 1; | |||
791 | } | |||
792 | } | |||
793 | ||||
794 | if(0 < m) { | |||
795 | lm = middle - m, rm = middle + m; | |||
796 | ss_blockswap(lm, middle, m); | |||
797 | l = r = middle, next = 0; | |||
798 | if(rm < last) { | |||
799 | if(*rm < 0) { | |||
800 | *rm = ~*rm; | |||
801 | if(first < lm) { for(; *--l < 0;) { } next |= 4; } | |||
802 | next |= 1; | |||
803 | } else if(first < lm) { | |||
804 | for(; *r < 0; ++r) { } | |||
805 | next |= 2; | |||
806 | } | |||
807 | } | |||
808 | ||||
809 | if((l - first) <= (last - r)) { | |||
810 | STACK_PUSH(r, rm, last, (next & 3) | (check & 4))do { ((void) (0)); stack[ssize].a = (r), stack[ssize].b = (rm ), stack[ssize].c = (last), stack[ssize++].d = ((next & 3 ) | (check & 4)); } while(0); | |||
811 | middle = lm, last = l, check = (check & 3) | (next & 4); | |||
812 | } else { | |||
813 | if((next & 2) && (r == middle)) { next ^= 6; } | |||
814 | STACK_PUSH(first, lm, l, (check & 3) | (next & 4))do { ((void) (0)); stack[ssize].a = (first), stack[ssize].b = (lm), stack[ssize].c = (l), stack[ssize++].d = ((check & 3) | (next & 4)); } while(0); | |||
815 | first = r, middle = rm, check = (next & 3) | (check & 4); | |||
816 | } | |||
817 | } else { | |||
818 | if(ss_compare(T, PA + GETIDX(*(middle - 1))((0 <= (*(middle - 1))) ? (*(middle - 1)) : (~(*(middle - 1 )))), PA + *middle, depth) == 0) { | |||
819 | *middle = ~*middle; | |||
820 | } | |||
821 | MERGE_CHECK(first, last, check)do { if(((check) & 1) || (((check) & 2) && (ss_compare (T, PA + ((0 <= (*((first) - 1))) ? (*((first) - 1)) : (~( *((first) - 1)))), PA + *(first), depth) == 0))) { *(first) = ~*(first); } if(((check) & 4) && ((ss_compare(T, PA + ((0 <= (*((last) - 1))) ? (*((last) - 1)) : (~(*((last ) - 1)))), PA + *(last), depth) == 0))) { *(last) = ~*(last); } } while(0); | |||
822 | STACK_POP(first, middle, last, check)do { ((void) (0)); if(ssize == 0) { return; } (first) = stack [--ssize].a, (middle) = stack[ssize].b, (last) = stack[ssize] .c, (check) = stack[ssize].d; } while(0); | |||
823 | } | |||
824 | } | |||
825 | #undef STACK_SIZE | |||
826 | } | |||
827 | ||||
828 | #endif /* SS_BLOCKSIZE != 0 */ | |||
829 | ||||
830 | ||||
831 | /*---------------------------------------------------------------------------*/ | |||
832 | ||||
833 | /* Substring sort */ | |||
834 | static | |||
835 | void | |||
836 | sssort(const unsigned char *T, const int *PA, | |||
837 | int *first, int *last, | |||
838 | int *buf, int bufsize, | |||
839 | int depth, int n, int lastsuffix) { | |||
840 | int *a; | |||
841 | #if SS_BLOCKSIZE(1024) != 0 | |||
842 | int *b, *middle, *curbuf; | |||
843 | int j, k, curbufsize, limit; | |||
844 | #endif | |||
845 | int i; | |||
846 | ||||
847 | if(lastsuffix != 0) { ++first; } | |||
848 | ||||
849 | #if SS_BLOCKSIZE(1024) == 0 | |||
850 | ss_mintrosort(T, PA, first, last, depth); | |||
851 | #else | |||
852 | if((bufsize < SS_BLOCKSIZE(1024)) && | |||
853 | (bufsize < (last - first)) && | |||
854 | (bufsize < (limit = ss_isqrt(last - first)))) { | |||
855 | if(SS_BLOCKSIZE(1024) < limit) { limit = SS_BLOCKSIZE(1024); } | |||
856 | buf = middle = last - limit, bufsize = limit; | |||
857 | } else { | |||
858 | middle = last, limit = 0; | |||
859 | } | |||
860 | for(a = first, i = 0; SS_BLOCKSIZE(1024) < (middle - a); a += SS_BLOCKSIZE(1024), ++i) { | |||
861 | #if SS_INSERTIONSORT_THRESHOLD(8) < SS_BLOCKSIZE(1024) | |||
862 | ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE(1024), depth); | |||
863 | #elif 1 < SS_BLOCKSIZE(1024) | |||
864 | ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE(1024), depth); | |||
865 | #endif | |||
866 | curbufsize = last - (a + SS_BLOCKSIZE(1024)); | |||
867 | curbuf = a + SS_BLOCKSIZE(1024); | |||
868 | if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; } | |||
869 | for(b = a, k = SS_BLOCKSIZE(1024), j = i; j & 1; b -= k, k <<= 1, j >>= 1) { | |||
870 | ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth); | |||
871 | } | |||
872 | } | |||
873 | #if SS_INSERTIONSORT_THRESHOLD(8) < SS_BLOCKSIZE(1024) | |||
874 | ss_mintrosort(T, PA, a, middle, depth); | |||
875 | #elif 1 < SS_BLOCKSIZE(1024) | |||
876 | ss_insertionsort(T, PA, a, middle, depth); | |||
877 | #endif | |||
878 | for(k = SS_BLOCKSIZE(1024); i != 0; k <<= 1, i >>= 1) { | |||
879 | if(i & 1) { | |||
880 | ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth); | |||
881 | a -= k; | |||
882 | } | |||
883 | } | |||
884 | if(limit != 0) { | |||
885 | #if SS_INSERTIONSORT_THRESHOLD(8) < SS_BLOCKSIZE(1024) | |||
886 | ss_mintrosort(T, PA, middle, last, depth); | |||
887 | #elif 1 < SS_BLOCKSIZE(1024) | |||
888 | ss_insertionsort(T, PA, middle, last, depth); | |||
889 | #endif | |||
890 | ss_inplacemerge(T, PA, first, middle, last, depth); | |||
891 | } | |||
892 | #endif | |||
893 | ||||
894 | if(lastsuffix != 0) { | |||
895 | /* Insert last type B* suffix. */ | |||
896 | int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2; | |||
897 | for(a = first, i = *(first - 1); | |||
898 | (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth))); | |||
899 | ++a) { | |||
900 | *(a - 1) = *a; | |||
901 | } | |||
902 | *(a - 1) = i; | |||
903 | } | |||
904 | } | |||
905 | ||||
906 | ||||
907 | /*---------------------------------------------------------------------------*/ | |||
908 | ||||
909 | static INLINE__inline | |||
910 | int | |||
911 | tr_ilg(int n) { | |||
912 | return (n & 0xffff0000) ? | |||
913 | ((n & 0xff000000) ? | |||
914 | 24 + lg_table[(n >> 24) & 0xff] : | |||
915 | 16 + lg_table[(n >> 16) & 0xff]) : | |||
916 | ((n & 0x0000ff00) ? | |||
917 | 8 + lg_table[(n >> 8) & 0xff] : | |||
918 | 0 + lg_table[(n >> 0) & 0xff]); | |||
919 | } | |||
920 | ||||
921 | ||||
922 | /*---------------------------------------------------------------------------*/ | |||
923 | ||||
924 | /* Simple insertionsort for small size groups. */ | |||
925 | static | |||
926 | void | |||
927 | tr_insertionsort(const int *ISAd, int *first, int *last) { | |||
928 | int *a, *b; | |||
929 | int t, r; | |||
930 | ||||
931 | for(a = first + 1; a < last; ++a) { | |||
932 | for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) { | |||
933 | do { *(b + 1) = *b; } while((first <= --b) && (*b < 0)); | |||
934 | if(b < first) { break; } | |||
935 | } | |||
936 | if(r == 0) { *b = ~*b; } | |||
937 | *(b + 1) = t; | |||
938 | } | |||
939 | } | |||
940 | ||||
941 | ||||
942 | /*---------------------------------------------------------------------------*/ | |||
943 | ||||
944 | static INLINE__inline | |||
945 | void | |||
946 | tr_fixdown(const int *ISAd, int *SA, int i, int size) { | |||
947 | int j, k; | |||
948 | int v; | |||
949 | int c, d, e; | |||
950 | ||||
951 | for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { | |||
952 | d = ISAd[SA[k = j++]]; | |||
953 | if(d < (e = ISAd[SA[j]])) { k = j; d = e; } | |||
954 | if(d <= c) { break; } | |||
955 | } | |||
956 | SA[i] = v; | |||
957 | } | |||
958 | ||||
959 | /* Simple top-down heapsort. */ | |||
960 | static | |||
961 | void | |||
962 | tr_heapsort(const int *ISAd, int *SA, int size) { | |||
963 | int i, m; | |||
964 | int t; | |||
965 | ||||
966 | m = size; | |||
967 | if((size % 2) == 0) { | |||
968 | m--; | |||
969 | if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2])do { t = (SA[m]); (SA[m]) = (SA[m / 2]); (SA[m / 2]) = t; } while (0); } | |||
970 | } | |||
971 | ||||
972 | for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); } | |||
973 | if((size % 2) == 0) { SWAP(SA[0], SA[m])do { t = (SA[0]); (SA[0]) = (SA[m]); (SA[m]) = t; } while(0); tr_fixdown(ISAd, SA, 0, m); } | |||
974 | for(i = m - 1; 0 < i; --i) { | |||
975 | t = SA[0], SA[0] = SA[i]; | |||
976 | tr_fixdown(ISAd, SA, 0, i); | |||
977 | SA[i] = t; | |||
978 | } | |||
979 | } | |||
980 | ||||
981 | ||||
982 | /*---------------------------------------------------------------------------*/ | |||
983 | ||||
984 | /* Returns the median of three elements. */ | |||
985 | static INLINE__inline | |||
986 | int * | |||
987 | tr_median3(const int *ISAd, int *v1, int *v2, int *v3) { | |||
988 | int *t; | |||
989 | if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2)do { t = (v1); (v1) = (v2); (v2) = t; } while(0); } | |||
990 | if(ISAd[*v2] > ISAd[*v3]) { | |||
991 | if(ISAd[*v1] > ISAd[*v3]) { return v1; } | |||
992 | else { return v3; } | |||
993 | } | |||
994 | return v2; | |||
995 | } | |||
996 | ||||
997 | /* Returns the median of five elements. */ | |||
998 | static INLINE__inline | |||
999 | int * | |||
1000 | tr_median5(const int *ISAd, | |||
1001 | int *v1, int *v2, int *v3, int *v4, int *v5) { | |||
1002 | int *t; | |||
1003 | if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3)do { t = (v2); (v2) = (v3); (v3) = t; } while(0); } | |||
1004 | if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5)do { t = (v4); (v4) = (v5); (v5) = t; } while(0); } | |||
1005 | if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4)do { t = (v2); (v2) = (v4); (v4) = t; } while(0); SWAP(v3, v5)do { t = (v3); (v3) = (v5); (v5) = t; } while(0); } | |||
1006 | if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3)do { t = (v1); (v1) = (v3); (v3) = t; } while(0); } | |||
1007 | if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4)do { t = (v1); (v1) = (v4); (v4) = t; } while(0); SWAP(v3, v5)do { t = (v3); (v3) = (v5); (v5) = t; } while(0); } | |||
1008 | if(ISAd[*v3] > ISAd[*v4]) { return v4; } | |||
1009 | return v3; | |||
1010 | } | |||
1011 | ||||
1012 | /* Returns the pivot element. */ | |||
1013 | static INLINE__inline | |||
1014 | int * | |||
1015 | tr_pivot(const int *ISAd, int *first, int *last) { | |||
1016 | int *middle; | |||
1017 | int t; | |||
1018 | ||||
1019 | t = last - first; | |||
1020 | middle = first + t / 2; | |||
1021 | ||||
1022 | if(t <= 512) { | |||
1023 | if(t <= 32) { | |||
1024 | return tr_median3(ISAd, first, middle, last - 1); | |||
1025 | } else { | |||
1026 | t >>= 2; | |||
1027 | return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1); | |||
1028 | } | |||
1029 | } | |||
1030 | t >>= 3; | |||
1031 | first = tr_median3(ISAd, first, first + t, first + (t << 1)); | |||
1032 | middle = tr_median3(ISAd, middle - t, middle, middle + t); | |||
1033 | last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1); | |||
1034 | return tr_median3(ISAd, first, middle, last); | |||
1035 | } | |||
1036 | ||||
1037 | ||||
1038 | /*---------------------------------------------------------------------------*/ | |||
1039 | ||||
1040 | typedef struct _trbudget_t trbudget_t; | |||
1041 | struct _trbudget_t { | |||
1042 | int chance; | |||
1043 | int remain; | |||
1044 | int incval; | |||
1045 | int count; | |||
1046 | }; | |||
1047 | ||||
1048 | static INLINE__inline | |||
1049 | void | |||
1050 | trbudget_init(trbudget_t *budget, int chance, int incval) { | |||
1051 | budget->chance = chance; | |||
1052 | budget->remain = budget->incval = incval; | |||
1053 | } | |||
1054 | ||||
1055 | static INLINE__inline | |||
1056 | int | |||
1057 | trbudget_check(trbudget_t *budget, int size) { | |||
1058 | if(size <= budget->remain) { budget->remain -= size; return 1; } | |||
1059 | if(budget->chance == 0) { budget->count += size; return 0; } | |||
1060 | budget->remain += budget->incval - size; | |||
1061 | budget->chance -= 1; | |||
1062 | return 1; | |||
1063 | } | |||
1064 | ||||
1065 | ||||
1066 | /*---------------------------------------------------------------------------*/ | |||
1067 | ||||
1068 | static INLINE__inline | |||
1069 | void | |||
1070 | tr_partition(const int *ISAd, | |||
1071 | int *first, int *middle, int *last, | |||
1072 | int **pa, int **pb, int v) { | |||
1073 | int *a, *b, *c, *d, *e, *f; | |||
1074 | int t, s; | |||
1075 | int x = 0; | |||
1076 | ||||
1077 | for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { } | |||
1078 | if(((a = b) < last) && (x < v)) { | |||
1079 | for(; (++b < last) && ((x = ISAd[*b]) <= v);) { | |||
1080 | if(x == v) { SWAP(*b, *a)do { t = (*b); (*b) = (*a); (*a) = t; } while(0); ++a; } | |||
1081 | } | |||
1082 | } | |||
1083 | for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { } | |||
1084 | if((b < (d = c)) && (x > v)) { | |||
1085 | for(; (b < --c) && ((x = ISAd[*c]) >= v);) { | |||
1086 | if(x == v) { SWAP(*c, *d)do { t = (*c); (*c) = (*d); (*d) = t; } while(0); --d; } | |||
1087 | } | |||
1088 | } | |||
1089 | for(; b < c;) { | |||
1090 | SWAP(*b, *c)do { t = (*b); (*b) = (*c); (*c) = t; } while(0); | |||
1091 | for(; (++b < c) && ((x = ISAd[*b]) <= v);) { | |||
1092 | if(x == v) { SWAP(*b, *a)do { t = (*b); (*b) = (*a); (*a) = t; } while(0); ++a; } | |||
1093 | } | |||
1094 | for(; (b < --c) && ((x = ISAd[*c]) >= v);) { | |||
1095 | if(x == v) { SWAP(*c, *d)do { t = (*c); (*c) = (*d); (*d) = t; } while(0); --d; } | |||
1096 | } | |||
1097 | } | |||
1098 | ||||
1099 | if(a <= d) { | |||
1100 | c = b - 1; | |||
1101 | if((s = a - first) > (t = b - a)) { s = t; } | |||
1102 | for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f)do { t = (*e); (*e) = (*f); (*f) = t; } while(0); } | |||
1103 | if((s = d - c) > (t = last - d - 1)) { s = t; } | |||
1104 | for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f)do { t = (*e); (*e) = (*f); (*f) = t; } while(0); } | |||
1105 | first += (b - a), last -= (d - c); | |||
1106 | } | |||
1107 | *pa = first, *pb = last; | |||
1108 | } | |||
1109 | ||||
1110 | static | |||
1111 | void | |||
1112 | tr_copy(int *ISA, const int *SA, | |||
1113 | int *first, int *a, int *b, int *last, | |||
1114 | int depth) { | |||
1115 | /* sort suffixes of middle partition | |||
1116 | by using sorted order of suffixes of left and right partition. */ | |||
1117 | int *c, *d, *e; | |||
1118 | int s, v; | |||
1119 | ||||
1120 | v = b - SA - 1; | |||
1121 | for(c = first, d = a - 1; c <= d; ++c) { | |||
1122 | if((0 <= (s = *c - depth)) && (ISA[s] == v)) { | |||
1123 | *++d = s; | |||
1124 | ISA[s] = d - SA; | |||
1125 | } | |||
1126 | } | |||
1127 | for(c = last - 1, e = d + 1, d = b; e < d; --c) { | |||
1128 | if((0 <= (s = *c - depth)) && (ISA[s] == v)) { | |||
1129 | *--d = s; | |||
1130 | ISA[s] = d - SA; | |||
1131 | } | |||
1132 | } | |||
1133 | } | |||
1134 | ||||
1135 | static | |||
1136 | void | |||
1137 | tr_partialcopy(int *ISA, const int *SA, | |||
1138 | int *first, int *a, int *b, int *last, | |||
1139 | int depth) { | |||
1140 | int *c, *d, *e; | |||
1141 | int s, v; | |||
1142 | int rank, lastrank, newrank = -1; | |||
1143 | ||||
1144 | v = b - SA - 1; | |||
1145 | lastrank = -1; | |||
1146 | for(c = first, d = a - 1; c <= d; ++c) { | |||
1147 | if((0 <= (s = *c - depth)) && (ISA[s] == v)) { | |||
1148 | *++d = s; | |||
1149 | rank = ISA[s + depth]; | |||
1150 | if(lastrank != rank) { lastrank = rank; newrank = d - SA; } | |||
1151 | ISA[s] = newrank; | |||
1152 | } | |||
1153 | } | |||
1154 | ||||
1155 | lastrank = -1; | |||
1156 | for(e = d; first <= e; --e) { | |||
1157 | rank = ISA[*e]; | |||
1158 | if(lastrank != rank) { lastrank = rank; newrank = e - SA; } | |||
1159 | if(newrank != rank) { ISA[*e] = newrank; } | |||
1160 | } | |||
1161 | ||||
1162 | lastrank = -1; | |||
1163 | for(c = last - 1, e = d + 1, d = b; e < d; --c) { | |||
1164 | if((0 <= (s = *c - depth)) && (ISA[s] == v)) { | |||
1165 | *--d = s; | |||
1166 | rank = ISA[s + depth]; | |||
1167 | if(lastrank != rank) { lastrank = rank; newrank = d - SA; } | |||
1168 | ISA[s] = newrank; | |||
1169 | } | |||
1170 | } | |||
1171 | } | |||
1172 | ||||
1173 | static | |||
1174 | void | |||
1175 | tr_introsort(int *ISA, const int *ISAd, | |||
1176 | int *SA, int *first, int *last, | |||
1177 | trbudget_t *budget) { | |||
1178 | #define STACK_SIZE TR_STACKSIZE(64) | |||
1179 | struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE]; | |||
1180 | int *a, *b, *c; | |||
1181 | int t; | |||
1182 | int v, x = 0; | |||
1183 | int incr = ISAd - ISA; | |||
1184 | int limit, next; | |||
1185 | int ssize, trlink = -1; | |||
1186 | ||||
1187 | for(ssize = 0, limit = tr_ilg(last - first);;) { | |||
1188 | ||||
1189 | if(limit < 0) { | |||
1190 | if(limit == -1) { | |||
1191 | /* tandem repeat partition */ | |||
1192 | tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1); | |||
1193 | ||||
1194 | /* update ranks */ | |||
1195 | if(a < last) { | |||
1196 | for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } | |||
1197 | } | |||
1198 | if(b < last) { | |||
1199 | for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } | |||
1200 | } | |||
1201 | ||||
1202 | /* push */ | |||
1203 | if(1 < (b - a)) { | |||
1204 | STACK_PUSH5(NULL, a, b, 0, 0)do { ((void) (0)); stack[ssize].a = (((void*)0)), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (0), stack[ ssize++].e = (0); } while(0); | |||
1205 | STACK_PUSH5(ISAd - incr, first, last, -2, trlink)do { ((void) (0)); stack[ssize].a = (ISAd - incr), stack[ssize ].b = (first), stack[ssize].c = (last), stack[ssize].d = (-2) , stack[ssize++].e = (trlink); } while(0); | |||
1206 | trlink = ssize - 2; | |||
1207 | } | |||
1208 | if((a - first) <= (last - b)) { | |||
1209 | if(1 < (a - first)) { | |||
1210 | STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (tr_ilg(last - b )), stack[ssize++].e = (trlink); } while(0); | |||
1211 | last = a, limit = tr_ilg(a - first); | |||
1212 | } else if(1 < (last - b)) { | |||
1213 | first = b, limit = tr_ilg(last - b); | |||
1214 | } else { | |||
1215 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1216 | } | |||
1217 | } else { | |||
1218 | if(1 < (last - b)) { | |||
1219 | STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (tr_ilg(a - first )), stack[ssize++].e = (trlink); } while(0); | |||
1220 | first = b, limit = tr_ilg(last - b); | |||
1221 | } else if(1 < (a - first)) { | |||
1222 | last = a, limit = tr_ilg(a - first); | |||
1223 | } else { | |||
1224 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1225 | } | |||
1226 | } | |||
1227 | } else if(limit == -2) { | |||
1228 | /* tandem repeat copy */ | |||
1229 | a = stack[--ssize].b, b = stack[ssize].c; | |||
1230 | if(stack[ssize].d == 0) { | |||
1231 | tr_copy(ISA, SA, first, a, b, last, ISAd - ISA); | |||
1232 | } else { | |||
1233 | if(0 <= trlink) { stack[trlink].d = -1; } | |||
1234 | tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA); | |||
1235 | } | |||
1236 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1237 | } else { | |||
1238 | /* sorted partition */ | |||
1239 | if(0 <= *first) { | |||
1240 | a = first; | |||
1241 | do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a)); | |||
1242 | first = a; | |||
1243 | } | |||
1244 | if(first < last) { | |||
1245 | a = first; do { *a = ~*a; } while(*++a < 0); | |||
1246 | next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1; | |||
1247 | if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } } | |||
1248 | ||||
1249 | /* push */ | |||
1250 | if(trbudget_check(budget, a - first)) { | |||
1251 | if((a - first) <= (last - a)) { | |||
1252 | STACK_PUSH5(ISAd, a, last, -3, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( a), stack[ssize].c = (last), stack[ssize].d = (-3), stack[ssize ++].e = (trlink); } while(0); | |||
1253 | ISAd += incr, last = a, limit = next; | |||
1254 | } else { | |||
1255 | if(1 < (last - a)) { | |||
1256 | STACK_PUSH5(ISAd + incr, first, a, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (first), stack[ssize].c = (a), stack[ssize].d = (next), stack[ssize++].e = (trlink); } while(0); | |||
1257 | first = a, limit = -3; | |||
1258 | } else { | |||
1259 | ISAd += incr, last = a, limit = next; | |||
1260 | } | |||
1261 | } | |||
1262 | } else { | |||
1263 | if(0 <= trlink) { stack[trlink].d = -1; } | |||
1264 | if(1 < (last - a)) { | |||
1265 | first = a, limit = -3; | |||
1266 | } else { | |||
1267 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1268 | } | |||
1269 | } | |||
1270 | } else { | |||
1271 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1272 | } | |||
1273 | } | |||
1274 | continue; | |||
1275 | } | |||
1276 | ||||
1277 | if((last - first) <= TR_INSERTIONSORT_THRESHOLD(8)) { | |||
1278 | tr_insertionsort(ISAd, first, last); | |||
1279 | limit = -3; | |||
1280 | continue; | |||
1281 | } | |||
1282 | ||||
1283 | if(limit-- == 0) { | |||
1284 | tr_heapsort(ISAd, first, last - first); | |||
1285 | for(a = last - 1; first < a; a = b) { | |||
1286 | for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; } | |||
1287 | } | |||
1288 | limit = -3; | |||
1289 | continue; | |||
1290 | } | |||
1291 | ||||
1292 | /* choose pivot */ | |||
1293 | a = tr_pivot(ISAd, first, last); | |||
1294 | SWAP(*first, *a)do { t = (*first); (*first) = (*a); (*a) = t; } while(0); | |||
1295 | v = ISAd[*first]; | |||
1296 | ||||
1297 | /* partition */ | |||
1298 | tr_partition(ISAd, first, first + 1, last, &a, &b, v); | |||
1299 | if((last - first) != (b - a)) { | |||
1300 | next = (ISA[*a] != v) ? tr_ilg(b - a) : -1; | |||
1301 | ||||
1302 | /* update ranks */ | |||
1303 | for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } | |||
1304 | if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } } | |||
1305 | ||||
1306 | /* push */ | |||
1307 | if((1 < (b - a)) && (trbudget_check(budget, b - a))) { | |||
1308 | if((a - first) <= (last - b)) { | |||
1309 | if((last - b) <= (b - a)) { | |||
1310 | if(1 < (a - first)) { | |||
1311 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1312 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1313 | last = a; | |||
1314 | } else if(1 < (last - b)) { | |||
1315 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1316 | first = b; | |||
1317 | } else { | |||
1318 | ISAd += incr, first = a, last = b, limit = next; | |||
1319 | } | |||
1320 | } else if((a - first) <= (b - a)) { | |||
1321 | if(1 < (a - first)) { | |||
1322 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1323 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1324 | last = a; | |||
1325 | } else { | |||
1326 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1327 | ISAd += incr, first = a, last = b, limit = next; | |||
1328 | } | |||
1329 | } else { | |||
1330 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1331 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1332 | ISAd += incr, first = a, last = b, limit = next; | |||
1333 | } | |||
1334 | } else { | |||
1335 | if((a - first) <= (b - a)) { | |||
1336 | if(1 < (last - b)) { | |||
1337 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1338 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1339 | first = b; | |||
1340 | } else if(1 < (a - first)) { | |||
1341 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1342 | last = a; | |||
1343 | } else { | |||
1344 | ISAd += incr, first = a, last = b, limit = next; | |||
1345 | } | |||
1346 | } else if((last - b) <= (b - a)) { | |||
1347 | if(1 < (last - b)) { | |||
1348 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1349 | STACK_PUSH5(ISAd + incr, a, b, next, trlink)do { ((void) (0)); stack[ssize].a = (ISAd + incr), stack[ssize ].b = (a), stack[ssize].c = (b), stack[ssize].d = (next), stack [ssize++].e = (trlink); } while(0); | |||
1350 | first = b; | |||
1351 | } else { | |||
1352 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1353 | ISAd += incr, first = a, last = b, limit = next; | |||
1354 | } | |||
1355 | } else { | |||
1356 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1357 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1358 | ISAd += incr, first = a, last = b, limit = next; | |||
1359 | } | |||
1360 | } | |||
1361 | } else { | |||
1362 | if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; } | |||
1363 | if((a - first) <= (last - b)) { | |||
1364 | if(1 < (a - first)) { | |||
1365 | STACK_PUSH5(ISAd, b, last, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( b), stack[ssize].c = (last), stack[ssize].d = (limit), stack[ ssize++].e = (trlink); } while(0); | |||
1366 | last = a; | |||
1367 | } else if(1 < (last - b)) { | |||
1368 | first = b; | |||
1369 | } else { | |||
1370 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1371 | } | |||
1372 | } else { | |||
1373 | if(1 < (last - b)) { | |||
1374 | STACK_PUSH5(ISAd, first, a, limit, trlink)do { ((void) (0)); stack[ssize].a = (ISAd), stack[ssize].b = ( first), stack[ssize].c = (a), stack[ssize].d = (limit), stack [ssize++].e = (trlink); } while(0); | |||
1375 | first = b; | |||
1376 | } else if(1 < (a - first)) { | |||
1377 | last = a; | |||
1378 | } else { | |||
1379 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1380 | } | |||
1381 | } | |||
1382 | } | |||
1383 | } else { | |||
1384 | if(trbudget_check(budget, last - first)) { | |||
1385 | limit = tr_ilg(last - first), ISAd += incr; | |||
1386 | } else { | |||
1387 | if(0 <= trlink) { stack[trlink].d = -1; } | |||
1388 | STACK_POP5(ISAd, first, last, limit, trlink)do { ((void) (0)); if(ssize == 0) { return; } (ISAd) = stack[ --ssize].a, (first) = stack[ssize].b, (last) = stack[ssize].c , (limit) = stack[ssize].d, (trlink) = stack[ssize].e; } while (0); | |||
1389 | } | |||
1390 | } | |||
1391 | } | |||
1392 | #undef STACK_SIZE | |||
1393 | } | |||
1394 | ||||
1395 | ||||
1396 | ||||
1397 | /*---------------------------------------------------------------------------*/ | |||
1398 | ||||
1399 | /* Tandem repeat sort */ | |||
1400 | static | |||
1401 | void | |||
1402 | trsort(int *ISA, int *SA, int n, int depth) { | |||
1403 | int *ISAd; | |||
1404 | int *first, *last; | |||
1405 | trbudget_t budget; | |||
1406 | int t, skip, unsorted; | |||
1407 | ||||
1408 | trbudget_init(&budget, tr_ilg(n) * 2 / 3, n); | |||
1409 | /* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */ | |||
1410 | for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) { | |||
1411 | first = SA; | |||
1412 | skip = 0; | |||
1413 | unsorted = 0; | |||
1414 | do { | |||
1415 | if((t = *first) < 0) { first -= t; skip += t; } | |||
1416 | else { | |||
1417 | if(skip != 0) { *(first + skip) = skip; skip = 0; } | |||
1418 | last = SA + ISA[t] + 1; | |||
1419 | if(1 < (last - first)) { | |||
1420 | budget.count = 0; | |||
1421 | tr_introsort(ISA, ISAd, SA, first, last, &budget); | |||
1422 | if(budget.count != 0) { unsorted += budget.count; } | |||
1423 | else { skip = first - last; } | |||
1424 | } else if((last - first) == 1) { | |||
1425 | skip = -1; | |||
1426 | } | |||
1427 | first = last; | |||
1428 | } | |||
1429 | } while(first < (SA + n)); | |||
1430 | if(skip != 0) { *(first + skip) = skip; } | |||
1431 | if(unsorted == 0) { break; } | |||
1432 | } | |||
1433 | } | |||
1434 | ||||
1435 | ||||
1436 | /*---------------------------------------------------------------------------*/ | |||
1437 | ||||
1438 | /* Sorts suffixes of type B*. */ | |||
1439 | static | |||
1440 | int | |||
1441 | sort_typeBstar(const unsigned char *T, int *SA, | |||
1442 | int *bucket_A, int *bucket_B, | |||
1443 | int n, int openMP) { | |||
1444 | int *PAb, *ISAb, *buf; | |||
1445 | #ifdef LIBBSC_OPENMP | |||
1446 | int *curbuf; | |||
1447 | int l; | |||
1448 | #endif | |||
1449 | int i, j, k, t, m, bufsize; | |||
1450 | int c0, c1; | |||
1451 | #ifdef LIBBSC_OPENMP | |||
1452 | int d0, d1; | |||
1453 | #endif | |||
1454 | (void)openMP; | |||
1455 | ||||
1456 | /* Initialize bucket arrays. */ | |||
1457 | for(i = 0; i < BUCKET_A_SIZE((256)); ++i) { bucket_A[i] = 0; } | |||
1458 | for(i = 0; i < BUCKET_B_SIZE((256) * (256)); ++i) { bucket_B[i] = 0; } | |||
1459 | ||||
1460 | /* Count the number of occurrences of the first one or two characters of each | |||
1461 | type A, B and B* suffix. Moreover, store the beginning position of all | |||
1462 | type B* suffixes into the array SA. */ | |||
1463 | for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) { | |||
1464 | /* type A suffix. */ | |||
1465 | do { ++BUCKET_A(c1 = c0)bucket_A[(c1 = c0)]; } while((0 <= --i) && ((c0 = T[i]) >= c1)); | |||
1466 | if(0 <= i) { | |||
1467 | /* type B* suffix. */ | |||
1468 | ++BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]); | |||
1469 | SA[--m] = i; | |||
1470 | /* type B suffix. */ | |||
1471 | for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { | |||
1472 | ++BUCKET_B(c0, c1)(bucket_B[((c1) << 8) | (c0)]); | |||
1473 | } | |||
1474 | } | |||
1475 | } | |||
1476 | m = n - m; | |||
1477 | /* | |||
1478 | note: | |||
1479 | A type B* suffix is lexicographically smaller than a type B suffix that | |||
1480 | begins with the same first two characters. | |||
1481 | */ | |||
1482 | ||||
1483 | /* Calculate the index of start/end point of each bucket. */ | |||
1484 | for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE(256); ++c0) { | |||
1485 | t = i + BUCKET_A(c0)bucket_A[(c0)]; | |||
1486 | BUCKET_A(c0)bucket_A[(c0)] = i + j; /* start point */ | |||
1487 | i = t + BUCKET_B(c0, c0)(bucket_B[((c0) << 8) | (c0)]); | |||
1488 | for(c1 = c0 + 1; c1 < ALPHABET_SIZE(256); ++c1) { | |||
1489 | j += BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]); | |||
1490 | BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]) = j; /* end point */ | |||
1491 | i += BUCKET_B(c0, c1)(bucket_B[((c1) << 8) | (c0)]); | |||
1492 | } | |||
1493 | } | |||
1494 | ||||
1495 | if(0 < m) { | |||
1496 | /* Sort the type B* suffixes by their first two characters. */ | |||
1497 | PAb = SA + n - m; ISAb = SA + m; | |||
1498 | for(i = m - 2; 0 <= i; --i) { | |||
1499 | t = PAb[i], c0 = T[t], c1 = T[t + 1]; | |||
1500 | SA[--BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)])] = i; | |||
1501 | } | |||
1502 | t = PAb[m - 1], c0 = T[t], c1 = T[t + 1]; | |||
1503 | SA[--BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)])] = m - 1; | |||
1504 | ||||
1505 | /* Sort the type B* substrings using sssort. */ | |||
1506 | #ifdef LIBBSC_OPENMP | |||
1507 | if (openMP) | |||
1508 | { | |||
1509 | buf = SA + m; | |||
1510 | c0 = ALPHABET_SIZE(256) - 2, c1 = ALPHABET_SIZE(256) - 1, j = m; | |||
1511 | #pragma omp parallel default(shared) private(bufsize, curbuf, k, l, d0, d1) | |||
1512 | { | |||
1513 | bufsize = (n - (2 * m)) / omp_get_num_threads(); | |||
1514 | curbuf = buf + omp_get_thread_num() * bufsize; | |||
1515 | k = 0; | |||
1516 | for(;;) { | |||
1517 | #pragma omp critical(sssort_lock) | |||
1518 | { | |||
1519 | if(0 < (l = j)) { | |||
1520 | d0 = c0, d1 = c1; | |||
1521 | do { | |||
1522 | k = BUCKET_BSTAR(d0, d1)(bucket_B[((d0) << 8) | (d1)]); | |||
1523 | if(--d1 <= d0) { | |||
1524 | d1 = ALPHABET_SIZE(256) - 1; | |||
1525 | if(--d0 < 0) { break; } | |||
1526 | } | |||
1527 | } while(((l - k) <= 1) && (0 < (l = k))); | |||
1528 | c0 = d0, c1 = d1, j = k; | |||
1529 | } | |||
1530 | } | |||
1531 | if(l == 0) { break; } | |||
1532 | sssort(T, PAb, SA + k, SA + l, | |||
1533 | curbuf, bufsize, 2, n, *(SA + k) == (m - 1)); | |||
1534 | } | |||
1535 | } | |||
1536 | } | |||
1537 | else | |||
1538 | { | |||
1539 | buf = SA + m, bufsize = n - (2 * m); | |||
1540 | for(c0 = ALPHABET_SIZE(256) - 2, j = m; 0 < j; --c0) { | |||
1541 | for(c1 = ALPHABET_SIZE(256) - 1; c0 < c1; j = i, --c1) { | |||
1542 | i = BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]); | |||
1543 | if(1 < (j - i)) { | |||
1544 | sssort(T, PAb, SA + i, SA + j, | |||
1545 | buf, bufsize, 2, n, *(SA + i) == (m - 1)); | |||
1546 | } | |||
1547 | } | |||
1548 | } | |||
1549 | } | |||
1550 | #else | |||
1551 | buf = SA + m, bufsize = n - (2 * m); | |||
1552 | for(c0 = ALPHABET_SIZE(256) - 2, j = m; 0 < j; --c0) { | |||
1553 | for(c1 = ALPHABET_SIZE(256) - 1; c0 < c1; j = i, --c1) { | |||
1554 | i = BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]); | |||
1555 | if(1 < (j - i)) { | |||
1556 | sssort(T, PAb, SA + i, SA + j, | |||
1557 | buf, bufsize, 2, n, *(SA + i) == (m - 1)); | |||
1558 | } | |||
1559 | } | |||
1560 | } | |||
1561 | #endif | |||
1562 | ||||
1563 | /* Compute ranks of type B* substrings. */ | |||
1564 | for(i = m - 1; 0 <= i; --i) { | |||
1565 | if(0 <= SA[i]) { | |||
1566 | j = i; | |||
1567 | do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i])); | |||
1568 | SA[i + 1] = i - j; | |||
1569 | if(i <= 0) { break; } | |||
1570 | } | |||
1571 | j = i; | |||
1572 | do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0); | |||
1573 | ISAb[SA[i]] = j; | |||
1574 | } | |||
1575 | ||||
1576 | /* Construct the inverse suffix array of type B* suffixes using trsort. */ | |||
1577 | trsort(ISAb, SA, m, 1); | |||
1578 | ||||
1579 | /* Set the sorted order of tyoe B* suffixes. */ | |||
1580 | for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) { | |||
1581 | for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { } | |||
1582 | if(0 <= i) { | |||
1583 | t = i; | |||
1584 | for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { } | |||
1585 | SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t; | |||
1586 | } | |||
1587 | } | |||
1588 | ||||
1589 | /* Calculate the index of start/end point of each bucket. */ | |||
1590 | BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1)(bucket_B[(((256) - 1) << 8) | ((256) - 1)]) = n; /* end point */ | |||
1591 | for(c0 = ALPHABET_SIZE(256) - 2, k = m - 1; 0 <= c0; --c0) { | |||
1592 | i = BUCKET_A(c0 + 1)bucket_A[(c0 + 1)] - 1; | |||
1593 | for(c1 = ALPHABET_SIZE(256) - 1; c0 < c1; --c1) { | |||
1594 | t = i - BUCKET_B(c0, c1)(bucket_B[((c1) << 8) | (c0)]); | |||
1595 | BUCKET_B(c0, c1)(bucket_B[((c1) << 8) | (c0)]) = i; /* end point */ | |||
1596 | ||||
1597 | /* Move all type B* suffixes to the correct position. */ | |||
1598 | for(i = t, j = BUCKET_BSTAR(c0, c1)(bucket_B[((c0) << 8) | (c1)]); | |||
1599 | j <= k; | |||
1600 | --i, --k) { SA[i] = SA[k]; } | |||
1601 | } | |||
1602 | BUCKET_BSTAR(c0, c0 + 1)(bucket_B[((c0) << 8) | (c0 + 1)]) = i - BUCKET_B(c0, c0)(bucket_B[((c0) << 8) | (c0)]) + 1; /* start point */ | |||
1603 | BUCKET_B(c0, c0)(bucket_B[((c0) << 8) | (c0)]) = i; /* end point */ | |||
1604 | } | |||
1605 | } | |||
1606 | ||||
1607 | return m; | |||
1608 | } | |||
1609 | ||||
1610 | /* Constructs the suffix array by using the sorted order of type B* suffixes. */ | |||
1611 | static | |||
1612 | void | |||
1613 | construct_SA(const unsigned char *T, int *SA, | |||
1614 | int *bucket_A, int *bucket_B, | |||
1615 | int n, int m) { | |||
1616 | int *i, *j, *k; | |||
1617 | int s; | |||
1618 | int c0, c1, c2; | |||
1619 | ||||
1620 | if(0 < m) { | |||
1621 | /* Construct the sorted order of type B suffixes by using | |||
1622 | the sorted order of type B* suffixes. */ | |||
1623 | for(c1 = ALPHABET_SIZE(256) - 2; 0 <= c1; --c1) { | |||
1624 | /* Scan the suffix array from right to left. */ | |||
1625 | for(i = SA + BUCKET_BSTAR(c1, c1 + 1)(bucket_B[((c1) << 8) | (c1 + 1)]), | |||
1626 | j = SA + BUCKET_A(c1 + 1)bucket_A[(c1 + 1)] - 1, k = NULL((void*)0), c2 = -1; | |||
1627 | i <= j; | |||
1628 | --j) { | |||
1629 | if(0 < (s = *j)) { | |||
1630 | assert(T[s] == c1)((void) (0)); | |||
1631 | assert(((s + 1) < n) && (T[s] <= T[s + 1]))((void) (0)); | |||
1632 | assert(T[s - 1] <= T[s])((void) (0)); | |||
1633 | *j = ~s; | |||
1634 | c0 = T[--s]; | |||
1635 | if((0 < s) && (T[s - 1] > c0)) { s = ~s; } | |||
1636 | if(c0 != c2) { | |||
1637 | if(0 <= c2) { BUCKET_B(c2, c1)(bucket_B[((c1) << 8) | (c2)]) = k - SA; } | |||
1638 | k = SA + BUCKET_B(c2 = c0, c1)(bucket_B[((c1) << 8) | (c2 = c0)]); | |||
1639 | } | |||
1640 | assert(k < j)((void) (0)); assert(k != NULL)((void) (0)); | |||
1641 | *k-- = s; | |||
1642 | } else { | |||
1643 | assert(((s == 0) && (T[s] == c1)) || (s < 0))((void) (0)); | |||
1644 | *j = ~s; | |||
1645 | } | |||
1646 | } | |||
1647 | } | |||
1648 | } | |||
1649 | ||||
1650 | /* Construct the suffix array by using | |||
1651 | the sorted order of type B suffixes. */ | |||
1652 | k = SA + BUCKET_A(c2 = T[n - 1])bucket_A[(c2 = T[n - 1])]; | |||
1653 | *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1); | |||
1654 | /* Scan the suffix array from left to right. */ | |||
1655 | for(i = SA, j = SA + n; i < j; ++i) { | |||
1656 | if(0 < (s = *i)) { | |||
1657 | assert(T[s - 1] >= T[s])((void) (0)); | |||
1658 | c0 = T[--s]; | |||
1659 | if((s == 0) || (T[s - 1] < c0)) { s = ~s; } | |||
1660 | if(c0 != c2) { | |||
1661 | BUCKET_A(c2)bucket_A[(c2)] = k - SA; | |||
1662 | k = SA + BUCKET_A(c2 = c0)bucket_A[(c2 = c0)]; | |||
1663 | } | |||
1664 | assert(i < k)((void) (0)); | |||
1665 | *k++ = s; | |||
1666 | } else { | |||
1667 | assert(s < 0)((void) (0)); | |||
1668 | *i = ~s; | |||
1669 | } | |||
1670 | } | |||
1671 | } | |||
1672 | ||||
1673 | /* Constructs the burrows-wheeler transformed string directly | |||
1674 | by using the sorted order of type B* suffixes. */ | |||
1675 | static | |||
1676 | int | |||
1677 | construct_BWT(const unsigned char *T, int *SA, | |||
1678 | int *bucket_A, int *bucket_B, | |||
1679 | int n, int m) { | |||
1680 | int *i, *j, *k, *orig; | |||
1681 | int s; | |||
1682 | int c0, c1, c2; | |||
1683 | ||||
1684 | if(0 < m) { | |||
1685 | /* Construct the sorted order of type B suffixes by using | |||
1686 | the sorted order of type B* suffixes. */ | |||
1687 | for(c1 = ALPHABET_SIZE(256) - 2; 0 <= c1; --c1) { | |||
1688 | /* Scan the suffix array from right to left. */ | |||
1689 | for(i = SA + BUCKET_BSTAR(c1, c1 + 1)(bucket_B[((c1) << 8) | (c1 + 1)]), | |||
1690 | j = SA + BUCKET_A(c1 + 1)bucket_A[(c1 + 1)] - 1, k = NULL((void*)0), c2 = -1; | |||
1691 | i <= j; | |||
1692 | --j) { | |||
1693 | if(0 < (s = *j)) { | |||
1694 | assert(T[s] == c1)((void) (0)); | |||
1695 | assert(((s + 1) < n) && (T[s] <= T[s + 1]))((void) (0)); | |||
1696 | assert(T[s - 1] <= T[s])((void) (0)); | |||
1697 | c0 = T[--s]; | |||
1698 | *j = ~((int)c0); | |||
1699 | if((0 < s) && (T[s - 1] > c0)) { s = ~s; } | |||
1700 | if(c0 != c2) { | |||
1701 | if(0 <= c2) { BUCKET_B(c2, c1)(bucket_B[((c1) << 8) | (c2)]) = k - SA; } | |||
1702 | k = SA + BUCKET_B(c2 = c0, c1)(bucket_B[((c1) << 8) | (c2 = c0)]); | |||
1703 | } | |||
1704 | assert(k < j)((void) (0)); assert(k != NULL)((void) (0)); | |||
1705 | *k-- = s; | |||
1706 | } else if(s != 0) { | |||
1707 | *j = ~s; | |||
1708 | #ifndef NDEBUG1 | |||
1709 | } else { | |||
1710 | assert(T[s] == c1)((void) (0)); | |||
1711 | #endif | |||
1712 | } | |||
1713 | } | |||
1714 | } | |||
1715 | } | |||
1716 | ||||
1717 | /* Construct the BWTed string by using | |||
1718 | the sorted order of type B suffixes. */ | |||
1719 | k = SA + BUCKET_A(c2 = T[n - 1])bucket_A[(c2 = T[n - 1])]; | |||
1720 | *k++ = (T[n - 2] < c2) ? ~((int)T[n - 2]) : (n - 1); | |||
1721 | /* Scan the suffix array from left to right. */ | |||
1722 | for(i = SA, j = SA + n, orig = SA; i < j; ++i) { | |||
1723 | if(0 < (s = *i)) { | |||
1724 | assert(T[s - 1] >= T[s])((void) (0)); | |||
1725 | c0 = T[--s]; | |||
1726 | *i = c0; | |||
1727 | if((0 < s) && (T[s - 1] < c0)) { s = ~((int)T[s - 1]); } | |||
1728 | if(c0 != c2) { | |||
1729 | BUCKET_A(c2)bucket_A[(c2)] = k - SA; | |||
1730 | k = SA + BUCKET_A(c2 = c0)bucket_A[(c2 = c0)]; | |||
1731 | } | |||
1732 | assert(i < k)((void) (0)); | |||
1733 | *k++ = s; | |||
1734 | } else if(s != 0) { | |||
1735 | *i = ~s; | |||
1736 | } else { | |||
1737 | orig = i; | |||
1738 | } | |||
1739 | } | |||
1740 | ||||
1741 | return orig - SA; | |||
1742 | } | |||
1743 | ||||
1744 | /* Constructs the burrows-wheeler transformed string directly | |||
1745 | by using the sorted order of type B* suffixes. */ | |||
1746 | static | |||
1747 | int | |||
1748 | construct_BWT_indexes(const unsigned char *T, int *SA, | |||
1749 | int *bucket_A, int *bucket_B, | |||
1750 | int n, int m, | |||
1751 | unsigned char * num_indexes, int * indexes) { | |||
1752 | int *i, *j, *k, *orig; | |||
1753 | int s; | |||
1754 | int c0, c1, c2; | |||
1755 | ||||
1756 | int mod = n / 8; | |||
1757 | { | |||
1758 | mod |= mod >> 1; mod |= mod >> 2; | |||
1759 | mod |= mod >> 4; mod |= mod >> 8; | |||
1760 | mod |= mod >> 16; mod >>= 1; | |||
1761 | ||||
1762 | *num_indexes = (unsigned char)((n - 1) / (mod + 1)); | |||
1763 | } | |||
1764 | ||||
1765 | if(0 < m) { | |||
1766 | /* Construct the sorted order of type B suffixes by using | |||
1767 | the sorted order of type B* suffixes. */ | |||
1768 | for(c1 = ALPHABET_SIZE(256) - 2; 0 <= c1; --c1) { | |||
1769 | /* Scan the suffix array from right to left. */ | |||
1770 | for(i = SA + BUCKET_BSTAR(c1, c1 + 1)(bucket_B[((c1) << 8) | (c1 + 1)]), | |||
1771 | j = SA + BUCKET_A(c1 + 1)bucket_A[(c1 + 1)] - 1, k = NULL((void*)0), c2 = -1; | |||
1772 | i <= j; | |||
1773 | --j) { | |||
1774 | if(0 < (s = *j)) { | |||
1775 | assert(T[s] == c1)((void) (0)); | |||
1776 | assert(((s + 1) < n) && (T[s] <= T[s + 1]))((void) (0)); | |||
1777 | assert(T[s - 1] <= T[s])((void) (0)); | |||
1778 | ||||
1779 | if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = j - SA; | |||
1780 | ||||
1781 | c0 = T[--s]; | |||
1782 | *j = ~((int)c0); | |||
1783 | if((0 < s) && (T[s - 1] > c0)) { s = ~s; } | |||
1784 | if(c0 != c2) { | |||
1785 | if(0 <= c2) { BUCKET_B(c2, c1)(bucket_B[((c1) << 8) | (c2)]) = k - SA; } | |||
1786 | k = SA + BUCKET_B(c2 = c0, c1)(bucket_B[((c1) << 8) | (c2 = c0)]); | |||
1787 | } | |||
1788 | assert(k < j)((void) (0)); assert(k != NULL)((void) (0)); | |||
1789 | *k-- = s; | |||
| ||||
1790 | } else if(s != 0) { | |||
1791 | *j = ~s; | |||
1792 | #ifndef NDEBUG1 | |||
1793 | } else { | |||
1794 | assert(T[s] == c1)((void) (0)); | |||
1795 | #endif | |||
1796 | } | |||
1797 | } | |||
1798 | } | |||
1799 | } | |||
1800 | ||||
1801 | /* Construct the BWTed string by using | |||
1802 | the sorted order of type B suffixes. */ | |||
1803 | k = SA + BUCKET_A(c2 = T[n - 1])bucket_A[(c2 = T[n - 1])]; | |||
1804 | if (T[n - 2] < c2) { | |||
1805 | if (((n - 1) & mod) == 0) indexes[(n - 1) / (mod + 1) - 1] = k - SA; | |||
1806 | *k++ = ~((int)T[n - 2]); | |||
1807 | } | |||
1808 | else { | |||
1809 | *k++ = n - 1; | |||
1810 | } | |||
1811 | ||||
1812 | /* Scan the suffix array from left to right. */ | |||
1813 | for(i = SA, j = SA + n, orig = SA; i < j; ++i) { | |||
1814 | if(0 < (s = *i)) { | |||
1815 | assert(T[s - 1] >= T[s])((void) (0)); | |||
1816 | ||||
1817 | if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = i - SA; | |||
1818 | ||||
1819 | c0 = T[--s]; | |||
1820 | *i = c0; | |||
1821 | if(c0 != c2) { | |||
1822 | BUCKET_A(c2)bucket_A[(c2)] = k - SA; | |||
1823 | k = SA + BUCKET_A(c2 = c0)bucket_A[(c2 = c0)]; | |||
1824 | } | |||
1825 | assert(i < k)((void) (0)); | |||
1826 | if((0 < s) && (T[s - 1] < c0)) { | |||
1827 | if ((s & mod) == 0) indexes[s / (mod + 1) - 1] = k - SA; | |||
1828 | *k++ = ~((int)T[s - 1]); | |||
1829 | } else | |||
1830 | *k++ = s; | |||
1831 | } else if(s != 0) { | |||
1832 | *i = ~s; | |||
1833 | } else { | |||
1834 | orig = i; | |||
1835 | } | |||
1836 | } | |||
1837 | ||||
1838 | return orig - SA; | |||
1839 | } | |||
1840 | ||||
1841 | ||||
1842 | /*---------------------------------------------------------------------------*/ | |||
1843 | ||||
1844 | /*- Function -*/ | |||
1845 | ||||
1846 | int | |||
1847 | divsufsort(const unsigned char *T, int *SA, int n, int openMP) { | |||
1848 | int *bucket_A, *bucket_B; | |||
1849 | int m; | |||
1850 | int err = 0; | |||
1851 | ||||
1852 | /* Check arguments. */ | |||
1853 | if((T == NULL((void*)0)) || (SA == NULL((void*)0)) || (n < 0)) { return -1; } | |||
1854 | else if(n == 0) { return 0; } | |||
1855 | else if(n == 1) { SA[0] = 0; return 0; } | |||
1856 | else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; } | |||
1857 | ||||
1858 | bucket_A = (int *)malloc(BUCKET_A_SIZE((256)) * sizeof(int)); | |||
1859 | bucket_B = (int *)malloc(BUCKET_B_SIZE((256) * (256)) * sizeof(int)); | |||
1860 | ||||
1861 | /* Suffixsort. */ | |||
1862 | if((bucket_A != NULL((void*)0)) && (bucket_B != NULL((void*)0))) { | |||
1863 | m = sort_typeBstar(T, SA, bucket_A, bucket_B, n, openMP); | |||
1864 | construct_SA(T, SA, bucket_A, bucket_B, n, m); | |||
1865 | } else { | |||
1866 | err = -2; | |||
1867 | } | |||
1868 | ||||
1869 | free(bucket_B); | |||
1870 | free(bucket_A); | |||
1871 | ||||
1872 | return err; | |||
1873 | } | |||
1874 | ||||
1875 | int | |||
1876 | divbwt(const unsigned char *T, unsigned char *U, int *A, int n, unsigned char * num_indexes, int * indexes, int openMP) { | |||
1877 | int *B; | |||
1878 | int *bucket_A, *bucket_B; | |||
1879 | int m, pidx, i; | |||
1880 | ||||
1881 | /* Check arguments. */ | |||
1882 | if((T == NULL((void*)0)) || (U == NULL((void*)0)) || (n < 0)) { return -1; } | |||
| ||||
1883 | else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; } | |||
1884 | ||||
1885 | if((B = A) == NULL((void*)0)) { B = (int *)malloc((size_t)(n + 1) * sizeof(int)); } | |||
1886 | bucket_A = (int *)malloc(BUCKET_A_SIZE((256)) * sizeof(int)); | |||
1887 | bucket_B = (int *)malloc(BUCKET_B_SIZE((256) * (256)) * sizeof(int)); | |||
1888 | ||||
1889 | /* Burrows-Wheeler Transform. */ | |||
1890 | if((B != NULL((void*)0)) && (bucket_A != NULL((void*)0)) && (bucket_B != NULL((void*)0))) { | |||
1891 | m = sort_typeBstar(T, B, bucket_A, bucket_B, n, openMP); | |||
1892 | ||||
1893 | if (num_indexes == NULL((void*)0) || indexes == NULL((void*)0)) { | |||
1894 | pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m); | |||
1895 | } else { | |||
1896 | pidx = construct_BWT_indexes(T, B, bucket_A, bucket_B, n, m, num_indexes, indexes); | |||
1897 | } | |||
1898 | ||||
1899 | /* Copy to output string. */ | |||
1900 | U[0] = T[n - 1]; | |||
1901 | for(i = 0; i < pidx; ++i) { U[i + 1] = (unsigned char)B[i]; } | |||
1902 | for(i += 1; i < n; ++i) { U[i] = (unsigned char)B[i]; } | |||
1903 | pidx += 1; | |||
1904 | } else { | |||
1905 | pidx = -2; | |||
1906 | } | |||
1907 | ||||
1908 | free(bucket_B); | |||
1909 | free(bucket_A); | |||
1910 | if(A == NULL((void*)0)) { free(B); } | |||
1911 | ||||
1912 | return pidx; | |||
1913 | } |