aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CMakeLists.txt22
-rw-r--r--CMakeModules/ProjectCPack.cmake38
-rw-r--r--COPYING5
-rw-r--r--examples/bwt.c4
-rw-r--r--examples/mksary.c6
-rw-r--r--examples/sasearch.c3
-rw-r--r--examples/suftest.c6
-rw-r--r--examples/unbwt.c4
-rw-r--r--include/config.h.cmake5
-rw-r--r--include/divsufsort.h.cmake6
-rw-r--r--include/divsufsort_private.h116
-rw-r--r--lib/CMakeLists.txt2
-rw-r--r--lib/divsufsort.c260
-rw-r--r--lib/sssort.c815
-rw-r--r--lib/substringsort.c619
-rw-r--r--lib/trsort.c662
-rw-r--r--lib/utils.c219
-rw-r--r--pkgconfig/CMakeLists.txt3
-rw-r--r--pkgconfig/libdivsufsort.pc.cmake11
19 files changed, 1525 insertions, 1281 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 142c4bc..db21f0b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,4 +1,4 @@
-### cmake file for building sais Package ###
+### cmake file for building libdivsufsort Package ###
cmake_minimum_required(VERSION 2.4.2)
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules")
include(AppendCompilerFlags)
@@ -17,6 +17,10 @@ endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.svn")
## Project information ##
project(libdivsufsort C)
+set(PROJECT_VENDOR "Yuta Mori")
+set(PROJECT_CONTACT "yuta.256@gmail.com")
+set(PROJECT_URL "http://libdivsufsort.googlecode.com/")
+set(PROJECT_DESCRIPTION "A lightweight suffix sorting library")
file(READ "${CMAKE_CURRENT_SOURCE_DIR}/VERSION" PROJECT_VERSION_FULL)
string(REGEX REPLACE "[\n\r]" "" PROJECT_VERSION_FULL "${PROJECT_VERSION_FULL}")
string(REGEX REPLACE "^([0-9]+)\\.[0-9]+\\.[0-9]+$" "\\1" PROJECT_VERSION_MAJOR "${PROJECT_VERSION_FULL}")
@@ -32,9 +36,15 @@ if(SVN_REVISION)
set(PROJECT_VERSION_FULL "${PROJECT_VERSION_FULL}-svn-r${SVN_REVISION}")
endif(SVN_REVISION)
+## CPack configuration ##
+set(CPACK_GENERATOR "TGZ;TBZ2;ZIP")
+set(CPACK_SOURCE_GENERATOR "TGZ;TBZ2;ZIP")
+include(ProjectCPack)
+
## Project options ##
option(BUILD_SHARED_LIBS "Set to OFF to build static libraries" ON)
option(BUILD_EXAMPLES "Build examples" ON)
+option(USE_OPENMP "Use OpenMP for parallelization" OFF)
## Build type ##
if(NOT CMAKE_BUILD_TYPE)
@@ -47,6 +57,9 @@ endif(NOT CMAKE_BUILD_TYPE)
if(MSVC)
append_c_compiler_flags("/W4" "VC" CMAKE_C_FLAGS)
append_c_compiler_flags("/Oi;/Ot;/Ox;/Oy" "VC" CMAKE_C_FLAGS_RELEASE)
+ if(USE_OPENMP)
+ append_c_compiler_flags("/openmp" "VC" CMAKE_C_FLAGS)
+ endif(USE_OPENMP)
elseif(BORLAND)
append_c_compiler_flags("-w" "BCC" CMAKE_C_FLAGS)
append_c_compiler_flags("-Oi;-Og;-Os;-Ov;-Ox" "BCC" CMAKE_C_FLAGS_RELEASE)
@@ -54,9 +67,15 @@ else(MSVC)
if(CMAKE_COMPILER_IS_GNUCC)
append_c_compiler_flags("-Wall" "GCC" CMAKE_C_FLAGS)
append_c_compiler_flags("-fomit-frame-pointer" "GCC" CMAKE_C_FLAGS_RELEASE)
+ if(USE_OPENMP)
+ append_c_compiler_flags("-fopenmp" "GCC" CMAKE_C_FLAGS)
+ endif(USE_OPENMP)
else(CMAKE_COMPILER_IS_GNUCC)
append_c_compiler_flags("-Wall" "UNKNOWN" CMAKE_C_FLAGS)
append_c_compiler_flags("-fomit-frame-pointer" "UNKNOWN" CMAKE_C_FLAGS_RELEASE)
+ if(USE_OPENMP)
+ append_c_compiler_flags("-fopenmp;-openmp;-omp" "UNKNOWN" CMAKE_C_FLAGS)
+ endif(USE_OPENMP)
endif(CMAKE_COMPILER_IS_GNUCC)
endif(MSVC)
@@ -64,6 +83,7 @@ endif(MSVC)
add_definitions(-DHAVE_CONFIG_H=1 -D__STDC_LIMIT_MACROS -D__STDC_CONSTANT_MACROS -D__STDC_FORMAT_MACROS)
## Add subdirectories ##
+add_subdirectory(pkgconfig)
add_subdirectory(include)
add_subdirectory(lib)
if(BUILD_EXAMPLES)
diff --git a/CMakeModules/ProjectCPack.cmake b/CMakeModules/ProjectCPack.cmake
new file mode 100644
index 0000000..35694b1
--- /dev/null
+++ b/CMakeModules/ProjectCPack.cmake
@@ -0,0 +1,38 @@
+# If the cmake version includes cpack, use it
+IF(EXISTS "${CMAKE_ROOT}/Modules/CPack.cmake")
+ SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY "${PROJECT_DESCRIPTION}")
+ SET(CPACK_PACKAGE_VENDOR "${PROJECT_VENDOR}")
+ SET(CPACK_PACKAGE_DESCRIPTION_FILE "${CMAKE_CURRENT_SOURCE_DIR}/COPYING")
+ SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/COPYING")
+ SET(CPACK_PACKAGE_VERSION_MAJOR "${PROJECT_VERSION_MAJOR}")
+ SET(CPACK_PACKAGE_VERSION_MINOR "${PROJECT_VERSION_MINOR}")
+ SET(CPACK_PACKAGE_VERSION_PATCH "${PROJECT_VERSION_PATCH}")
+# SET(CPACK_PACKAGE_INSTALL_DIRECTORY "${PROJECT_NAME} ${PROJECT_VERSION}")
+ SET(CPACK_SOURCE_PACKAGE_FILE_NAME "${PROJECT_NAME}-${PROJECT_VERSION_FULL}")
+
+ IF(NOT DEFINED CPACK_SYSTEM_NAME)
+ SET(CPACK_SYSTEM_NAME "${CMAKE_SYSTEM_NAME}-${CMAKE_SYSTEM_PROCESSOR}")
+ ENDIF(NOT DEFINED CPACK_SYSTEM_NAME)
+
+ IF(${CPACK_SYSTEM_NAME} MATCHES Windows)
+ IF(CMAKE_CL_64)
+ SET(CPACK_SYSTEM_NAME win64-${CMAKE_SYSTEM_PROCESSOR})
+ ELSE(CMAKE_CL_64)
+ SET(CPACK_SYSTEM_NAME win32-${CMAKE_SYSTEM_PROCESSOR})
+ ENDIF(CMAKE_CL_64)
+ ENDIF(${CPACK_SYSTEM_NAME} MATCHES Windows)
+
+ IF(NOT DEFINED CPACK_PACKAGE_FILE_NAME)
+ SET(CPACK_PACKAGE_FILE_NAME "${CPACK_SOURCE_PACKAGE_FILE_NAME}-${CPACK_SYSTEM_NAME}")
+ ENDIF(NOT DEFINED CPACK_PACKAGE_FILE_NAME)
+
+ SET(CPACK_PACKAGE_CONTACT "${PROJECT_CONTACT}")
+ IF(UNIX)
+ SET(CPACK_STRIP_FILES "")
+ SET(CPACK_SOURCE_STRIP_FILES "")
+# SET(CPACK_PACKAGE_EXECUTABLES "ccmake" "CMake")
+ ENDIF(UNIX)
+ SET(CPACK_SOURCE_IGNORE_FILES "/CVS/" "/build/" "/\\\\.build/" "/\\\\.svn/" "~$")
+ # include CPack model once all variables are set
+ INCLUDE(CPack)
+ENDIF(EXISTS "${CMAKE_ROOT}/Modules/CPack.cmake")
diff --git a/COPYING b/COPYING
index 3bad2dc..a0a6432 100644
--- a/COPYING
+++ b/COPYING
@@ -1,3 +1,5 @@
+The libdivsufsort copyright is as follows:
+
Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
Permission is hereby granted, free of charge, to any person
@@ -20,3 +22,6 @@ HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.
+
+See also the libdivsufsort web site:
+ http://libdivsufsort.googlecode.com/ for more information.
diff --git a/examples/bwt.c b/examples/bwt.c
index 41591f8..842712d 100644
--- a/examples/bwt.c
+++ b/examples/bwt.c
@@ -101,9 +101,9 @@ main(int argc, const char *argv[]) {
rewind(stdin);
}
- /* Allocate blocksize+4*(blocksize+1) bytes of memory. */
+ /* Allocate 5blocksize bytes of memory. */
if(((T = malloc(blocksize * sizeof(sauchar_t))) == NULL) ||
- ((SA = malloc((blocksize + 1) * sizeof(saidx_t))) == NULL)) {
+ ((SA = malloc(blocksize * sizeof(saidx_t))) == NULL)) {
fprintf(stderr, "%s: Cannot allocate memory.\n", argv[0]);
exit(EXIT_FAILURE);
}
diff --git a/examples/mksary.c b/examples/mksary.c
index a548582..cdb6d35 100644
--- a/examples/mksary.c
+++ b/examples/mksary.c
@@ -90,9 +90,9 @@ main(int argc, const char *argv[]) {
rewind(ifp);
#endif
- /* Allocate n+4(n+1) bytes of memory. */
+ /* Allocate 5n bytes of memory. */
if(((T = malloc(n * sizeof(sauchar_t))) == NULL) ||
- ((SA = malloc((n + 1) * sizeof(saidx_t))) == NULL)) {
+ ((SA = malloc(n * sizeof(saidx_t))) == NULL)) {
fprintf(stderr, "%s: Cannot allocate memory.\n", argv[0]);
exit(EXIT_FAILURE);
}
@@ -117,7 +117,7 @@ main(int argc, const char *argv[]) {
(double)(finish - start) / (double)CLOCKS_PER_SEC);
/* Write the suffix array. */
- if(fwrite(SA, sizeof(saidx_t), n + 1, ofp) != (n + 1)) {
+ if(fwrite(SA, sizeof(saidx_t), n, ofp) != n) {
fprintf(stderr, "%s: Cannot write to `%s': ", argv[0], argv[2]);
perror(NULL);
exit(EXIT_FAILURE);
diff --git a/examples/sasearch.c b/examples/sasearch.c
index bf8f7cc..026aff9 100644
--- a/examples/sasearch.c
+++ b/examples/sasearch.c
@@ -380,7 +380,8 @@ main(int argc, const char *argv[]) {
exit(EXIT_FAILURE);
}
/* Read (n + 1) * sizeof(saidx_t) bytes of data. */
- if(fread(SA, sizeof(saidx_t), size + 1, fp) != (size + 1)) {
+ SA[0] = size;
+ if(fread(SA + 1, sizeof(saidx_t), size, fp) != size) {
fprintf(stderr, "%s: %s `%s': ",
argv[0],
(ferror(fp) || !feof(fp)) ? "Cannot read from" : "Unexpected EOF in",
diff --git a/examples/suftest.c b/examples/suftest.c
index 9d9e2d4..90b7c83 100644
--- a/examples/suftest.c
+++ b/examples/suftest.c
@@ -94,9 +94,9 @@ main(int argc, const char *argv[]) {
rewind(fp);
#endif
- /* Allocate n+4(n+1) bytes of memory. */
+ /* Allocate 5n bytes of memory. */
if(((T = malloc(n * sizeof(sauchar_t))) == NULL) ||
- ((SA = malloc((n + 1) * sizeof(saidx_t))) == NULL)) {
+ ((SA = malloc(n * sizeof(saidx_t))) == NULL)) {
fprintf(stderr, "%s: Cannot allocate memory.\n", argv[0]);
exit(EXIT_FAILURE);
}
@@ -121,7 +121,7 @@ main(int argc, const char *argv[]) {
(double)(finish - start) / (double)CLOCKS_PER_SEC);
/* Check the suffix array. */
- if(sufcheck(T, SA, n, 3) != 0) {
+ if(sufcheck(T, SA, n, 1) != 0) {
exit(EXIT_FAILURE);
}
diff --git a/examples/unbwt.c b/examples/unbwt.c
index e989053..a0aae45 100644
--- a/examples/unbwt.c
+++ b/examples/unbwt.c
@@ -72,9 +72,9 @@ main(int argc, const char *argv[]) {
exit(EXIT_FAILURE);
}
- /* Allocate blocksize+4(blocksize+1) bytes of memory. */
+ /* Allocate 5blocksize bytes of memory. */
if(((T = malloc(blocksize * sizeof(sauchar_t))) == NULL) ||
- ((A = malloc((blocksize + 1) * sizeof(saidx_t))) == NULL)) {
+ ((A = malloc(blocksize * sizeof(saidx_t))) == NULL)) {
fprintf(stderr, "%s: Cannot allocate memory.\n", argv[0]);
exit(EXIT_FAILURE);
}
diff --git a/include/config.h.cmake b/include/config.h.cmake
index a2223da..f0a0159 100644
--- a/include/config.h.cmake
+++ b/include/config.h.cmake
@@ -49,6 +49,11 @@ extern "C" {
# define INLINE @INLINE@
#endif
+/** for VC++ warning **/
+#ifdef _MSC_VER
+#pragma warning(disable: 4127)
+#endif
+
#ifdef __cplusplus
} /* extern "C" */
diff --git a/include/divsufsort.h.cmake b/include/divsufsort.h.cmake
index 73feefe..45c7b04 100644
--- a/include/divsufsort.h.cmake
+++ b/include/divsufsort.h.cmake
@@ -67,7 +67,7 @@ typedef @SAINDEX_TYPE@ saidx_t;
/**
* Constructs the suffix array of a given string.
* @param T[0..n-1] The input string.
- * @param SA[0..n] The output array of suffixes.
+ * @param SA[0..n-1] The output array of suffixes.
* @param n The length of the given string.
* @return 0 if no error occurred, -1 or -2 otherwise.
*/
@@ -79,7 +79,7 @@ divsufsort(const sauchar_t *T, saidx_t *SA, saidx_t n);
* Constructs the burrows-wheeler transformed string of a given string.
* @param T[0..n-1] The input string.
* @param U[0..n-1] The output string. (can be T)
- * @param A[0..n] The temporary array. (can be NULL)
+ * @param A[0..n-1] The temporary array. (can be NULL)
* @param n The length of the given string.
* @return The primary index if no error occurred, -1 or -2 otherwise.
*/
@@ -113,7 +113,7 @@ inverse_bw_transform(const sauchar_t *T, sauchar_t *U,
/**
* Checks the correctness of a given suffix array.
* @param T[0..n-1] The input string.
- * @param SA[0..n] The input suffix array.
+ * @param SA[0..n-1] The input suffix array.
* @param n The length of the given string.
* @param verbose The verbose mode.
* @return 0 if no error occurred.
diff --git a/include/divsufsort_private.h b/include/divsufsort_private.h
index 538e2c0..0d7bd28 100644
--- a/include/divsufsort_private.h
+++ b/include/divsufsort_private.h
@@ -36,16 +36,14 @@ extern "C" {
#endif
#include <assert.h>
#include <stdio.h>
-#if STDC_HEADERS
-# include <stdlib.h>
+#if HAVE_STRING_H
# include <string.h>
-#else
-# if HAVE_STDLIB_H
-# include <stdlib.h>
-# endif
-# if HAVE_MEMORY_H
-# include <memory.h>
-# endif
+#endif
+#if HAVE_STDLIB_H
+# include <stdlib.h>
+#endif
+#if HAVE_MEMORY_H
+# include <memory.h>
#endif
#if HAVE_STDDEF_H
# include <stddef.h>
@@ -73,83 +71,109 @@ extern "C" {
#if !defined(ALPHABET_SIZE)
# define ALPHABET_SIZE (UINT8_MAX + 1)
#endif
-#if defined(STACK_SIZE) && (STACK_SIZE < 32)
-# undef STACK_SIZE
-#endif
-#if !defined(STACK_SIZE)
-# define STACK_SIZE (64)
-#endif
-#if defined(LOCALMERGE_BUFFERSIZE) && (LOCALMERGE_BUFFERSIZE < 32)
-# undef LOCALMERGE_BUFFERSIZE
-#endif
-#if !defined(LOCALMERGE_BUFFERSIZE)
-# define LOCALMERGE_BUFFERSIZE (256)
-#endif
/* for divsufsort.c */
#define BUCKET_A_SIZE (ALPHABET_SIZE)
#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE)
/* for sssort.c */
-#define SS_INSERTIONSORT_THRESHOLD (8)
-#define SS_BLOCKSIZE (1024)
+#if defined(SS_INSERTIONSORT_THRESHOLD)
+# if SS_INSERTIONSORT_THRESHOLD < 1
+# undef SS_INSERTIONSORT_THRESHOLD
+# define SS_INSERTIONSORT_THRESHOLD 1
+# endif
+#else
+# define SS_INSERTIONSORT_THRESHOLD (8)
+#endif
+#if defined(SS_BLOCKSIZE)
+# if SS_BLOCKSIZE < 0
+# undef SS_BLOCKSIZE
+# define SS_BLOCKSIZE 0
+# elif 65536 <= SS_BLOCKSIZE
+# undef SS_BLOCKSIZE
+# define SS_BLOCKSIZE 65535
+# endif
+#else
+# define SS_BLOCKSIZE (1024)
+#endif
+/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */
+#if SS_BLOCKSIZE == 0
+# if defined(BUILD_DIVSUFSORT64)
+# define SS_MISORT_STACKSIZE (96)
+# else
+# define SS_MISORT_STACKSIZE (64)
+# endif
+#elif SS_BLOCKSIZE <= 4096
+# define SS_MISORT_STACKSIZE (16)
+#else
+# define SS_MISORT_STACKSIZE (24)
+#endif
+#if defined(BUILD_DIVSUFSORT64)
+# define SS_SMERGE_STACKSIZE (64)
+#else
+# define SS_SMERGE_STACKSIZE (32)
+#endif
/* for trsort.c */
-#define LS_INSERTIONSORT_THRESHOLD (8)
#define TR_INSERTIONSORT_THRESHOLD (8)
+#if defined(BUILD_DIVSUFSORT64)
+# define TR_STACKSIZE (96)
+#else
+# define TR_STACKSIZE (64)
+#endif
/*- Macros -*/
#ifndef SWAP
-# define SWAP(a,b) do { t = (a); (a) = (b); (b) = t; } while(0)
+# define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0)
#endif /* SWAP */
#ifndef MIN
-# define MIN(a,b) (((a) < (b)) ? (a) : (b))
+# define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b))
#endif /* MIN */
#ifndef MAX
-# define MAX(a,b) (((a) > (b)) ? (a) : (b))
+# define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b))
#endif /* MAX */
-#define STACK_PUSH3(_a, _b, _c)\
+#define STACK_PUSH(_a, _b, _c, _d)\
do {\
assert(ssize < STACK_SIZE);\
stack[ssize].a = (_a), stack[ssize].b = (_b),\
- stack[ssize++].c = (_c);\
+ stack[ssize].c = (_c), stack[ssize++].d = (_d);\
} while(0)
-#define STACK_PUSH(_a, _b, _c, _d)\
+#define STACK_PUSH5(_a, _b, _c, _d, _e)\
do {\
assert(ssize < STACK_SIZE);\
stack[ssize].a = (_a), stack[ssize].b = (_b),\
- stack[ssize].c = (_c), stack[ssize++].d = (_d);\
+ stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\
} while(0)
-#define STACK_POP3(_a, _b, _c)\
+#define STACK_POP(_a, _b, _c, _d)\
do {\
assert(0 <= ssize);\
if(ssize == 0) { return; }\
(_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
- (_c) = stack[ssize].c;\
+ (_c) = stack[ssize].c, (_d) = stack[ssize].d;\
} while(0)
-#define STACK_POP(_a, _b, _c, _d)\
+#define STACK_POP5(_a, _b, _c, _d, _e)\
do {\
assert(0 <= ssize);\
if(ssize == 0) { return; }\
(_a) = stack[--ssize].a, (_b) = stack[ssize].b,\
- (_c) = stack[ssize].c, (_d) = stack[ssize].d;\
+ (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\
} while(0)
/* for divsufsort.c */
-#define BUCKET_A(c0) bucket_A[(c0)]
+#define BUCKET_A(_c0) bucket_A[(_c0)]
#if ALPHABET_SIZE == 256
-#define BUCKET_B(c0, c1) (bucket_B[((c1) << 8) | (c0)])
-#define BUCKET_BSTAR(c0, c1) (bucket_B[((c0) << 8) | (c1)])
+#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)])
+#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)])
#else
-#define BUCKET_B(c0, c1) (bucket_B[(c1) * ALPHABET_SIZE + (c0)])
-#define BUCKET_BSTAR(c0, c1) (bucket_B[(c0) * ALPHABET_SIZE + (c1)])
+#define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)])
+#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)])
#endif
/*- Private Prototypes -*/
-/* substringsort.c */
+/* sssort.c */
void
-substringsort(const sauchar_t *Td, const saidx_t *PA,
- saidx_t *first, saidx_t *last,
- saidx_t *buf, saidx_t bufsize,
- saidx_t depth, saint_t lastsuffix);
+sssort(const sauchar_t *Td, const saidx_t *PA,
+ saidx_t *first, saidx_t *last,
+ saidx_t *buf, saidx_t bufsize,
+ saidx_t depth, saidx_t n, saint_t lastsuffix);
/* trsort.c */
void
trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth);
@@ -159,4 +183,4 @@ trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth);
} /* extern "C" */
#endif /* __cplusplus */
-#endif /* _DIVSUFSORT_PRIVATE_H_ */
+#endif /* _DIVSUFSORT_PRIVATE_H */
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
index e8c5295..cbcb427 100644
--- a/lib/CMakeLists.txt
+++ b/lib/CMakeLists.txt
@@ -2,7 +2,7 @@ include_directories("${CMAKE_CURRENT_SOURCE_DIR}/../include"
"${CMAKE_CURRENT_BINARY_DIR}/../include")
## libdivsufsort ##
-add_library(divsufsort divsufsort.c substringsort.c trsort.c utils.c)
+add_library(divsufsort divsufsort.c sssort.c trsort.c utils.c)
install(TARGETS divsufsort
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
diff --git a/lib/divsufsort.c b/lib/divsufsort.c
index ac14ea3..568803f 100644
--- a/lib/divsufsort.c
+++ b/lib/divsufsort.c
@@ -25,6 +25,9 @@
*/
#include "divsufsort_private.h"
+#ifdef _OPENMP
+# include <omp.h>
+#endif
/*- Private Functions -*/
@@ -32,12 +35,20 @@
/* Sorts suffixes of type B*. */
static
saidx_t
-_sort_typeBstar(const sauchar_t *T, saidx_t *SA,
- saidx_t *bucket_A, saidx_t *bucket_B,
- saidx_t n) {
+sort_typeBstar(const sauchar_t *T, saidx_t *SA,
+ saidx_t *bucket_A, saidx_t *bucket_B,
+ saidx_t n) {
saidx_t *PAb, *ISAb, *buf;
+#ifdef _OPENMP
+ saidx_t *curbuf;
+ saidx_t l;
+#endif
saidx_t i, j, k, t, m, bufsize;
saint_t c0, c1;
+#ifdef _OPENMP
+ saint_t d0, d1;
+ int tmp;
+#endif
/* Initialize bucket arrays. */
for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; }
@@ -46,16 +57,16 @@ _sort_typeBstar(const sauchar_t *T, saidx_t *SA,
/* Count the number of occurrences of the first one or two characters of each
type A, B and B* suffix. Moreover, store the beginning position of all
type B* suffixes into the array SA. */
- for(i = n - 1, m = n; 0 <= i;) {
+ for(i = n - 1, m = n, c0 = T[n - 1]; 0 <= i;) {
/* type A suffix. */
- do { ++BUCKET_A(T[i]); } while((0 <= --i) && (T[i] >= T[i + 1]));
+ do { ++BUCKET_A(c1 = c0); } while((0 <= --i) && ((c0 = T[i]) >= c1));
if(0 <= i) {
/* type B* suffix. */
- ++BUCKET_BSTAR(T[i], T[i + 1]);
+ ++BUCKET_BSTAR(c0, c1);
SA[--m] = i;
/* type B suffix. */
- for(--i; (0 <= i) && (T[i] <= T[i + 1]); --i) {
- ++BUCKET_B(T[i], T[i + 1]);
+ for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) {
+ ++BUCKET_B(c0, c1);
}
}
}
@@ -81,7 +92,6 @@ note:
if(0 < m) {
/* Sort the type B* suffixes by their first two characters. */
PAb = SA + n - m; ISAb = SA + m;
- PAb[m] = n - 2; /* for sentinel. */
for(i = m - 2; 0 <= i; --i) {
t = PAb[i], c0 = T[t], c1 = T[t + 1];
SA[--BUCKET_BSTAR(c0, c1)] = i;
@@ -90,23 +100,47 @@ note:
SA[--BUCKET_BSTAR(c0, c1)] = m - 1;
/* Sort the type B* substrings using sssort. */
- buf = SA + m, bufsize = n - (2 * m);
- if(bufsize <= LOCALMERGE_BUFFERSIZE) {
- if((buf = malloc(LOCALMERGE_BUFFERSIZE * sizeof(saidx_t))) == NULL) {
- return -1;
+#ifdef _OPENMP
+ tmp = omp_get_max_threads();
+ buf = SA + m, bufsize = (n - (2 * m)) / tmp;
+ c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m;
+#pragma omp parallel default(shared) private(curbuf, k, l, d0, d1, tmp)
+ {
+ tmp = omp_get_thread_num();
+ curbuf = buf + tmp * bufsize;
+ k = 0;
+ for(;;) {
+ #pragma omp critical(sssort_lock)
+ {
+ if(0 < (l = j)) {
+ d0 = c0, d1 = c1;
+ do {
+ k = BUCKET_BSTAR(d0, d1);
+ if(--d1 <= d0) {
+ d1 = ALPHABET_SIZE - 1;
+ if(--d0 < 0) { break; }
+ }
+ } while(((l - k) <= 1) && (0 < (l = k)));
+ c0 = d0, c1 = d1, j = k;
+ }
+ }
+ if(l == 0) { break; }
+ sssort(T, PAb, SA + k, SA + l,
+ curbuf, bufsize, 2, n, *(SA + k) == (m - 1));
}
- bufsize = LOCALMERGE_BUFFERSIZE;
}
- for(c0 = ALPHABET_SIZE - 1, j = m; 0 < j; --c0) {
+#else
+ buf = SA + m, bufsize = n - (2 * m);
+ for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) {
for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) {
i = BUCKET_BSTAR(c0, c1);
if(1 < (j - i)) {
- substringsort(T, PAb, SA + i, SA + j,
- buf, bufsize, 2, *(SA + i) == (m - 1));
+ sssort(T, PAb, SA + i, SA + j,
+ buf, bufsize, 2, n, *(SA + i) == (m - 1));
}
}
}
- if(bufsize == LOCALMERGE_BUFFERSIZE) { free(buf); }
+#endif
/* Compute ranks of type B* substrings. */
for(i = m - 1; 0 <= i; --i) {
@@ -125,31 +159,30 @@ note:
trsort(ISAb, SA, m, 1);
/* Set the sorted order of tyoe B* suffixes. */
- for(i = n - 1, j = m; 0 <= i;) {
- for(--i; (0 <= i) && (T[i] >= T[i + 1]); --i) { }
+ for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) {
+ for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { }
if(0 <= i) {
- SA[ISAb[--j]] = i;
- for(--i; (0 <= i) && (T[i] <= T[i + 1]); --i) { }
+ t = i;
+ for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { }
+ SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t;
}
}
/* Calculate the index of start/end point of each bucket. */
- for(c0 = ALPHABET_SIZE - 1, i = n, k = m - 1; 0 <= c0; --c0) {
+ BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */
+ for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) {
+ i = BUCKET_A(c0 + 1) - 1;
for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) {
t = i - BUCKET_B(c0, c1);
- BUCKET_B(c0, c1) = i + 1; /* end point */
+ BUCKET_B(c0, c1) = i; /* end point */
/* Move all type B* suffixes to the correct position. */
for(i = t, j = BUCKET_BSTAR(c0, c1);
j <= k;
--i, --k) { SA[i] = SA[k]; }
}
- t = i - BUCKET_B(c0, c0);
- BUCKET_B(c0, c0) = i + 1; /* end point */
- if(c0 < (ALPHABET_SIZE - 1)) {
- BUCKET_BSTAR(c0, c0 + 1) = t + 1; /* start point */
- }
- i = BUCKET_A(c0);
+ BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */
+ BUCKET_B(c0, c0) = i; /* end point */
}
}
@@ -159,35 +192,37 @@ note:
/* Constructs the suffix array by using the sorted order of type B* suffixes. */
static
void
-_construct_SA(const sauchar_t *T, saidx_t *SA,
- saidx_t *bucket_A, saidx_t *bucket_B,
- saidx_t n, saidx_t m) {
- saidx_t *i, *j, *t;
+construct_SA(const sauchar_t *T, saidx_t *SA,
+ saidx_t *bucket_A, saidx_t *bucket_B,
+ saidx_t n, saidx_t m) {
+ saidx_t *i, *j, *k;
saidx_t s;
saint_t c0, c1, c2;
- /** An implementation version of MSufSort3's second stage. **/
-
if(0 < m) {
/* Construct the sorted order of type B suffixes by using
the sorted order of type B* suffixes. */
for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
/* Scan the suffix array from right to left. */
for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
- j = SA + BUCKET_A(c1 + 1), t = NULL, c2 = -1;
+ j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
i <= j;
--j) {
- if(0 <= (s = *j)) {
- if((0 <= --s) && ((c0 = T[s]) <= c1)) {
- *j = ~(s + 1);
- if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
- if(c2 == c0) { *--t = s; }
- else {
- if(0 <= c2) { BUCKET_B(c2, c1) = t - SA; }
- *(t = SA + BUCKET_B(c2 = c0, c1) - 1) = s;
- }
+ if(0 < (s = *j)) {
+ assert(T[s] == c1);
+ assert(((s + 1) < n) && (T[s] <= T[s + 1]));
+ assert(T[s - 1] <= T[s]);
+ *j = ~s;
+ c0 = T[--s];
+ if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
+ if(c0 != c2) {
+ if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
+ k = SA + BUCKET_B(c2 = c0, c1);
}
+ assert(k < j);
+ *k-- = s;
} else {
+ assert(((s == 0) && (T[s] == c1)) || (s < 0));
*j = ~s;
}
}
@@ -196,20 +231,22 @@ _construct_SA(const sauchar_t *T, saidx_t *SA,
/* Construct the suffix array by using
the sorted order of type B suffixes. */
- SA[0] = n;
- *(t = SA + BUCKET_A(c2 = T[n - 1]) + 1) = n - 1;
+ k = SA + BUCKET_A(c2 = T[n - 1]);
+ *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1);
/* Scan the suffix array from left to right. */
- for(i = SA + 1, j = SA + n; i <= j; ++i) {
- if(0 <= (s = *i)) {
- if((0 <= --s) && ((c0 = T[s]) >= T[s + 1])) {
- if((0 < s) && (T[s - 1] < c0)) { s = ~s; }
- if(c0 == c2) { *++t = s; }
- else {
- BUCKET_A(c2) = t - SA;
- *(t = SA + BUCKET_A(c2 = c0) + 1) = s;
- }
+ for(i = SA, j = SA + n; i < j; ++i) {
+ if(0 < (s = *i)) {
+ assert(T[s - 1] >= T[s]);
+ c0 = T[--s];
+ if((s == 0) || (T[s - 1] < c0)) { s = ~s; }
+ if(c0 != c2) {
+ BUCKET_A(c2) = k - SA;
+ k = SA + BUCKET_A(c2 = c0);
}
+ assert(i < k);
+ *k++ = s;
} else {
+ assert(s < 0);
*i = ~s;
}
}
@@ -219,36 +256,41 @@ _construct_SA(const sauchar_t *T, saidx_t *SA,
by using the sorted order of type B* suffixes. */
static
saidx_t
-_construct_BWT(const sauchar_t *T, saidx_t *SA,
- saidx_t *bucket_A, saidx_t *bucket_B,
- saidx_t n, saidx_t m) {
- saidx_t *i, *j, *t, *orig;
+construct_BWT(const sauchar_t *T, saidx_t *SA,
+ saidx_t *bucket_A, saidx_t *bucket_B,
+ saidx_t n, saidx_t m) {
+ saidx_t *i, *j, *k, *orig;
saidx_t s;
saint_t c0, c1, c2;
- /** An implementation version of MSufSort3's semidirect BWT. **/
-
if(0 < m) {
/* Construct the sorted order of type B suffixes by using
the sorted order of type B* suffixes. */
for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) {
/* Scan the suffix array from right to left. */
for(i = SA + BUCKET_BSTAR(c1, c1 + 1),
- j = SA + BUCKET_A(c1 + 1), t = NULL, c2 = -1;
+ j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1;
i <= j;
--j) {
- if(0 <= (s = *j)) {
- if((0 <= --s) && ((c0 = T[s]) <= c1)) {
- *j = ~((saidx_t)c0);
- if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
- if(c0 == c2) { *--t = s; }
- else {
- if(0 <= c2) { BUCKET_B(c2, c1) = t - SA; }
- *(t = SA + BUCKET_B(c2 = c0, c1) - 1) = s;
- }
+ if(0 < (s = *j)) {
+ assert(T[s] == c1);
+ assert(((s + 1) < n) && (T[s] <= T[s + 1]));
+ assert(T[s - 1] <= T[s]);
+ c0 = T[--s];
+ *j = ~((saidx_t)c0);
+ if((0 < s) && (T[s - 1] > c0)) { s = ~s; }
+ if(c0 != c2) {
+ if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; }
+ k = SA + BUCKET_B(c2 = c0, c1);
}
- } else {
+ assert(k < j);
+ *k-- = s;
+ } else if(s != 0) {
*j = ~s;
+#ifndef NDEBUG
+ } else {
+ assert(T[s] == c1);
+#endif
}
}
}
@@ -256,24 +298,25 @@ _construct_BWT(const sauchar_t *T, saidx_t *SA,
/* Construct the BWTed string by using
the sorted order of type B suffixes. */
- c0 = T[s = n - 1];
- SA[0] = c0;
- if(T[s - 1] < c0) { s = ~((saidx_t)T[s - 1]); }
- *(t = SA + BUCKET_A(c2 = c0) + 1) = s;
+ k = SA + BUCKET_A(c2 = T[n - 1]);
+ *k++ = (T[n - 2] < c2) ? ~((saidx_t)T[n - 2]) : (n - 1);
/* Scan the suffix array from left to right. */
- for(i = SA + 1, j = SA + n, orig = SA - 1; i <= j; ++i) {
- if(0 <= (s = *i)) {
- if((0 <= --s) && ((c0 = T[s]) >= T[s + 1])) {
- *i = c0;
- if((0 < s) && (T[s - 1] < c0)) { s = ~((saidx_t)T[s - 1]); }
- if(c0 == c2) { *++t = s; }
- else {
- BUCKET_A(c2) = t - SA;
- *(t = SA + BUCKET_A(c2 = c0) + 1) = s;
- }
- } else if(s < 0) { orig = i; }
- } else {
+ for(i = SA, j = SA + n, orig = SA; i < j; ++i) {
+ if(0 < (s = *i)) {
+ assert(T[s - 1] >= T[s]);
+ c0 = T[--s];
+ *i = c0;
+ if((0 < s) && (T[s - 1] < c0)) { s = ~((saidx_t)T[s - 1]); }
+ if(c0 != c2) {
+ BUCKET_A(c2) = k - SA;
+ k = SA + BUCKET_A(c2 = c0);
+ }
+ assert(i < k);
+ *k++ = s;
+ } else if(s != 0) {
*i = ~s;
+ } else {
+ orig = i;
}
}
@@ -293,20 +336,17 @@ divsufsort(const sauchar_t *T, saidx_t *SA, saidx_t n) {
/* Check arguments. */
if((T == NULL) || (SA == NULL) || (n < 0)) { return -1; }
- else if(n == 0) { SA[0] = 0; return 0; }
- else if(n == 1) { SA[0] = 1; SA[1] = 0; return 0; }
+ else if(n == 0) { return 0; }
+ else if(n == 1) { SA[0] = 0; return 0; }
+ else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; }
bucket_A = malloc(BUCKET_A_SIZE * sizeof(saidx_t));
bucket_B = malloc(BUCKET_B_SIZE * sizeof(saidx_t));
/* Suffixsort. */
if((bucket_A != NULL) && (bucket_B != NULL)) {
- m = _sort_typeBstar(T, SA, bucket_A, bucket_B, n);
- if(0 <= m) {
- _construct_SA(T, SA, bucket_A, bucket_B, n, m);
- } else {
- err = -3;
- }
+ m = sort_typeBstar(T, SA, bucket_A, bucket_B, n);
+ construct_SA(T, SA, bucket_A, bucket_B, n, m);
} else {
err = -2;
}
@@ -321,8 +361,7 @@ saidx_t
divbwt(const sauchar_t *T, sauchar_t *U, saidx_t *A, saidx_t n) {
saidx_t *B;
saidx_t *bucket_A, *bucket_B;
- saidx_t m, pidx = -1, i;
- saint_t err = 0;
+ saidx_t m, pidx, i;
/* Check arguments. */
if((T == NULL) || (U == NULL) || (n < 0)) { return -1; }
@@ -334,18 +373,16 @@ divbwt(const sauchar_t *T, sauchar_t *U, saidx_t *A, saidx_t n) {
/* Burrows-Wheeler Transform. */
if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) {
- m = _sort_typeBstar(T, B, bucket_A, bucket_B, n);
- if(0 <= m) {
- pidx = _construct_BWT(T, B, bucket_A, bucket_B, n, m);
-
- /* Copy to output string. */
- for(i = 0; i < pidx; ++i) { U[i] = B[i]; }
- for(; i < n; ++i) { U[i] = B[i + 1]; }
- } else {
- err = -3;
- }
+ m = sort_typeBstar(T, B, bucket_A, bucket_B, n);
+ pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m);
+
+ /* Copy to output string. */
+ U[0] = T[n - 1];
+ for(i = 0; i < pidx; ++i) { U[i + 1] = (sauchar_t)A[i]; }
+ for(i += 1; i < n; ++i) { U[i] = (sauchar_t)A[i]; }
+ pidx += 1;
} else {
- err = -2;
+ pidx = -2;
}
free(bucket_B);
@@ -355,7 +392,6 @@ divbwt(const sauchar_t *T, sauchar_t *U, saidx_t *A, saidx_t n) {
return pidx;
}
-
const char *
divsufsort_version(void) {
return PROJECT_VERSION_FULL;
diff --git a/lib/sssort.c b/lib/sssort.c
new file mode 100644
index 0000000..037ac97
--- /dev/null
+++ b/lib/sssort.c
@@ -0,0 +1,815 @@
+/*
+ * sssort.c for libdivsufsort
+ * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person
+ * obtaining a copy of this software and associated documentation
+ * files (the "Software"), to deal in the Software without
+ * restriction, including without limitation the rights to use,
+ * copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following
+ * conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+ * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+ * OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "divsufsort_private.h"
+
+
+/*- Private Functions -*/
+
+static const saint_t lg_table[256]= {
+ -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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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
+};
+
+#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
+
+static INLINE
+saint_t
+ss_ilg(saidx_t n) {
+#if SS_BLOCKSIZE == 0
+# if defined(BUILD_DIVSUFSORT64)
+ return (n & 0xffffffff00000000) ?
+ ((n & 0xffff000000000000) ?
+ ((n & 0xff00000000000000) ?
+ 56 + lg_table[(n >> 56) & 0xff] :
+ 48 + lg_table[(n >> 48) & 0xff]) :
+ ((n & 0x0000ff0000000000) ?
+ 40 + lg_table[(n >> 40) & 0xff] :
+ 32 + lg_table[(n >> 32) & 0xff])) :
+ ((n & 0x00000000ffff0000) ?
+ ((n & 0x00000000ff000000) ?
+ 24 + lg_table[(n >> 24) & 0xff] :
+ 16 + lg_table[(n >> 16) & 0xff]) :
+ ((n & 0x000000000000ff00) ?
+ 8 + lg_table[(n >> 8) & 0xff] :
+ 0 + lg_table[(n >> 0) & 0xff]));
+# else
+ return (n & 0xffff0000) ?
+ ((n & 0xff000000) ?
+ 24 + lg_table[(n >> 24) & 0xff] :
+ 16 + lg_table[(n >> 16) & 0xff]) :
+ ((n & 0x0000ff00) ?
+ 8 + lg_table[(n >> 8) & 0xff] :
+ 0 + lg_table[(n >> 0) & 0xff]);
+# endif
+#elif SS_BLOCKSIZE < 256
+ return lg_table[n];
+#else
+ return (n & 0xff00) ?
+ 8 + lg_table[(n >> 8) & 0xff] :
+ 0 + lg_table[(n >> 0) & 0xff];
+#endif
+}
+
+#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
+
+#if SS_BLOCKSIZE != 0
+
+static const int sqq_table[256] = {
+ 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61,
+ 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89,
+ 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109,
+110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
+128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
+143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155,
+156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
+169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180,
+181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191,
+192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201,
+202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211,
+212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221,
+221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230,
+230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
+239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247,
+247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255
+};
+
+static INLINE
+saidx_t
+ss_isqrt(saidx_t x) {
+ saidx_t y, e;
+
+ if(x >= (65535UL * 65535UL)) { return 65535; }
+ e = (x & 0xffff0000) ?
+ ((x & 0xff000000) ?
+ 24 + lg_table[(x >> 24) & 0xff] :
+ 16 + lg_table[(x >> 16) & 0xff]) :
+ ((x & 0x0000ff00) ?
+ 8 + lg_table[(x >> 8) & 0xff] :
+ 0 + lg_table[(x >> 0) & 0xff]);
+
+ if(e >= 16) {
+ y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7);
+ if(e >= 24) { y = (y + 1 + x / y) >> 1; }
+ y = (y + 1 + x / y) >> 1;
+ } else if(e >= 8) {
+ y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1;
+ } else {
+ return sqq_table[x] >> 4;
+ }
+
+ return (x < (y * y)) ? y - 1 : y;
+}
+
+#endif /* SS_BLOCKSIZE != 0 */
+
+
+/*---------------------------------------------------------------------------*/
+
+/* Compares two suffixes. */
+static INLINE
+saint_t
+ss_compare(const sauchar_t *T,
+ const saidx_t *p1, const saidx_t *p2,
+ saidx_t depth) {
+ const sauchar_t *U1, *U2, *U1n, *U2n;
+
+ for(U1 = T + depth + *p1,
+ U2 = T + depth + *p2,
+ U1n = T + *(p1 + 1) + 2,
+ U2n = T + *(p2 + 1) + 2;
+ (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
+ ++U1, ++U2) {
+ }
+
+ return U1 < U1n ?
+ (U2 < U2n ? *U1 - *U2 : 1) :
+ (U2 < U2n ? -1 : 0);
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1)
+
+/* Insertionsort for small size groups */
+static
+void
+ss_insertionsort(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *last, saidx_t depth) {
+ saidx_t *i, *j;
+ saidx_t t;
+ saint_t r;
+
+ for(i = last - 2; first <= i; --i) {
+ for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) {
+ do { *(j - 1) = *j; } while((++j < last) && (*j < 0));
+ if(last <= j) { break; }
+ }
+ if(r == 0) { *j = ~*j; }
+ *(j - 1) = t;
+ }
+}
+
+#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */
+
+
+/*---------------------------------------------------------------------------*/
+
+#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE)
+
+static INLINE
+void
+ss_fixdown(const sauchar_t *Td, const saidx_t *PA,
+ saidx_t *SA, saidx_t i, saidx_t size) {
+ saidx_t j, k;
+ saidx_t v;
+ saint_t c, d, e;
+
+ for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
+ d = Td[PA[SA[k = j++]]];
+ if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; }
+ if(d <= c) { break; }
+ }
+ SA[i] = v;
+}
+
+/* Simple top-down heapsort. */
+static
+void
+ss_heapsort(const sauchar_t *Td, const saidx_t *PA, saidx_t *SA, saidx_t size) {
+ saidx_t i, m;
+ saidx_t t;
+
+ m = size;
+ if((size % 2) == 0) {
+ m--;
+ if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); }
+ }
+
+ for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); }
+ if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); }
+ for(i = m - 1; 0 < i; --i) {
+ t = SA[0], SA[0] = SA[i];
+ ss_fixdown(Td, PA, SA, 0, i);
+ SA[i] = t;
+ }
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+/* Returns the median of three elements. */
+static INLINE
+saidx_t *
+ss_median3(const sauchar_t *Td, const saidx_t *PA,
+ saidx_t *v1, saidx_t *v2, saidx_t *v3) {
+ saidx_t *t;
+ if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); }
+ if(Td[PA[*v2]] > Td[PA[*v3]]) {
+ if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; }
+ else { return v3; }
+ }
+ return v2;
+}
+
+/* Returns the median of five elements. */
+static INLINE
+saidx_t *
+ss_median5(const sauchar_t *Td, const saidx_t *PA,
+ saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5) {
+ saidx_t *t;
+ if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); }
+ if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); }
+ if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); }
+ if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); }
+ if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); }
+ if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; }
+ return v3;
+}
+
+/* Returns the pivot element. */
+static INLINE
+saidx_t *
+ss_pivot(const sauchar_t *Td, const saidx_t *PA, saidx_t *first, saidx_t *last) {
+ saidx_t *middle;
+ saidx_t t;
+
+ t = last - first;
+ middle = first + t / 2;
+
+ if(t <= 512) {
+ if(t <= 32) {
+ return ss_median3(Td, PA, first, middle, last - 1);
+ } else {
+ t >>= 2;
+ return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1);
+ }
+ }
+ t >>= 3;
+ first = ss_median3(Td, PA, first, first + t, first + (t << 1));
+ middle = ss_median3(Td, PA, middle - t, middle, middle + t);
+ last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1);
+ return ss_median3(Td, PA, first, middle, last);
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+/* Binary partition for substrings. */
+static INLINE
+saidx_t *
+ss_partition(const saidx_t *PA,
+ saidx_t *first, saidx_t *last, saidx_t depth) {
+ saidx_t *a, *b;
+ saidx_t t;
+ for(a = first - 1, b = last;;) {
+ for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
+ for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { }
+ if(b <= a) { break; }
+ t = ~*b;
+ *b = *a;
+ *a = t;
+ }
+ if(first < a) { *first = ~*first; }
+ return a;
+}
+
+/* Multikey introsort for medium size groups. */
+static
+void
+ss_mintrosort(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *last,
+ saidx_t depth) {
+#define STACK_SIZE SS_MISORT_STACKSIZE
+ struct { saidx_t *a, *b, c; saint_t d; } stack[STACK_SIZE];
+ const sauchar_t *Td;
+ saidx_t *a, *b, *c, *d, *e, *f;
+ saidx_t s, t;
+ saint_t ssize;
+ saint_t limit;
+ saint_t v, x = 0;
+
+ for(ssize = 0, limit = ss_ilg(last - first);;) {
+
+ if((last - first) <= SS_INSERTIONSORT_THRESHOLD) {
+#if 1 < SS_INSERTIONSORT_THRESHOLD
+ if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); }
+#endif
+ STACK_POP(first, last, depth, limit);
+ continue;
+ }
+
+ Td = T + depth;
+ if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); }
+ if(limit < 0) {
+ for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) {
+ if((x = Td[PA[*a]]) != v) {
+ if(1 < (a - first)) { break; }
+ v = x;
+ first = a;
+ }
+ }
+ if(Td[PA[*first] - 1] < v) {
+ first = ss_partition(PA, first, a, depth);
+ }
+ if((a - first) <= (last - a)) {
+ if(1 < (a - first)) {
+ STACK_PUSH(a, last, depth, -1);
+ last = a, depth += 1, limit = ss_ilg(a - first);
+ } else {
+ first = a, limit = -1;
+ }
+ } else {
+ if(1 < (last - a)) {
+ STACK_PUSH(first, a, depth + 1, ss_ilg(a - first));
+ first = a, limit = -1;
+ } else {
+ last = a, depth += 1, limit = ss_ilg(a - first);
+ }
+ }
+ continue;
+ }
+
+ /* choose pivot */
+ a = ss_pivot(Td, PA, first, last);
+ v = Td[PA[*a]];
+ SWAP(*first, *a);
+
+ /* partition */
+ for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { }
+ if(((a = b) < last) && (x < v)) {
+ for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) {
+ if(x == v) { SWAP(*b, *a); ++a; }
+ }
+ }
+ for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { }
+ if((b < (d = c)) && (x > v)) {
+ for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
+ if(x == v) { SWAP(*c, *d); --d; }
+ }
+ }
+ for(; b < c;) {
+ SWAP(*b, *c);
+ for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) {
+ if(x == v) { SWAP(*b, *a); ++a; }
+ }
+ for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
+ if(x == v) { SWAP(*c, *d); --d; }
+ }
+ }
+
+ if(a <= d) {
+ c = b - 1;
+
+ if((s = a - first) > (t = b - a)) { s = t; }
+ for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
+ if((s = d - c) > (t = last - d - 1)) { s = t; }
+ for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
+
+ a = first + (b - a), c = last - (d - c);
+ b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth);
+
+ if((a - first) <= (last - c)) {
+ if((last - c) <= (c - b)) {
+ STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
+ STACK_PUSH(c, last, depth, limit);
+ last = a;
+ } else if((a - first) <= (c - b)) {
+ STACK_PUSH(c, last, depth, limit);
+ STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
+ last = a;
+ } else {
+ STACK_PUSH(c, last, depth, limit);
+ STACK_PUSH(first, a, depth, limit);
+ first = b, last = c, depth += 1, limit = ss_ilg(c - b);
+ }
+ } else {
+ if((a - first) <= (c - b)) {
+ STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
+ STACK_PUSH(first, a, depth, limit);
+ first = c;
+ } else if((last - c) <= (c - b)) {
+ STACK_PUSH(first, a, depth, limit);
+ STACK_PUSH(b, c, depth + 1, ss_ilg(c - b));
+ first = c;
+ } else {
+ STACK_PUSH(first, a, depth, limit);
+ STACK_PUSH(c, last, depth, limit);
+ first = b, last = c, depth += 1, limit = ss_ilg(c - b);
+ }
+ }
+ } else {
+ limit += 1;
+ if(Td[PA[*first] - 1] < v) {
+ first = ss_partition(PA, first, last, depth);
+ limit = ss_ilg(last - first);
+ }
+ depth += 1;
+ }
+ }
+#undef STACK_SIZE
+}
+
+#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */
+
+
+/*---------------------------------------------------------------------------*/
+
+#if SS_BLOCKSIZE != 0
+
+static INLINE
+void
+ss_blockswap(saidx_t *a, saidx_t *b, saidx_t n) {
+ saidx_t t;
+ for(; 0 < n; --n, ++a, ++b) {
+ t = *a, *a = *b, *b = t;
+ }
+}
+
+static INLINE
+void
+ss_rotate(saidx_t *first, saidx_t *middle, saidx_t *last) {
+ saidx_t *a, *b, t;
+ saidx_t l, r;
+ l = middle - first, r = last - middle;
+ for(; (0 < l) && (0 < r);) {
+ if(l == r) { ss_blockswap(first, middle, l); break; }
+ if(l < r) {
+ a = last - 1, b = middle - 1;
+ t = *a;
+ do {
+ *a-- = *b, *b-- = *a;
+ if(b < first) {
+ *a = t;
+ last = a;
+ if((r -= l + 1) <= l) { break; }
+ a -= 1, b = middle - 1;
+ t = *a;
+ }
+ } while(1);
+ } else {
+ a = first, b = middle;
+ t = *a;
+ do {
+ *a++ = *b, *b++ = *a;
+ if(last <= b) {
+ *a = t;
+ first = a + 1;
+ if((l -= r + 1) <= r) { break; }
+ a += 1, b = middle;
+ t = *a;
+ }
+ } while(1);
+ }
+ }
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+static
+void
+ss_inplacemerge(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *middle, saidx_t *last,
+ saidx_t depth) {
+ const saidx_t *p;
+ saidx_t *a, *b;
+ saidx_t len, half;
+ saint_t q, r;
+ saint_t x;
+
+ for(;;) {
+ if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); }
+ else { x = 0; p = PA + *(last - 1); }
+ for(a = first, len = middle - first, half = len >> 1, r = -1;
+ 0 < len;
+ len = half, half >>= 1) {
+ b = a + half;
+ q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth);
+ if(q < 0) {
+ a = b + 1;
+ half -= (len & 1) ^ 1;
+ } else {
+ r = q;
+ }
+ }
+ if(a < middle) {
+ if(r == 0) { *a = ~*a; }
+ ss_rotate(a, middle, last);
+ last -= middle - a;
+ middle = a;
+ if(first == middle) { break; }
+ }
+ --last;
+ if(x != 0) { while(*--last < 0) { } }
+ if(middle == last) { break; }
+ }
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+/* Merge-forward with internal buffer. */
+static
+void
+ss_mergeforward(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *middle, saidx_t *last,
+ saidx_t *buf, saidx_t depth) {
+ saidx_t *a, *b, *c, *bufend;
+ saidx_t t;
+ saint_t r;
+
+ bufend = buf + (middle - first) - 1;
+ ss_blockswap(buf, first, middle - first);
+
+ for(t = *(a = first), b = buf, c = middle;;) {
+ r = ss_compare(T, PA + *b, PA + *c, depth);
+ if(r < 0) {
+ do {
+ *a++ = *b;
+ if(bufend <= b) { *bufend = t; return; }
+ *b++ = *a;
+ } while(*b < 0);
+ } else if(r > 0) {
+ do {
+ *a++ = *c, *c++ = *a;
+ if(last <= c) {
+ while(b < bufend) { *a++ = *b, *b++ = *a; }
+ *a = *b, *b = t;
+ return;
+ }
+ } while(*c < 0);
+ } else {
+ *c = ~*c;
+ do {
+ *a++ = *b;
+ if(bufend <= b) { *bufend = t; return; }
+ *b++ = *a;
+ } while(*b < 0);
+
+ do {
+ *a++ = *c, *c++ = *a;
+ if(last <= c) {
+ while(b < bufend) { *a++ = *b, *b++ = *a; }
+ *a = *b, *b = t;
+ return;
+ }
+ } while(*c < 0);
+ }
+ }
+}
+
+/* Merge-backward with internal buffer. */
+static
+void
+ss_mergebackward(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *middle, saidx_t *last,
+ saidx_t *buf, saidx_t depth) {
+ const saidx_t *p1, *p2;
+ saidx_t *a, *b, *c, *bufend;
+ saidx_t t;
+ saint_t r;
+ saint_t x;
+
+ bufend = buf + (last - middle) - 1;
+ ss_blockswap(buf, middle, last - middle);
+
+ x = 0;
+ if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; }
+ else { p1 = PA + *bufend; }
+ if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; }
+ else { p2 = PA + *(middle - 1); }
+ for(t = *(a = last - 1), b = bufend, c = middle - 1;;) {
+ r = ss_compare(T, p1, p2, depth);
+ if(0 < r) {
+ if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
+ *a-- = *b;
+ if(b <= buf) { *buf = t; break; }
+ *b-- = *a;
+ if(*b < 0) { p1 = PA + ~*b; x |= 1; }
+ else { p1 = PA + *b; }
+ } else if(r < 0) {
+ if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
+ *a-- = *c, *c-- = *a;
+ if(c < first) {
+ while(buf < b) { *a-- = *b, *b-- = *a; }
+ *a = *b, *b = t;
+ break;
+ }
+ if(*c < 0) { p2 = PA + ~*c; x |= 2; }
+ else { p2 = PA + *c; }
+ } else {
+ if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; }
+ *a-- = ~*b;
+ if(b <= buf) { *buf = t; break; }
+ *b-- = *a;
+ if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; }
+ *a-- = *c, *c-- = *a;
+ if(c < first) {
+ while(buf < b) { *a-- = *b, *b-- = *a; }
+ *a = *b, *b = t;
+ break;
+ }
+ if(*b < 0) { p1 = PA + ~*b; x |= 1; }
+ else { p1 = PA + *b; }
+ if(*c < 0) { p2 = PA + ~*c; x |= 2; }
+ else { p2 = PA + *c; }
+ }
+ }
+}
+
+/* D&C based merge. */
+static
+void
+ss_swapmerge(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *middle, saidx_t *last,
+ saidx_t *buf, saidx_t bufsize, saidx_t depth) {
+#define STACK_SIZE SS_SMERGE_STACKSIZE
+#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
+#define MERGE_CHECK(a, b, c)\
+ do {\
+ if(((c) & 1) ||\
+ (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\
+ *(a) = ~*(a);\
+ }\
+ if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\
+ *(b) = ~*(b);\
+ }\
+ } while(0)
+ struct { saidx_t *a, *b, *c; saint_t d; } stack[STACK_SIZE];
+ saidx_t *l, *r, *lm, *rm;
+ saidx_t m, len, half;
+ saint_t ssize;
+ saint_t check, next;
+
+ for(check = 0, ssize = 0;;) {
+ if((last - middle) <= bufsize) {
+ if((first < middle) && (middle < last)) {
+ ss_mergebackward(T, PA, first, middle, last, buf, depth);
+ }
+ MERGE_CHECK(first, last, check);
+ STACK_POP(first, middle, last, check);
+ continue;
+ }
+
+ if((middle - first) <= bufsize) {
+ if(first < middle) {
+ ss_mergeforward(T, PA, first, middle, last, buf, depth);
+ }
+ MERGE_CHECK(first, last, check);
+ STACK_POP(first, middle, last, check);
+ continue;
+ }
+
+ for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1;
+ 0 < len;
+ len = half, half >>= 1) {
+ if(ss_compare(T, PA + GETIDX(*(middle + m + half)),
+ PA + GETIDX(*(middle - m - half - 1)), depth) < 0) {
+ m += half + 1;
+ half -= (len & 1) ^ 1;
+ }
+ }
+
+ if(0 < m) {
+ lm = middle - m, rm = middle + m;
+ ss_blockswap(lm, middle, m);
+ l = r = middle, next = 0;
+ if(rm < last) {
+ if(*rm < 0) {
+ *rm = ~*rm;
+ if(first < lm) { for(; *--l < 0;) { } next |= 4; }
+ next |= 1;
+ } else if(first < lm) {
+ for(; *r < 0; ++r) { }
+ next |= 2;
+ }
+ }
+
+ if((l - first) <= (last - r)) {
+ STACK_PUSH(r, rm, last, (next & 3) | (check & 4));
+ middle = lm, last = l, check = (check & 3) | (next & 4);
+ } else {
+ if((next & 2) && (r == middle)) { next ^= 6; }
+ STACK_PUSH(first, lm, l, (check & 3) | (next & 4));
+ first = r, middle = rm, check = (next & 3) | (check & 4);
+ }
+ } else {
+ if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) {
+ *middle = ~*middle;
+ }
+ MERGE_CHECK(first, last, check);
+ STACK_POP(first, middle, last, check);
+ }
+ }
+#undef STACK_SIZE
+}
+
+#endif /* SS_BLOCKSIZE != 0 */
+
+
+/*---------------------------------------------------------------------------*/
+
+/*- Function -*/
+
+/* Substring sort */
+void
+sssort(const sauchar_t *T, const saidx_t *PA,
+ saidx_t *first, saidx_t *last,
+ saidx_t *buf, saidx_t bufsize,
+ saidx_t depth, saidx_t n, saint_t lastsuffix) {
+ saidx_t *a;
+#if SS_BLOCKSIZE != 0
+ saidx_t *b, *middle, *curbuf;
+ saidx_t j, k, curbufsize, limit;
+#endif
+ saidx_t i;
+
+ if(lastsuffix != 0) { ++first; }
+
+#if SS_BLOCKSIZE == 0
+ ss_mintrosort(T, PA, first, last, depth);
+#else
+ if((bufsize < SS_BLOCKSIZE) &&
+ (bufsize < (last - first)) &&
+ (bufsize < (limit = ss_isqrt(last - first)))) {
+ if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; }
+ buf = middle = last - limit, bufsize = limit;
+ } else {
+ middle = last, limit = 0;
+ }
+ for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) {
+#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
+ ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth);
+#elif 1 < SS_BLOCKSIZE
+ ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth);
+#endif
+ curbufsize = last - (a + SS_BLOCKSIZE);
+ curbuf = a + SS_BLOCKSIZE;
+ if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
+ for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) {
+ ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
+ }
+ }
+#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
+ ss_mintrosort(T, PA, a, middle, depth);
+#elif 1 < SS_BLOCKSIZE
+ ss_insertionsort(T, PA, a, middle, depth);
+#endif
+ for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) {
+ if(i & 1) {
+ ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth);
+ a -= k;
+ }
+ }
+ if(limit != 0) {
+#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE
+ ss_mintrosort(T, PA, middle, last, depth);
+#elif 1 < SS_BLOCKSIZE
+ ss_insertionsort(T, PA, middle, last, depth);
+#endif
+ ss_inplacemerge(T, PA, first, middle, last, depth);
+ }
+#endif
+
+ if(lastsuffix != 0) {
+ /* Insert last type B* suffix. */
+ saidx_t PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2;
+ for(a = first, i = *(first - 1);
+ (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth)));
+ ++a) {
+ *(a - 1) = *a;
+ }
+ *(a - 1) = i;
+ }
+}
diff --git a/lib/substringsort.c b/lib/substringsort.c
deleted file mode 100644
index d97ebb0..0000000
--- a/lib/substringsort.c
+++ /dev/null
@@ -1,619 +0,0 @@
-/*
- * substringsort.c for libdivsufsort
- * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved.
- *
- * Permission is hereby granted, free of charge, to any person
- * obtaining a copy of this software and associated documentation
- * files (the "Software"), to deal in the Software without
- * restriction, including without limitation the rights to use,
- * copy, modify, merge, publish, distribute, sublicense, and/or sell
- * copies of the Software, and to permit persons to whom the
- * Software is furnished to do so, subject to the following
- * conditions:
- *
- * The above copyright notice and this permission notice shall be
- * included in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
- * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
- * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
- * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
- * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
- * OTHER DEALINGS IN THE SOFTWARE.
- */
-
-#include "divsufsort_private.h"
-
-
-/*- Private Functions -*/
-
-/* Compares two suffixes. */
-static INLINE
-saint_t
-_compare(const sauchar_t *T,
- const saidx_t *p1, const saidx_t *p2,
- saidx_t depth) {
- const sauchar_t *U1, *U2, *U1n, *U2n;
-
- for(U1 = T + depth + *p1,
- U2 = T + depth + *p2,
- U1n = T + *(p1 + 1) + 2,
- U2n = T + *(p2 + 1) + 2;
- (U1 < U1n) && (U2 < U2n) && (*U1 == *U2);
- ++U1, ++U2) {
- }
-
- return U1 < U1n ?
- (U2 < U2n ? *U1 - *U2 : 1) :
- (U2 < U2n ? -1 : 0);
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-/* Insertionsort for small size groups */
-static
-void
-_insertionsort(const sauchar_t *T, const saidx_t *PA,
- saidx_t *first, saidx_t *last, saidx_t depth) {
- saidx_t *i, *j;
- saidx_t t;
- saint_t r;
-
- for(i = last - 2; first <= i; --i) {
- for(t = *i, j = i + 1; 0 < (r = _compare(T, PA + t, PA + *j, depth));) {
- do { *(j - 1) = *j; } while((++j < last) && (*j < 0));
- if(last <= j) { break; }
- }
- if(r == 0) { *j = ~*j; }
- *(j - 1) = t;
- }
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-static INLINE
-void
-_fixdown(const sauchar_t *Td, const saidx_t *PA,
- saidx_t *SA, saidx_t i, saidx_t size) {
- saidx_t j, k;
- saidx_t v;
- saint_t c, d, e;
-
- for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) {
- d = Td[PA[SA[k = j++]]];
- if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; }
- if(d <= c) { break; }
- }
- SA[i] = v;
-}
-
-/* Simple top-down heapsort. */
-static
-void
-_heapsort(const sauchar_t *Td, const saidx_t *PA, saidx_t *SA, saidx_t size) {
- saidx_t i, m;
- saidx_t t;
-
- m = size;
- if((size % 2) == 0) {
- m--;
- if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); }
- }
-
- for(i = m / 2 - 1; 0 <= i; --i) { _fixdown(Td, PA, SA, i, m); }
-
- if((size % 2) == 0) {
- SWAP(SA[0], SA[m]);
- _fixdown(Td, PA, SA, 0, m);
- }
-
- for(i = m - 1; 0 < i; --i) {
- t = SA[0];
- SA[0] = SA[i];
- _fixdown(Td, PA, SA, 0, i);
- SA[i] = t;
- }
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-/* Returns the median of three elements. */
-static INLINE
-saidx_t *
-_median3(const sauchar_t *Td, const saidx_t *PA,
- saidx_t *v1, saidx_t *v2, saidx_t *v3) {
- saidx_t *t;
- if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); }
- if(Td[PA[*v2]] > Td[PA[*v3]]) {
- if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; }
- else { return v3; }
- }
- return v2;
-}
-
-/* Returns the median of five elements. */
-static INLINE
-saidx_t *
-_median5(const sauchar_t *Td, const saidx_t *PA,
- saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5) {
- saidx_t *t;
- if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); }
- if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); }
- if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); }
- if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); }
- if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); }
- if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; }
- return v3;
-}
-
-/* Returns the pivot element. */
-static INLINE
-saidx_t *
-_pivot(const sauchar_t *Td, const saidx_t *PA, saidx_t *first, saidx_t *last) {
- saidx_t *middle;
- saidx_t t;
-
- t = last - first;
- middle = first + t / 2;
-
- if(t <= 512) {
- if(t <= 32) {
- return _median3(Td, PA, first, middle, last - 1);
- } else {
- t >>= 2;
- return _median5(Td, PA,
- first, first + t, middle,
- last - 1 - t, last - 1);
- }
- }
- t >>= 3;
- return _median3(Td, PA,
- _median3(Td, PA, first, first + t, first + (t << 1)),
- _median3(Td, PA, middle - t, middle, middle + t),
- _median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1));
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-static INLINE
-saidx_t
-_lg(saidx_t n) {
-static const int log2table[256]= {
- -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,
- 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,
- 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,
- 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,
- 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,
- 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,
- 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,
- 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
-};
- /* for 16 bits */
- return ((n & 0xff00) != 0) ?
- 8 + log2table[(n >> 8) & 0xff] :
- log2table[n & 0xff];
-
- /* for 32 bits */
-/*
- return (n & 0xffff0000) ?
- ((n & 0xff000000) ?
- 24 + log2table[(n >> 24) & 0xff] :
- 16 + log2table[(n >> 16) & 0xff]) :
- ((n & 0x0000ff00) ?
- 8 + log2table[(n >> 8) & 0xff] :
- 0 + log2table[(n >> 0) & 0xff]);
-*/
-}
-
-/* Binary partition for substrings. */
-static INLINE
-saidx_t *
-_substring_partition(const saidx_t *PA,
- saidx_t *first, saidx_t *last, saidx_t depth) {
- saidx_t *a, *b;
- saidx_t t;
- for(a = first - 1, b = last;;) {
- for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; }
- for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { }
- if(b <= a) { break; }
- t = ~*b;
- *b = *a;
- *a = t;
- }
- if(first < a) { *first = ~*first; }
- return a;
-}
-
-/* Multikey introsort for medium size groups. */
-static
-void
-_multikey_introsort(const sauchar_t *T, const saidx_t *PA,
- saidx_t *first, saidx_t *last,
- saidx_t depth) {
- struct { saidx_t *a, *b, c, d; } stack[STACK_SIZE];
- const sauchar_t *Td;
- saidx_t *a, *b, *c, *d, *e, *f;
- saidx_t s, t;
- saidx_t ssize;
- saidx_t limit;
- saint_t v, x = 0;
-
- for(ssize = 0, limit = _lg(last - first);;) {
-
- if((last - first) <= SS_INSERTIONSORT_THRESHOLD) {
- if(1 < (last - first)) { _insertionsort(T, PA, first, last, depth); }
- STACK_POP(first, last, depth, limit);
- continue;
- }
-
- Td = T + depth;
- if(limit-- == 0) { _heapsort(Td, PA, first, last - first); }
- if(limit < 0) {
- for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) {
- if((x = Td[PA[*a]]) != v) {
- if(1 < (a - first)) { break; }
- v = x;
- first = a;
- }
- }
- if(Td[PA[*first] - 1] < v) {
- first = _substring_partition(PA, first, a, depth);
- }
- if((a - first) <= (last - a)) {
- if(1 < (a - first)) {
- STACK_PUSH(a, last, depth, -1);
- last = a, depth += 1, limit = _lg(a - first);
- } else {
- first = a, limit = -1;
- }
- } else {
- if(1 < (last - a)) {
- STACK_PUSH(first, a, depth + 1, _lg(a - first));
- first = a, limit = -1;
- } else {
- last = a, depth += 1, limit = _lg(a - first);
- }
- }
- continue;
- }
-
- /* choose pivot */
- a = _pivot(Td, PA, first, last);
- v = Td[PA[*a]];
- SWAP(*first, *a);
-
- /* partition */
- for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { }
- if(((a = b) < last) && (x < v)) {
- for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- }
- for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { }
- if((b < (d = c)) && (x > v)) {
- for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
- for(; b < c;) {
- SWAP(*b, *c);
- for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
-
- if(a <= d) {
- c = b - 1;
-
- if((s = a - first) > (t = b - a)) { s = t; }
- for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
- if((s = d - c) > (t = last - d - 1)) { s = t; }
- for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
-
- a = first + (b - a), c = last - (d - c);
- b = (v <= Td[PA[*a] - 1]) ? a : _substring_partition(PA, a, c, depth);
-
- if((a - first) <= (last - c)) {
- if((last - c) <= (c - b)) {
- STACK_PUSH(b, c, depth + 1, _lg(c - b));
- STACK_PUSH(c, last, depth, limit);
- last = a;
- } else if((a - first) <= (c - b)) {
- STACK_PUSH(c, last, depth, limit);
- STACK_PUSH(b, c, depth + 1, _lg(c - b));
- last = a;
- } else {
- STACK_PUSH(c, last, depth, limit);
- STACK_PUSH(first, a, depth, limit);
- first = b, last = c, depth += 1, limit = _lg(c - b);
- }
- } else {
- if((a - first) <= (c - b)) {
- STACK_PUSH(b, c, depth + 1, _lg(c - b));
- STACK_PUSH(first, a, depth, limit);
- first = c;
- } else if((last - c) <= (c - b)) {
- STACK_PUSH(first, a, depth, limit);
- STACK_PUSH(b, c, depth + 1, _lg(c - b));
- first = c;
- } else {
- STACK_PUSH(first, a, depth, limit);
- STACK_PUSH(c, last, depth, limit);
- first = b, last = c, depth += 1, limit = _lg(c - b);
- }
- }
- } else {
- limit += 1;
- if(Td[PA[*first] - 1] < v) {
- first = _substring_partition(PA, first, last, depth);
- limit = _lg(last - first);
- }
- depth += 1;
- }
- }
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-/* Block swapping */
-static INLINE
-void
-_block_swap(saidx_t *first1, saidx_t *first2, saidx_t size) {
- saidx_t *a, *b;
- saidx_t i, t;
- for(i = size, a = first1, b = first2; 0 < i; --i, ++a, ++b) {
- SWAP(*a, *b);
- }
-}
-
-/* Merge-forward with internal buffer. */
-static
-void
-_merge_forward(const sauchar_t *T, const saidx_t *PA,
- saidx_t *buf, saidx_t *first, saidx_t *middle, saidx_t *last,
- saidx_t depth) {
- saidx_t *bufend;
- saidx_t *i, *j, *k;
- saidx_t t;
- saint_t r;
-
- bufend = buf + (middle - first) - 1;
- _block_swap(buf, first, middle - first);
-
- for(t = *first, i = first, j = buf, k = middle;;) {
- r = _compare(T, PA + *j, PA + *k, depth);
- if(r < 0) {
- do {
- *i++ = *j;
- if(bufend <= j) { *j = t; return; }
- *j++ = *i;
- } while(*j < 0);
- } else if(r > 0) {
- do {
- *i++ = *k, *k++ = *i;
- if(last <= k) {
- while(j < bufend) { *i++ = *j, *j++ = *i; }
- *i = *j, *j = t;
- return;
- }
- } while(*k < 0);
- } else {
- *k = ~*k;
- do {
- *i++ = *j;
- if(bufend <= j) { *j = t; return; }
- *j++ = *i;
- } while(*j < 0);
-
- do {
- *i++ = *k, *k++ = *i;
- if(last <= k) {
- while(j < bufend) { *i++ = *j, *j++ = *i; }
- *i = *j, *j = t;
- return;
- }
- } while(*k < 0);
- }
- }
-}
-
-/* Merge-backward with internal buffer. */
-static
-void
-_merge_backward(const sauchar_t *T, const saidx_t *PA, saidx_t *buf,
- saidx_t *first, saidx_t *middle, saidx_t *last,
- saidx_t depth) {
- const saidx_t *p1, *p2;
- saidx_t *bufend;
- saidx_t *i, *j, *k;
- saidx_t t;
- saint_t r;
- saint_t x;
-
- bufend = buf + (last - middle);
- _block_swap(buf, middle, last - middle);
-
- x = 0;
- if(*(bufend - 1) < 0) { x |= 1; p1 = PA + ~*(bufend - 1); }
- else { p1 = PA + *(bufend - 1); }
- if(*(middle - 1) < 0) { x |= 2; p2 = PA + ~*(middle - 1); }
- else { p2 = PA + *(middle - 1); }
- for(t = *(last - 1), i = last - 1, j = bufend - 1, k = middle - 1;;) {
-
- r = _compare(T, p1, p2, depth);
- if(r > 0) {
- if(x & 1) { do { *i-- = *j; *j-- = *i; } while(*j < 0); x ^= 1; }
- *i-- = *j;
- if(j <= buf) { *j = t; return; }
- *j-- = *i;
-
- if(*j < 0) { x |= 1; p1 = PA + ~*j; }
- else { p1 = PA + *j; }
- } else if(r < 0) {
- if(x & 2) { do { *i-- = *k; *k-- = *i; } while(*k < 0); x ^= 2; }
- *i-- = *k; *k-- = *i;
- if(k < first) {
- while(buf < j) { *i-- = *j; *j-- = *i; }
- *i = *j, *j = t;
- return;
- }
-
- if(*k < 0) { x |= 2; p2 = PA + ~*k; }
- else { p2 = PA + *k; }
- } else {
- if(x & 1) { do { *i-- = *j; *j-- = *i; } while(*j < 0); x ^= 1; }
- *i-- = ~*j;
- if(j <= buf) { *j = t; return; }
- *j-- = *i;
-
- if(x & 2) { do { *i-- = *k; *k-- = *i; } while(*k < 0); x ^= 2; }
- *i-- = *k; *k-- = *i;
- if(k < first) {
- while(buf < j) { *i-- = *j; *j-- = *i; }
- *i = *j, *j = t;
- return;
- }
-
- if(*j < 0) { x |= 1; p1 = PA + ~*j; }
- else { p1 = PA + *j; }
- if(*k < 0) { x |= 2; p2 = PA + ~*k; }
- else { p2 = PA + *k; }
- }
- }
-}
-
-/* Faster merge (based on divide and conquer technique). */
-static
-void
-_merge(const sauchar_t *T, const saidx_t *PA,
- saidx_t *first, saidx_t *middle, saidx_t *last,
- saidx_t *buf, saidx_t bufsize,
- saidx_t depth) {
-#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a)))
-#define MERGE_CHECK_EQUAL(a)\
- do {\
- if((0 <= *(a)) &&\
- (_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0)) {\
- *(a) = ~*(a);\
- }\
- } while(0)
- struct { saidx_t *a, *b, *c; int d; } stack[STACK_SIZE];
- saidx_t *i, *j;
- saidx_t m, len, half;
- saidx_t ssize;
- int check, next;
-
- for(check = 0, ssize = 0;;) {
-
- if((last - middle) <= bufsize) {
- if((first < middle) && (middle < last)) {
- _merge_backward(T, PA, buf, first, middle, last, depth);
- }
- if(check & 1) { MERGE_CHECK_EQUAL(first); }
- if(check & 2) { MERGE_CHECK_EQUAL(last); }
- STACK_POP(first, middle, last, check);
- continue;
- }
-
- if((middle - first) <= bufsize) {
- if(first < middle) {
- _merge_forward(T, PA, buf, first, middle, last, depth);
- }
- if(check & 1) { MERGE_CHECK_EQUAL(first); }
- if(check & 2) { MERGE_CHECK_EQUAL(last); }
- STACK_POP(first, middle, last, check);
- continue;
- }
-
- for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1;
- 0 < len;
- len = half, half >>= 1) {
- if(_compare(T, PA + GETIDX(*(middle + m + half)),
- PA + GETIDX(*(middle - m - half - 1)), depth) < 0) {
- m += half + 1;
- half -= (len & 1) ^ 1;
- }
- }
-
- if(0 < m) {
- _block_swap(middle - m, middle, m);
- i = j = middle, next = 0;
- if((middle + m) < last) {
- if(*(middle + m) < 0) {
- for(; *(i - 1) < 0; --i) { }
- *(middle + m) = ~*(middle + m);
- }
- for(j = middle; *j < 0; ++j) { }
- next = 1;
- }
- if((i - first) <= (last - j)) {
- STACK_PUSH(j, middle + m, last, (check & 2) | (next & 1));
- middle -= m, last = i, check = (check & 1);
- } else {
- if((i == middle) && (middle == j)) { next <<= 1; }
- STACK_PUSH(first, middle - m, i, (check & 1) | (next & 2));
- first = j, middle += m, check = (check & 2) | (next & 1);
- }
- } else {
- if(check & 1) { MERGE_CHECK_EQUAL(first); }
- MERGE_CHECK_EQUAL(middle);
- if(check & 2) { MERGE_CHECK_EQUAL(last); }
- STACK_POP(first, middle, last, check);
- }
- }
-}
-
-
-/*---------------------------------------------------------------------------*/
-
-/*- Function -*/
-
-/* Substring sort */
-void
-substringsort(const sauchar_t *T, const saidx_t *PA,
- saidx_t *first, saidx_t *last,
- saidx_t *buf, saidx_t bufsize,
- saidx_t depth, saint_t lastsuffix) {
- saidx_t *a, *b;
- saidx_t *curbuf;
- saidx_t i, j, k;
- saidx_t curbufsize;
-
- if(lastsuffix != 0) { ++first; }
- for(a = first, i = 0; (a + SS_BLOCKSIZE) < last; a += SS_BLOCKSIZE, ++i) {
- _multikey_introsort(T, PA, a, a + SS_BLOCKSIZE, depth);
- curbuf = a + SS_BLOCKSIZE;
- curbufsize = last - (a + SS_BLOCKSIZE);
- if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; }
- for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) {
- _merge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth);
- }
- }
- _multikey_introsort(T, PA, a, last, depth);
- for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) {
- if(i & 1) {
- _merge(T, PA, a - k, a, last, buf, bufsize, depth);
- a -= k;
- }
- }
-
- if(lastsuffix != 0) {
- /* Insert last type B* suffix. */
- for(a = first, i = *(first - 1);
- (a < last) && ((*a < 0) || (0 < _compare(T, PA + i, PA + *a, depth)));
- ++a) {
- *(a - 1) = *a;
- }
- *(a - 1) = i;
- }
-}
diff --git a/lib/trsort.c b/lib/trsort.c
index 6f646f1..4ec6a17 100644
--- a/lib/trsort.c
+++ b/lib/trsort.c
@@ -29,9 +29,73 @@
/*- Private Functions -*/
+static const saint_t lg_table[256]= {
+ -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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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,
+ 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
+};
+
static INLINE
+saint_t
+tr_ilg(saidx_t n) {
+#if defined(BUILD_DIVSUFSORT64)
+ return (n & 0xffffffff00000000) ?
+ ((n & 0xffff000000000000) ?
+ ((n & 0xff00000000000000) ?
+ 56 + lg_table[(n >> 56) & 0xff] :
+ 48 + lg_table[(n >> 48) & 0xff]) :
+ ((n & 0x0000ff0000000000) ?
+ 40 + lg_table[(n >> 40) & 0xff] :
+ 32 + lg_table[(n >> 32) & 0xff])) :
+ ((n & 0x00000000ffff0000) ?
+ ((n & 0x00000000ff000000) ?
+ 24 + lg_table[(n >> 24) & 0xff] :
+ 16 + lg_table[(n >> 16) & 0xff]) :
+ ((n & 0x000000000000ff00) ?
+ 8 + lg_table[(n >> 8) & 0xff] :
+ 0 + lg_table[(n >> 0) & 0xff]));
+#else
+ return (n & 0xffff0000) ?
+ ((n & 0xff000000) ?
+ 24 + lg_table[(n >> 24) & 0xff] :
+ 16 + lg_table[(n >> 16) & 0xff]) :
+ ((n & 0x0000ff00) ?
+ 8 + lg_table[(n >> 8) & 0xff] :
+ 0 + lg_table[(n >> 0) & 0xff]);
+#endif
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+/* Simple insertionsort for small size groups. */
+static
void
-_fixdown(const saidx_t *ISAd, saidx_t *SA, saidx_t i, saidx_t size) {
+tr_insertionsort(const saidx_t *ISAd, saidx_t *first, saidx_t *last) {
+ saidx_t *a, *b;
+ saidx_t t, r;
+
+ for(a = first + 1; a < last; ++a) {
+ for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) {
+ do { *(b + 1) = *b; } while((first <= --b) && (*b < 0));
+ if(b < first) { break; }
+ }
+ if(r == 0) { *b = ~*b; }
+ *(b + 1) = t;
+ }
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+static INLINE
+void
+tr_fixdown(const saidx_t *ISAd, saidx_t *SA, saidx_t i, saidx_t size) {
saidx_t j, k;
saidx_t v;
saidx_t c, d, e;
@@ -47,88 +111,32 @@ _fixdown(const saidx_t *ISAd, saidx_t *SA, saidx_t i, saidx_t size) {
/* Simple top-down heapsort. */
static
void
-_heapsort(const saidx_t *ISAd, saidx_t *SA, saidx_t size) {
+tr_heapsort(const saidx_t *ISAd, saidx_t *SA, saidx_t size) {
saidx_t i, m;
saidx_t t;
m = size;
if((size % 2) == 0) {
m--;
- if(ISAd[SA[m / 2]] < ISAd[SA[m]]) {
- SWAP(SA[m], SA[m / 2]);
- }
- }
-
- for(i = m / 2 - 1; 0 <= i; --i) {
- _fixdown(ISAd, SA, i, m);
- }
-
- if((size % 2) == 0) {
- SWAP(SA[0], SA[m]);
- _fixdown(ISAd, SA, 0, m);
+ if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); }
}
+ for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); }
+ if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); }
for(i = m - 1; 0 < i; --i) {
- t = SA[0];
- SA[0] = SA[i];
- _fixdown(ISAd, SA, 0, i);
+ t = SA[0], SA[0] = SA[i];
+ tr_fixdown(ISAd, SA, 0, i);
SA[i] = t;
}
}
-/* Simple insertionsort for small size groups. */
-static
-void
-_insertionsort(const saidx_t *ISAd, saidx_t *first, saidx_t *last) {
- saidx_t *a, *b;
- saidx_t t, r;
-
- for(a = first + 1; a < last; ++a) {
- for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) {
- do { *(b + 1) = *b; } while((first <= --b) && (*b < 0));
- if(b < first) { break; }
- }
- if(r == 0) { *b = ~*b; }
- *(b + 1) = t;
- }
-}
-
-static INLINE
-saidx_t
-_lg(saidx_t n) {
-static const int log2table[256]= {
- -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,
- 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,
- 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,
- 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,
- 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,
- 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,
- 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,
- 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
-};
- /* for 16 bits */
-/*
- return ((n & 0xff00) != 0) ?
- 8 + log2table[(n >> 8) & 0xff] :
- log2table[n & 0xff];
-*/
- /* for 32 bits */
- return (n & 0xffff0000) ?
- ((n & 0xff000000) ?
- 24 + log2table[(n >> 24) & 0xff] :
- 16 + log2table[(n >> 16) & 0xff]) :
- ((n & 0x0000ff00) ?
- 8 + log2table[(n >> 8) & 0xff] :
- 0 + log2table[(n >> 0) & 0xff]);
-}
-
/*---------------------------------------------------------------------------*/
/* Returns the median of three elements. */
static INLINE
saidx_t *
-_median3(const saidx_t *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3) {
+tr_median3(const saidx_t *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3) {
saidx_t *t;
if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); }
if(ISAd[*v2] > ISAd[*v3]) {
@@ -141,8 +149,8 @@ _median3(const saidx_t *ISAd, saidx_t *v1, saidx_t *v2, saidx_t *v3) {
/* Returns the median of five elements. */
static INLINE
saidx_t *
-_median5(const saidx_t *ISAd,
- saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5) {
+tr_median5(const saidx_t *ISAd,
+ saidx_t *v1, saidx_t *v2, saidx_t *v3, saidx_t *v4, saidx_t *v5) {
saidx_t *t;
if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); }
if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); }
@@ -156,7 +164,7 @@ _median5(const saidx_t *ISAd,
/* Returns the pivot element. */
static INLINE
saidx_t *
-_pivot(const saidx_t *ISAd, saidx_t *first, saidx_t *last) {
+tr_pivot(const saidx_t *ISAd, saidx_t *first, saidx_t *last) {
saidx_t *middle;
saidx_t t;
@@ -165,184 +173,60 @@ _pivot(const saidx_t *ISAd, saidx_t *first, saidx_t *last) {
if(t <= 512) {
if(t <= 32) {
- return _median3(ISAd, first, middle, last - 1);
+ return tr_median3(ISAd, first, middle, last - 1);
} else {
t >>= 2;
- return _median5(ISAd,
- first, first + t,
- middle,
- last - 1 - t, last - 1);
+ return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1);
}
}
t >>= 3;
- return _median3(ISAd,
- _median3(ISAd, first, first + t, first + (t << 1)),
- _median3(ISAd, middle - t, middle, middle + t),
- _median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1));
+ first = tr_median3(ISAd, first, first + t, first + (t << 1));
+ middle = tr_median3(ISAd, middle - t, middle, middle + t);
+ last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1);
+ return tr_median3(ISAd, first, middle, last);
}
/*---------------------------------------------------------------------------*/
-/* Update group numbers. */
-static
-void
-_ls_updategroup(saidx_t *ISA, const saidx_t *SA,
- saidx_t *first, saidx_t *last) {
- saidx_t *a, *b;
- saidx_t t;
-
- for(a = first; a < last; ++a) {
- if(0 <= *a) {
- b = a;
- do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a));
- *b = b - a;
- if(last <= a) { break; }
- }
- b = a;
- do { *a = ~*a; } while(*++a < 0);
- t = a - SA;
- do { ISA[*b] = t; } while(++b <= a);
- }
-}
+typedef struct _trbudget_t trbudget_t;
+struct _trbudget_t {
+ saidx_t chance;
+ saidx_t remain;
+ saidx_t incval;
+ saidx_t count;
+};
-/* Introspective sort. */
-static
+static INLINE
void
-_ls_introsort(saidx_t *ISA, saidx_t *ISAd, const saidx_t *SA,
- saidx_t *first, saidx_t *last) {
- struct { saidx_t *a, *b, c; } stack[STACK_SIZE];
- saidx_t *a, *b, *c, *d, *e, *f;
- saidx_t s, t;
- saidx_t limit;
- saidx_t v, x = 0;
- int ssize;
-
- for(ssize = 0, limit = _lg(last - first);;) {
-
- if((last - first) <= LS_INSERTIONSORT_THRESHOLD) {
- if(1 < (last - first)) {
- _insertionsort(ISAd, first, last);
- _ls_updategroup(ISA, SA, first, last);
- } else if((last - first) == 1) { *first = -1; }
- STACK_POP3(first, last, limit);
- continue;
- }
-
- if(limit-- == 0) {
- _heapsort(ISAd, first, last - first);
- for(a = last - 1; first < a; a = b) {
- for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
- }
- _ls_updategroup(ISA, SA, first, last);
- STACK_POP3(first, last, limit);
- continue;
- }
-
- /* choose pivot */
- a = _pivot(ISAd, first, last);
- SWAP(*first, *a);
- v = ISAd[*first];
-
- /* Two-stage double-index controlled ternary partition */
- for(b = first; (++b < last) && ((x = ISAd[*b]) == v);) { }
- if(((a = b) < last) && (x < v)) {
- for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- }
- for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { }
- if((b < (d = c)) && (x > v)) {
- for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
- for(; b < c;) {
- SWAP(*b, *c);
- for(; (++b < c) && ((x = ISAd[*b]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
-
- if(a <= d) {
- c = b - 1;
-
- if((s = a - first) > (t = b - a)) { s = t; }
- for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
- if((s = d - c) > (t = last - d - 1)) { s = t; }
- for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
-
- a = first + (b - a), b = last - (d - c);
-
- /* update ranks */
- for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
- if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } }
- if((b - a) == 1) { *a = - 1; }
-
- if((a - first) <= (last - b)) {
- if(first < a) {
- STACK_PUSH3(b, last, limit);
- last = a;
- } else {
- first = b;
- }
- } else {
- if(b < last) {
- STACK_PUSH3(first, a, limit);
- first = b;
- } else {
- last = a;
- }
- }
- } else {
- STACK_POP3(first, last, limit);
- }
- }
+trbudget_init(trbudget_t *budget, saidx_t chance, saidx_t incval) {
+ budget->chance = chance;
+ budget->remain = budget->incval = incval;
}
-
-/*---------------------------------------------------------------------------*/
-
-/* Larsson-Sadakane sort */
-static
-void
-_lssort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth) {
- saidx_t *ISAd;
- saidx_t *first, *last;
- saidx_t t, skip;
-
- for(ISAd = ISA + depth; -n < *SA; ISAd += (ISAd - ISA)) {
- first = SA;
- skip = 0;
- do {
- if((t = *first) < 0) { first -= t; skip += t; }
- else {
- if(skip != 0) { *(first + skip) = skip; skip = 0; }
- last = SA + ISA[t] + 1;
- _ls_introsort(ISA, ISAd, SA, first, last);
- first = last;
- }
- } while(first < (SA + n));
- if(skip != 0) { *(first + skip) = skip; }
- }
+static INLINE
+saint_t
+trbudget_check(trbudget_t *budget, saidx_t size) {
+ if(size <= budget->remain) { budget->remain -= size; return 1; }
+ if(budget->chance == 0) { budget->count += size; return 0; }
+ budget->remain += budget->incval - size;
+ budget->chance -= 1;
+ return 1;
}
/*---------------------------------------------------------------------------*/
-static
+static INLINE
void
-_tr_partition(const saidx_t *ISAd, const saidx_t *SA,
- saidx_t *first, saidx_t *last,
- saidx_t **pa, saidx_t **pb, saidx_t v) {
+tr_partition(const saidx_t *ISAd,
+ saidx_t *first, saidx_t *middle, saidx_t *last,
+ saidx_t **pa, saidx_t **pb, saidx_t v) {
saidx_t *a, *b, *c, *d, *e, *f;
saidx_t t, s;
saidx_t x = 0;
- for(b = first - 1; (++b < last) && ((x = ISAd[*b]) == v);) { }
+ for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { }
if(((a = b) < last) && (x < v)) {
for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
if(x == v) { SWAP(*b, *a); ++a; }
@@ -377,9 +261,9 @@ _tr_partition(const saidx_t *ISAd, const saidx_t *SA,
static
void
-_tr_copy(saidx_t *ISA, const saidx_t *SA,
- saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last,
- saidx_t depth) {
+tr_copy(saidx_t *ISA, const saidx_t *SA,
+ saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last,
+ saidx_t depth) {
/* sort suffixes of middle partition
by using sorted order of suffixes of left and right partition. */
saidx_t *c, *d, *e;
@@ -400,34 +284,64 @@ _tr_copy(saidx_t *ISA, const saidx_t *SA,
}
}
-/* Multikey introsort. */
static
void
-_tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
- saidx_t *SA, saidx_t *first, saidx_t *last,
- saidx_t *budget, saidx_t *chance, saidx_t size) {
-#define UPDATE_BUDGET(n)\
- {\
- (*budget) -= (n);\
- if(*budget <= 0) {\
- if(--(*chance) == 0) { break; }\
- (*budget) += size;\
- }\
+tr_partialcopy(saidx_t *ISA, const saidx_t *SA,
+ saidx_t *first, saidx_t *a, saidx_t *b, saidx_t *last,
+ saidx_t depth) {
+ saidx_t *c, *d, *e;
+ saidx_t s, v;
+ saidx_t rank, lastrank, newrank = -1;
+
+ v = b - SA - 1;
+ lastrank = -1;
+ for(c = first, d = a - 1; c <= d; ++c) {
+ if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
+ *++d = s;
+ rank = ISA[s + depth];
+ if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
+ ISA[s] = newrank;
+ }
}
- struct { const saidx_t *a; saidx_t *b, *c, d; }stack[STACK_SIZE];
- saidx_t *a, *b, *c, *d, *e, *f;
- saidx_t s, t;
+
+ lastrank = -1;
+ for(e = d; first <= e; --e) {
+ rank = ISA[*e];
+ if(lastrank != rank) { lastrank = rank; newrank = e - SA; }
+ if(newrank != rank) { ISA[*e] = newrank; }
+ }
+
+ lastrank = -1;
+ for(c = last - 1, e = d + 1, d = b; e < d; --c) {
+ if((0 <= (s = *c - depth)) && (ISA[s] == v)) {
+ *--d = s;
+ rank = ISA[s + depth];
+ if(lastrank != rank) { lastrank = rank; newrank = d - SA; }
+ ISA[s] = newrank;
+ }
+ }
+}
+
+static
+void
+tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
+ saidx_t *SA, saidx_t *first, saidx_t *last,
+ trbudget_t *budget) {
+#define STACK_SIZE TR_STACKSIZE
+ struct { const saidx_t *a; saidx_t *b, *c; saint_t d, e; }stack[STACK_SIZE];
+ saidx_t *a, *b, *c;
+ saidx_t t;
saidx_t v, x = 0;
- saidx_t limit, next;
- int ssize;
+ saidx_t incr = ISAd - ISA;
+ saint_t limit, next;
+ saint_t ssize, trlink = -1;
- for(ssize = 0, limit = _lg(last - first);;) {
+ for(ssize = 0, limit = tr_ilg(last - first);;) {
if(limit < 0) {
if(limit == -1) {
/* tandem repeat partition */
- UPDATE_BUDGET(last - first);
- _tr_partition(ISAd - 1, SA, first, last, &a, &b, last - SA - 1);
+ tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1);
/* update ranks */
if(a < last) {
@@ -438,32 +352,40 @@ _tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
}
/* push */
- STACK_PUSH(NULL, a, b, 0);
- STACK_PUSH(ISAd - 1, first, last, -2);
+ if(1 < (b - a)) {
+ STACK_PUSH5(NULL, a, b, 0, 0);
+ STACK_PUSH5(ISAd - incr, first, last, -2, trlink);
+ trlink = ssize - 2;
+ }
if((a - first) <= (last - b)) {
if(1 < (a - first)) {
- STACK_PUSH(ISAd, b, last, _lg(last - b));
- last = a, limit = _lg(a - first);
+ STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink);
+ last = a, limit = tr_ilg(a - first);
} else if(1 < (last - b)) {
- first = b, limit = _lg(last - b);
+ first = b, limit = tr_ilg(last - b);
} else {
- STACK_POP(ISAd, first, last, limit);
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
} else {
if(1 < (last - b)) {
- STACK_PUSH(ISAd, first, a, _lg(a - first));
- first = b, limit = _lg(last - b);
+ STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink);
+ first = b, limit = tr_ilg(last - b);
} else if(1 < (a - first)) {
- last = a, limit = _lg(a - first);
+ last = a, limit = tr_ilg(a - first);
} else {
- STACK_POP(ISAd, first, last, limit);
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
}
} else if(limit == -2) {
/* tandem repeat copy */
a = stack[--ssize].b, b = stack[ssize].c;
- _tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
- STACK_POP(ISAd, first, last, limit);
+ if(stack[ssize].d == 0) {
+ tr_copy(ISA, SA, first, a, b, last, ISAd - ISA);
+ } else {
+ if(0 <= trlink) { stack[trlink].d = -1; }
+ tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA);
+ }
+ STACK_POP5(ISAd, first, last, limit, trlink);
} else {
/* sorted partition */
if(0 <= *first) {
@@ -473,40 +395,45 @@ _tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
}
if(first < last) {
a = first; do { *a = ~*a; } while(*++a < 0);
- next = (ISA[*a] != ISAd[*a]) ? _lg(a - first + 1) : -1;
+ next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1;
if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } }
/* push */
- if((a - first) <= (last - a)) {
- STACK_PUSH(ISAd, a, last, -3);
- ISAd += 1, last = a, limit = next;
+ if(trbudget_check(budget, a - first)) {
+ if((a - first) <= (last - a)) {
+ STACK_PUSH5(ISAd, a, last, -3, trlink);
+ ISAd += incr, last = a, limit = next;
+ } else {
+ if(1 < (last - a)) {
+ STACK_PUSH5(ISAd + incr, first, a, next, trlink);
+ first = a, limit = -3;
+ } else {
+ ISAd += incr, last = a, limit = next;
+ }
+ }
} else {
+ if(0 <= trlink) { stack[trlink].d = -1; }
if(1 < (last - a)) {
- STACK_PUSH(ISAd + 1, first, a, next);
first = a, limit = -3;
} else {
- ISAd += 1, last = a, limit = next;
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
}
} else {
- STACK_POP(ISAd, first, last, limit);
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
}
continue;
}
if((last - first) <= TR_INSERTIONSORT_THRESHOLD) {
- UPDATE_BUDGET(last - first);
-// UPDATE_BUDGET((last - first) * (last - first) / 2);
- _insertionsort(ISAd, first, last);
+ tr_insertionsort(ISAd, first, last);
limit = -3;
continue;
}
if(limit-- == 0) {
- UPDATE_BUDGET(last - first);
-// UPDATE_BUDGET((last - first) * _lg(last - first));
- _heapsort(ISAd, first, last - first);
+ tr_heapsort(ISAd, first, last - first);
for(a = last - 1; first < a; a = b) {
for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; }
}
@@ -514,136 +441,111 @@ _tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
continue;
}
- UPDATE_BUDGET(last - first);
-
/* choose pivot */
- a = _pivot(ISAd, first, last);
+ a = tr_pivot(ISAd, first, last);
SWAP(*first, *a);
v = ISAd[*first];
/* partition */
- for(b = first; (++b < last) && ((x = ISAd[*b]) == v);) { }
- if(((a = b) < last) && (x < v)) {
- for(; (++b < last) && ((x = ISAd[*b]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- }
- for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { }
- if((b < (d = c)) && (x > v)) {
- for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
- for(; b < c;) {
- SWAP(*b, *c);
- for(; (++b < c) && ((x = ISAd[*b]) <= v);) {
- if(x == v) { SWAP(*b, *a); ++a; }
- }
- for(; (b < --c) && ((x = ISAd[*c]) >= v);) {
- if(x == v) { SWAP(*c, *d); --d; }
- }
- }
-
- if(a <= d) {
- c = b - 1;
-
- if((s = a - first) > (t = b - a)) { s = t; }
- for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
- if((s = d - c) > (t = last - d - 1)) { s = t; }
- for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); }
-
- a = first + (b - a), b = last - (d - c);
- next = (ISA[*a] != v) ? _lg(b - a) : -1;
+ tr_partition(ISAd, first, first + 1, last, &a, &b, v);
+ if((last - first) != (b - a)) {
+ next = (ISA[*a] != v) ? tr_ilg(b - a) : -1;
/* update ranks */
for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; }
if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } }
/* push */
- if((a - first) <= (last - b)) {
- if((last - b) <= (b - a)) {
- if(1 < (a - first)) {
- STACK_PUSH(ISAd + 1, a, b, next);
- STACK_PUSH(ISAd, b, last, limit);
- last = a;
- } else if(1 < (last - b)) {
- STACK_PUSH(ISAd + 1, a, b, next);
- first = b;
- } else if(1 < (b - a)) {
- ISAd += 1, first = a, last = b, limit = next;
- } else {
- STACK_POP(ISAd, first, last, limit);
- }
- } else if((a - first) <= (b - a)) {
- if(1 < (a - first)) {
- STACK_PUSH(ISAd, b, last, limit);
- STACK_PUSH(ISAd + 1, a, b, next);
- last = a;
- } else if(1 < (b - a)) {
- STACK_PUSH(ISAd, b, last, limit);
- ISAd += 1, first = a, last = b, limit = next;
+ if((1 < (b - a)) && (trbudget_check(budget, b - a))) {
+ if((a - first) <= (last - b)) {
+ if((last - b) <= (b - a)) {
+ if(1 < (a - first)) {
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
+ last = a;
+ } else if(1 < (last - b)) {
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ first = b;
+ } else {
+ ISAd += incr, first = a, last = b, limit = next;
+ }
+ } else if((a - first) <= (b - a)) {
+ if(1 < (a - first)) {
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ last = a;
+ } else {
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
+ ISAd += incr, first = a, last = b, limit = next;
+ }
} else {
- first = b;
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
+ ISAd += incr, first = a, last = b, limit = next;
}
} else {
- if(1 < (b - a)) {
- STACK_PUSH(ISAd, b, last, limit);
- STACK_PUSH(ISAd, first, a, limit);
- ISAd += 1, first = a, last = b, limit = next;
+ if((a - first) <= (b - a)) {
+ if(1 < (last - b)) {
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
+ first = b;
+ } else if(1 < (a - first)) {
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ last = a;
+ } else {
+ ISAd += incr, first = a, last = b, limit = next;
+ }
+ } else if((last - b) <= (b - a)) {
+ if(1 < (last - b)) {
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
+ STACK_PUSH5(ISAd + incr, a, b, next, trlink);
+ first = b;
+ } else {
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
+ ISAd += incr, first = a, last = b, limit = next;
+ }
} else {
- STACK_PUSH(ISAd, b, last, limit);
- last = a;
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
+ ISAd += incr, first = a, last = b, limit = next;
}
}
} else {
- if((a - first) <= (b - a)) {
- if(1 < (last - b)) {
- STACK_PUSH(ISAd + 1, a, b, next);
- STACK_PUSH(ISAd, first, a, limit);
- first = b;
- } else if(1 < (a - first)) {
- STACK_PUSH(ISAd + 1, a, b, next);
+ if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; }
+ if((a - first) <= (last - b)) {
+ if(1 < (a - first)) {
+ STACK_PUSH5(ISAd, b, last, limit, trlink);
last = a;
- } else if(1 < (b - a)) {
- ISAd += 1, first = a, last = b, limit = next;
+ } else if(1 < (last - b)) {
+ first = b;
} else {
- STACK_POP(ISAd, first, last, limit);
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
- } else if((last - b) <= (b - a)) {
+ } else {
if(1 < (last - b)) {
- STACK_PUSH(ISAd, first, a, limit);
- STACK_PUSH(ISAd + 1, a, b, next);
+ STACK_PUSH5(ISAd, first, a, limit, trlink);
first = b;
- } else if(1 < (b - a)) {
- STACK_PUSH(ISAd, first, a, limit);
- ISAd += 1, first = a, last = b, limit = next;
- } else {
+ } else if(1 < (a - first)) {
last = a;
- }
- } else {
- if(1 < (b - a)) {
- STACK_PUSH(ISAd, first, a, limit);
- STACK_PUSH(ISAd, b, last, limit);
- ISAd += 1, first = a, last = b, limit = next;
} else {
- STACK_PUSH(ISAd, first, a, limit);
- first = b;
+ STACK_POP5(ISAd, first, last, limit, trlink);
}
}
}
} else {
- limit += 1, ISAd += 1;
- }
- }
-
- for(s = 0; s < ssize; ++s) {
- if(stack[s].d == -3) {
- _ls_updategroup(ISA, SA, stack[s].b, stack[s].c);
+ if(trbudget_check(budget, last - first)) {
+ limit = tr_ilg(last - first), ISAd += incr;
+ } else {
+ if(0 <= trlink) { stack[trlink].d = -1; }
+ STACK_POP5(ISAd, first, last, limit, trlink);
+ }
}
}
+#undef STACK_SIZE
}
+
/*---------------------------------------------------------------------------*/
/*- Function -*/
@@ -651,34 +553,34 @@ _tr_introsort(saidx_t *ISA, const saidx_t *ISAd,
/* Tandem repeat sort */
void
trsort(saidx_t *ISA, saidx_t *SA, saidx_t n, saidx_t depth) {
+ saidx_t *ISAd;
saidx_t *first, *last;
- saidx_t t, skip;
- saidx_t budget;
- saidx_t chance;
+ trbudget_t budget;
+ saidx_t t, skip, unsorted;
- if(-n < *SA) {
+ trbudget_init(&budget, tr_ilg(n) * 2 / 3, n);
+/* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */
+ for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) {
first = SA;
skip = 0;
- budget = n;
-/* chance = _lg(n); */
-/* chance = _lg(n) / 2 + 1; */
- chance = _lg(n) * 2 / 3 + 1;
+ unsorted = 0;
do {
if((t = *first) < 0) { first -= t; skip += t; }
else {
- skip = 0;
+ if(skip != 0) { *(first + skip) = skip; skip = 0; }
last = SA + ISA[t] + 1;
if(1 < (last - first)) {
- _tr_introsort(ISA, ISA + depth, SA, first, last, &budget, &chance, n);
- if(chance == 0) {
- /* Switch to Larsson-Sadakane sorting algorithm. */
- if(SA < first) { *SA = -(first - SA); }
- _lssort(ISA, SA, n, depth);
- break;
- }
+ budget.count = 0;
+ tr_introsort(ISA, ISAd, SA, first, last, &budget);
+ if(budget.count != 0) { unsorted += budget.count; }
+ else { skip = first - last; }
+ } else if((last - first) == 1) {
+ skip = -1;
}
first = last;
}
} while(first < (SA + n));
+ if(skip != 0) { *(first + skip) = skip; }
+ if(unsorted == 0) { break; }
}
}
diff --git a/lib/utils.c b/lib/utils.c
index 086b1ca..75bd274 100644
--- a/lib/utils.c
+++ b/lib/utils.c
@@ -32,17 +32,17 @@
/* Binary search for inverse bwt. */
static
saidx_t
-_binarysearch(const saidx_t *A, saidx_t len, saidx_t val) {
- saidx_t half, m;
- for(m = 0, half = len >> 1;
- 0 < len;
- len = half, half >>= 1) {
- if(A[m + half] < val) {
- m += half + 1;
- half -= ((len & 1) == 0);
+binarysearch_lower(const saidx_t *A, saidx_t size, saidx_t value) {
+ saidx_t half, i;
+ for(i = 0, half = size >> 1;
+ 0 < size;
+ size = half, half >>= 1) {
+ if(A[i + half] < value) {
+ i += half + 1;
+ half -= (size & 1) ^ 1;
}
}
- return m;
+ return i;
}
@@ -52,12 +52,13 @@ _binarysearch(const saidx_t *A, saidx_t len, saidx_t val) {
saint_t
bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *SA,
saidx_t n, saidx_t *idx) {
- saidx_t *A, i, p;
+ saidx_t *A, i, j, p, t;
saint_t c;
/* Check arguments. */
if((T == NULL) || (U == NULL) || (n < 0) || (idx == NULL)) { return -1; }
if(n <= 1) {
+ if(n == 1) { U[0] = T[0]; }
*idx = n;
return 0;
}
@@ -70,22 +71,32 @@ bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *SA,
/* BW transform. */
if(T == U) {
- for(i = 0; 0 <= (p = A[i] - 1); ++i) {
- c = T[i];
- U[i] = (i <= p) ? T[p] : A[p];
- A[i] = c;
+ t = n;
+ for(i = 0, j = 0; i < n; ++i) {
+ p = t - 1;
+ t = A[i];
+ if(0 <= p) {
+ c = T[j];
+ U[j] = (j <= p) ? T[p] : (sauchar_t)A[p];
+ A[j] = c;
+ j++;
+ } else {
+ *idx = i;
+ }
}
- *idx = i;
- for( ; i < n; ++i) {
- p = A[i + 1] - 1;
- c = T[i];
- U[i] = (i <= p) ? T[p] : A[p];
- A[i] = c;
+ p = t - 1;
+ if(0 <= p) {
+ c = T[j];
+ U[j] = (j <= p) ? T[p] : (sauchar_t)A[p];
+ A[j] = c;
+ } else {
+ *idx = i;
}
} else {
- for(i = 0; A[i] != 0; ++i) { U[i] = T[A[i] - 1]; }
- *idx = i;
- for( ; i < n; ++i) { U[i] = T[A[i + 1] - 1]; }
+ U[0] = T[n - 1];
+ for(i = 0; A[i] != 0; ++i) { U[i + 1] = T[A[i] - 1]; }
+ *idx = i + 1;
+ for(++i; i < n; ++i) { U[i] = T[A[i] - 1]; }
}
if(SA == NULL) {
@@ -101,9 +112,10 @@ saint_t
inverse_bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *A,
saidx_t n, saidx_t idx) {
saidx_t C[ALPHABET_SIZE];
- saidx_t D[ALPHABET_SIZE];
+ sauchar_t D[ALPHABET_SIZE];
saidx_t *B;
- saidx_t i, k, t, sum;
+ saidx_t i, p;
+ saint_t c, d;
/* Check arguments. */
if((T == NULL) || (U == NULL) || (n < 0) || (idx < 0) ||
@@ -113,30 +125,27 @@ inverse_bw_transform(const sauchar_t *T, sauchar_t *U, saidx_t *A,
if(n <= 1) { return 0; }
if((B = A) == NULL) {
- /* Allocate (n+1)*sizeof(saidx_t) bytes of memory. */
- if((B = malloc((n + 1) * sizeof(saidx_t))) == NULL) { return -2; }
+ /* Allocate n*sizeof(saidx_t) bytes of memory. */
+ if((B = malloc(n * sizeof(saidx_t))) == NULL) { return -2; }
}
/* Inverse BW transform. */
- for(i = 0; i < ALPHABET_SIZE; ++i) { C[i] = 0; }
+ for(c = 0; c < ALPHABET_SIZE; ++c) { C[c] = 0; }
for(i = 0; i < n; ++i) { ++C[T[i]]; }
- for(i = 0, sum = 0; i < ALPHABET_SIZE; ++i) {
- t = C[i];
- C[i] = sum;
- sum += t;
- }
- B[0] = idx;
- for(i = 0; i < idx; ++i) { B[++C[T[i]]] = i; }
- for( ; i < n; ++i) { B[++C[T[i]]] = i + 1; }
- for(i = 0, k = 0, t = 0; i < ALPHABET_SIZE; ++i) {
- if(t != C[i]) {
- D[k] = i;
- t = C[k] = C[i];
- ++k;
+ for(c = 0, d = 0, i = 0; c < ALPHABET_SIZE; ++c) {
+ p = C[c];
+ if(0 < p) {
+ C[c] = i;
+ D[d++] = (sauchar_t)c;
+ i += p;
}
}
- for(i = 0, t = 0; i < n; ++i) {
- U[i] = D[_binarysearch(C, k, t = B[t])];
+ for(i = 0; i < idx; ++i) { B[C[T[i]]++] = i; }
+ for( ; i < n; ++i) { B[C[T[i]]++] = i + 1; }
+ for(c = 0; c < d; ++c) { C[c] = C[D[c]]; }
+ for(i = 0, p = idx; i < n; ++i) {
+ U[i] = D[binarysearch_lower(C, d, p)];
+ p = B[p - 1];
}
if(A == NULL) {
@@ -152,89 +161,83 @@ saint_t
sufcheck(const sauchar_t *T, const saidx_t *SA,
saidx_t n, saint_t verbose) {
saidx_t C[ALPHABET_SIZE];
- saidx_t i = 0, p, t = 0;
+ saidx_t i, p, q, t;
saint_t c;
- saint_t err = 0;
- if(1 <= verbose) { fprintf(stderr, "sufchecker: "); }
+ if(verbose) { fprintf(stderr, "sufcheck: "); }
/* Check arguments. */
- if((T == NULL) || (SA == NULL) || (n < 0)) { err = -1; }
-
- /* ranges. */
- if(err == 0) {
- for(i = 0; i <= n; ++i) {
- if((SA[i] < 0) || (n < SA[i])) {
- err = -2;
- break;
+ if((T == NULL) || (SA == NULL) || (n < 0)) {
+ if(verbose) { fprintf(stderr, "Invalid arguments.\n"); }
+ return -1;
+ }
+ if(n == 0) {
+ if(verbose) { fprintf(stderr, "Done.\n"); }
+ return 0;
+ }
+
+ /* check range: [0..n-1] */
+ for(i = 0; i < n; ++i) {
+ if((SA[i] < 0) || (n <= SA[i])) {
+ if(verbose) {
+ fprintf(stderr, "Out of the range [0,%" PRIdSAIDX_T "].\n"
+ " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "\n",
+ n - 1, i, SA[i]);
}
+ return -2;
}
}
- /* first characters. */
- if(err == 0) {
- for(i = 1; i < n; ++i) {
- if(T[SA[i]] > T[SA[i + 1]]) {
- err = -3;
- break;
+ /* check first characters. */
+ for(i = 1; i < n; ++i) {
+ if(T[SA[i - 1]] > T[SA[i]]) {
+ if(verbose) {
+ fprintf(stderr, "Suffixes in wrong order.\n"
+ " T[SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "]=%d"
+ " > T[SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "]=%d\n",
+ i - 1, SA[i - 1], T[SA[i - 1]], i, SA[i], T[SA[i]]);
}
+ return -3;
}
}
- /* suffixes. */
- if(err == 0) {
- for(i = 0; i < ALPHABET_SIZE; ++i) { C[i] = 0; }
- for(i = 0; i < n; ++i) { ++C[T[i]]; }
- for(i = 0, p = 1; i < ALPHABET_SIZE; ++i) {
- t = C[i];
- C[i] = p;
- p += t;
- }
+ /* check suffixes. */
+ for(i = 0; i < ALPHABET_SIZE; ++i) { C[i] = 0; }
+ for(i = 0; i < n; ++i) { ++C[T[i]]; }
+ for(i = 0, p = 0; i < ALPHABET_SIZE; ++i) {
+ t = C[i];
+ C[i] = p;
+ p += t;
+ }
- for(i = 0; i <= n; ++i) {
- p = SA[i];
- if(0 < p) {
- c = T[--p];
- t = C[c];
- } else {
- p = n;
- c = -1;
- t = 0;
- }
- if(p != SA[t]) {
- err = -4;
- break;
- }
- if(0 <= c) {
- ++C[c];
- if((n < C[c]) || (T[SA[C[c]]] != c)) { C[c] = -1; }
+ q = C[T[n - 1]];
+ C[T[n - 1]] += 1;
+ for(i = 0; i < n; ++i) {
+ p = SA[i];
+ if(0 < p) {
+ c = T[--p];
+ t = C[c];
+ } else {
+ c = T[p = n - 1];
+ t = q;
+ }
+ if((t < 0) || (p != SA[t])) {
+ if(verbose) {
+ fprintf(stderr, "Suffix in wrong position.\n"
+ " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T " or\n"
+ " SA[%" PRIdSAIDX_T "]=%" PRIdSAIDX_T "\n",
+ t, (0 <= t) ? SA[t] : -1, i, SA[i]);
}
+ return -4;
}
- }
-
- if(1 <= verbose) {
- if(err == 0) {
- fprintf(stderr, "Done.\n");
- } else if(verbose == 1) {
- fprintf(stderr, "Error.\n");
- } else if(err == -1) {
- fprintf(stderr, "Invalid arguments.\n");
- } else if(err == -2) {
- fprintf(stderr, "Out of the range [0,%d].\n SA[%d]=%d\n",
- (int)n, (int)i, (int)SA[i]);
- } else if(err == -3) {
- fprintf(stderr, "Suffixes in wrong order.\n"
- " T[SA[%d]=%d]=%d > T[SA[%d]=%d]=%d\n",
- (int)i, (int)SA[i], (int)T[SA[i]],
- (int)i + 1, (int)SA[i + 1], (int)T[SA[i + 1]]);
- } else if(err == -4) {
- fprintf(stderr, "Suffix in wrong position.\n");
- if(0 <= t) { fprintf(stderr, " SA[%d]=%d or\n", (int)t, (int)SA[t]); }
- fprintf(stderr, " SA[%d]=%d\n", (int)i, (int)SA[i]);
+ if(t != q) {
+ ++C[c];
+ if((n <= C[c]) || (T[SA[C[c]]] != c)) { C[c] = -1; }
}
}
- return err;
+ if(1 <= verbose) { fprintf(stderr, "Done.\n"); }
+ return 0;
}
diff --git a/pkgconfig/CMakeLists.txt b/pkgconfig/CMakeLists.txt
new file mode 100644
index 0000000..c9e57ab
--- /dev/null
+++ b/pkgconfig/CMakeLists.txt
@@ -0,0 +1,3 @@
+## generate libdivsufsort.pc ##
+configure_file("${CMAKE_CURRENT_SOURCE_DIR}/libdivsufsort.pc.cmake" "${CMAKE_CURRENT_BINARY_DIR}/libdivsufsort.pc" @ONLY)
+install(FILES "${CMAKE_CURRENT_BINARY_DIR}/libdivsufsort.pc" DESTINATION lib/pkgconfig)
diff --git a/pkgconfig/libdivsufsort.pc.cmake b/pkgconfig/libdivsufsort.pc.cmake
new file mode 100644
index 0000000..878e47e
--- /dev/null
+++ b/pkgconfig/libdivsufsort.pc.cmake
@@ -0,0 +1,11 @@
+prefix=@CMAKE_INSTALL_PREFIX@
+exec_prefix=${prefix}
+libdir=${exec_prefix}/lib
+includedir=${prefix}/include
+
+Name: @PROJECT_NAME@
+Description: @PROJECT_DESCRIPTION@
+Version: @PROJECT_VERSION_FULL@
+URL: @PROJECT_URL@
+Libs: -L${libdir} -ldivsufsort
+Cflags: -I${includedir}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-26 15:35:17 +0000