2 * utils.c for libdivsufsort
3 * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
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
14 * The above copyright notice and this permission notice shall be
15 * included in all copies or substantial portions of the Software.
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.
27 #include "divsufsort_private.h"
30 /*- Private Function -*/
32 /* Binary search for inverse bwt. */
35 binarysearch_lower(const saidx_t *A, saidx_t size, saidx_t value) {
37 for(i = 0, half = size >> 1;
39 size = half, half >>= 1) {
40 if(A[i + half] < value) {
42 half -= (size & 1) ^ 1;
51 /* Burrows-Wheeler transform. */
53 bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *SA,
54 saidx_t n, saidx_t *idx) {
55 saidx_t *A, i, j, p, t;
58 /* Check arguments. */
59 if((T == NULL) || (U == NULL) || (n < 0) || (idx == NULL)) { return -1; }
61 if(n == 1) { U[0] = T[0]; }
66 if((A = SA) == NULL) {
67 i = divbwt(T, U, NULL, n);
68 if(0 <= i) { *idx = i; i = 0; }
75 for(i = 0, j = 0; i < n; ++i) {
80 U[j] = (j <= p) ? T[p] : (sauchar_t)A[p];
90 U[j] = (j <= p) ? T[p] : (sauchar_t)A[p];
97 for(i = 0; A[i] != 0; ++i) { U[i + 1] = T[A[i] - 1]; }
99 for(++i; i < n; ++i) { U[i] = T[A[i] - 1]; }
103 /* Deallocate memory. */
110 /* Inverse Burrows-Wheeler transform. */
112 inverse_bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *A,
113 saidx_t n, saidx_t idx) {
114 saidx_t C[ALPHABET_SIZE];
115 sauchar_t D[ALPHABET_SIZE];
120 /* Check arguments. */
121 if((T == NULL) || (U == NULL) || (n < 0) || (idx < 0) ||
122 (n < idx) || ((0 < n) && (idx == 0))) {
125 if(n <= 1) { return 0; }
127 if((B = A) == NULL) {
128 /* Allocate n*sizeof(saidx_t) bytes of memory. */
129 if((B = (saidx_t *)malloc((size_t)n * sizeof(saidx_t))) == NULL) { return -2; }
132 /* Inverse BW transform. */
133 for(c = 0; c < ALPHABET_SIZE; ++c) { C[c] = 0; }
134 for(i = 0; i < n; ++i) { ++C[T[i]]; }
135 for(c = 0, d = 0, i = 0; c < ALPHABET_SIZE; ++c) {
139 D[d++] = (sauchar_t)c;
143 for(i = 0; i < idx; ++i) { B[C[T[i]]++] = i; }
144 for( ; i < n; ++i) { B[C[T[i]]++] = i + 1; }
145 for(c = 0; c < d; ++c) { C[c] = C[D[c]]; }
146 for(i = 0, p = idx; i < n; ++i) {
147 U[i] = D[binarysearch_lower(C, d, p)];
152 /* Deallocate memory. */
159 /* Checks the suffix array SA of the string T. */
161 sufcheck(const sauchar_t *T, const saidx_t *SA,
162 saidx_t n, saint_t verbose) {
163 saidx_t C[ALPHABET_SIZE];
167 if(verbose) { fprintf(stderr, "sufcheck: "); }
169 /* Check arguments. */
170 if((T == NULL) || (SA == NULL) || (n < 0)) {
171 if(verbose) { fprintf(stderr, "Invalid arguments.\n"); }
175 if(verbose) { fprintf(stderr, "Done.\n"); }
179 /* check range: [0..n-1] */
180 for(i = 0; i < n; ++i) {
181 if((SA[i] < 0) || (n <= SA[i])) {
183 fprintf(stderr, "Out of the range [0,%" PRIdSAIDX_T "].\n"
184 " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "\n",
191 /* check first characters. */
192 for(i = 1; i < n; ++i) {
193 if(T[SA[i - 1]] > T[SA[i]]) {
195 fprintf(stderr, "Suffixes in wrong order.\n"
196 " T[SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "]=%d"
197 " > T[SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "]=%d\n",
198 i - 1, SA[i - 1], T[SA[i - 1]], i, SA[i], T[SA[i]]);
204 /* check suffixes. */
205 for(i = 0; i < ALPHABET_SIZE; ++i) { C[i] = 0; }
206 for(i = 0; i < n; ++i) { ++C[T[i]]; }
207 for(i = 0, p = 0; i < ALPHABET_SIZE; ++i) {
215 for(i = 0; i < n; ++i) {
224 if((t < 0) || (p != SA[t])) {
226 fprintf(stderr, "Suffix in wrong position.\n"
227 " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T " or\n"
228 " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "\n",
229 t, (0 <= t) ? SA[t] : -1, i, SA[i]);
235 if((n <= C[c]) || (T[SA[C[c]]] != c)) { C[c] = -1; }
239 if(1 <= verbose) { fprintf(stderr, "Done.\n"); }
246 _compare(const sauchar_t *T, saidx_t Tsize,
247 const sauchar_t *P, saidx_t Psize,
248 saidx_t suf, saidx_t *match) {
251 for(i = suf + *match, j = *match, r = 0;
252 (i < Tsize) && (j < Psize) && ((r = T[i] - P[j]) == 0); ++i, ++j) { }
254 return (r == 0) ? -(j != Psize) : r;
257 /* Search for the pattern P in the string T. */
259 sa_search(const sauchar_t *T, saidx_t Tsize,
260 const sauchar_t *P, saidx_t Psize,
261 const saidx_t *SA, saidx_t SAsize,
263 saidx_t size, lsize, rsize, half;
264 saidx_t match, lmatch, rmatch;
265 saidx_t llmatch, lrmatch, rlmatch, rrmatch;
269 if(idx != NULL) { *idx = -1; }
270 if((T == NULL) || (P == NULL) || (SA == NULL) ||
271 (Tsize < 0) || (Psize < 0) || (SAsize < 0)) { return -1; }
272 if((Tsize == 0) || (SAsize == 0)) { return 0; }
273 if(Psize == 0) { if(idx != NULL) { *idx = 0; } return SAsize; }
275 for(i = j = k = 0, lmatch = rmatch = 0, size = SAsize, half = size >> 1;
277 size = half, half >>= 1) {
278 match = MIN(lmatch, rmatch);
279 r = _compare(T, Tsize, P, Psize, SA[i + half], &match);
282 half -= (size & 1) ^ 1;
287 lsize = half, j = i, rsize = size - half - 1, k = i + half + 1;
290 for(llmatch = lmatch, lrmatch = match, half = lsize >> 1;
292 lsize = half, half >>= 1) {
293 lmatch = MIN(llmatch, lrmatch);
294 r = _compare(T, Tsize, P, Psize, SA[j + half], &lmatch);
297 half -= (lsize & 1) ^ 1;
305 for(rlmatch = match, rrmatch = rmatch, half = rsize >> 1;
307 rsize = half, half >>= 1) {
308 rmatch = MIN(rlmatch, rrmatch);
309 r = _compare(T, Tsize, P, Psize, SA[k + half], &rmatch);
312 half -= (rsize & 1) ^ 1;
323 if(idx != NULL) { *idx = (0 < (k - j)) ? j : i; }
327 /* Search for the character c in the string T. */
329 sa_simplesearch(const sauchar_t *T, saidx_t Tsize,
330 const saidx_t *SA, saidx_t SAsize,
331 saint_t c, saidx_t *idx) {
332 saidx_t size, lsize, rsize, half;
336 if(idx != NULL) { *idx = -1; }
337 if((T == NULL) || (SA == NULL) || (Tsize < 0) || (SAsize < 0)) { return -1; }
338 if((Tsize == 0) || (SAsize == 0)) { return 0; }
340 for(i = j = k = 0, size = SAsize, half = size >> 1;
342 size = half, half >>= 1) {
344 r = (p < Tsize) ? T[p] - c : -1;
347 half -= (size & 1) ^ 1;
349 lsize = half, j = i, rsize = size - half - 1, k = i + half + 1;
352 for(half = lsize >> 1;
354 lsize = half, half >>= 1) {
356 r = (p < Tsize) ? T[p] - c : -1;
359 half -= (lsize & 1) ^ 1;
364 for(half = rsize >> 1;
366 rsize = half, half >>= 1) {
368 r = (p < Tsize) ? T[p] - c : -1;
371 half -= (rsize & 1) ^ 1;
379 if(idx != NULL) { *idx = (0 < (k - j)) ? j : i; }