Leahnaya/Flips
1
0
Fork 0
mirror of https://github.com/Alcaro/Flips.git synced 2026-09-07 10:06:09 -05:00
Code Issues Actions Packages Projects Releases Wiki Activity

Rewrite this one to Arlib API, and wipe out moremem

Browse Source
This commit is contained in:
Alcaro
2016-12-20 22:19:58 +01:00
parent 8d1780bb29
commit 0e31df85fd
6 changed files with 101 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files

206
patch/bps-create.cppx → patch/bps-create.cpp
View File

@@ -1,10 +1,4 @@
#include "libbps.h"
#include "arlib/crc32.h"
#include "arlib/file.h"
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "patch.h"
//These two give minor performance penalties and will print some random stuff to stdout.
//The former will verify the correctness of the output patch, the latter will print some performance data.
@@ -19,10 +13,6 @@
//If it's something else, get a non-broken array calculator.
#define EOF_IS_LAST false
#if defined(TEST_CORRECT) || defined(TEST_PERF)
#include <stdio.h>
#endif
//Algorithm description:
//
//This is heavily built upon suffix sorting; the implementation I use, libdivsufsort, claims
@@ -31,10 +21,14 @@
//
//The program starts by taking an equal amount of the source file and target file, concatenates that
// with target first, and suffix sorts it.
//It also calculates a reverse index, such that reverse[sorted[i]]==i.
//
//To find a match, it goes to reverse[outpos], and scans sorted[] up and down for the closest entry
// that either starts before the current output position, or is somewhere in the source file.
//To find a match, it finds the sortpos where sorted[sortpos]==outpos. This is a binary search; it's
// called O(n) times, with O(log n) comparisons per iteration. Each comparison is potentially O(n),
// but for each matched byte, another iteration is removed from the outer loop, so the comparisons
// can be considered O(1) each; the sum is O(n log n).
//
//After it's found sortpos, it scans sorted[] up and down for the closest entry that either starts
// before the current output position, or is somewhere in the source file.
//As the source file comes last, the end-of-file marker (whose value is outside the range of a byte)
// is guaranteed to not be in the way for a better match.
//This is called O(n) times, and averages O(1) as at least 50% of sorted[] is in range. However, it
@@ -42,8 +36,8 @@
//
//It then checks which of the two candidates are superior, by checking how far they match each
// other, and then checking if the upper one has another correct byte.
//This is potentially O(n), but for each matched byte, another iteration is removed from the outer
// loop, so the sum of all calls is O(n).
//This is potentially O(n), but like the binary search, long matches reduce the outer loop. The sum
// is O(n).
//
//When the program approaches the end of the sorted area, it re-sorts twice as much as last time.
// This gives O(log n) calls to the suffix sorter.
@@ -53,15 +47,15 @@
//
//Many details were omitted from the above, but that's the basic setup.
//
//Thus, the program is O(max(n log n, n, n) = n log n) average and O(max(n log n, n^2, n) = n^2)
// worst case.
//Thus, the program is O(max(n log n, n log n, n, n) = n log n) average and O(max(n log n, n log n,
// n^2, n) = n^2) worst case.
//
//I conclude that the task of finding, understanding and implementing a sub-O(n^2) algorithm for
// delta patching is resolved.
//Known cases where this function does not emit the optimal encoding:
//If a match in the target file would extend further than target_search_size, it is often skipped.
//If a match in the target file would extend further than target_search_size, it's cut off.
// Penalty: O(log n), with extremely low constants (it'd require a >256B match to be exactly there).
// Even for big files, the penalty is very likely to remain zero; even hitting double-digit bytes
// would require a file designed exactly for that.
@@ -70,20 +64,23 @@
//However, due to better heuristics and others' performance optimizations, this one still beats its
// competitors.
//TODO: test multiple same-length matches
// but only for lengths <= 64,
//Possible optimizations:
//divsufsort() takes approximately 2/3 of the total time. create_reverse_index() takes roughly a third of the remainder.
//divsufsort() takes approximately 1/2 of the total time.
//Each iteration takes four times as long as the previous one.
//If each iteration takes 4 times as long as the previous one, then the last one takes 3/4 of the total time.
//Since divsufsort+create_reverse_index doesn't depend on anything else, the last iteration can be split off to its own thread.
//Since divsufsort doesn't depend on anything else, the last iteration can be split off to its own thread.
//This would split it to
//Search, non-final: 2/9 * 1/4 = 2/36
//Search, final: 2/9 * 3/4 = 6/36
//Sort+rev, non-final: 7/9 * 1/4 = 7/36
//Sort+rev, final: 7/9 * 3/4 = 21/36
//Search, non-final: 1/2 * 1/4 = 1/8
//Search, final: 1/2 * 3/4 = 3/8
//Sort+rev, non-final: 1/2 * 1/4 = 1/8
//Sort+rev, final: 1/2 * 3/4 = 3/8
//All non-final must be done sequentially. Both Sort Final and non-final must be done before Search Final can start.
//This means the final time, if Sort Final is split off, is
//max(7/36+2/36, 21/36) + 6/36 = 27/36 = 3/4
//max(1/8+1/8, 3/8) + 3/8 = 6/8 = 3/4
//of the original time.
//Due to
//- the considerable complexity costs (OpenMP doesn't seem able to represent the "insert a wait in
@@ -94,13 +91,13 @@
// and that the small ones are not, as that'd starve the big one
//I deem a possible 25% boost not worthwhile.
//Both sorting algorithms claim O(1) memory use (in addition to the bytes and the output). In
// addition to that, this algorithm uses (source.len*target.len)*(sizeof(uint8_t)+2*sizeof(off_t))
// bytes of memory, plus the input and output files, plus the patch.
//For most hardware, this is 9*(source.len+target.len), or 5*(source+target) for the slim one.
//Both sorting algorithms claim O(1) memory use, in addition to the in/outputs. For most hardware,
// this is 5*(source.len+target.len).
//If the output is stored to disk, that's all this algorithm needs as well.
namespace patch { namespace bps {
//TODO: HEAVY cleanups needed here
#include "sais.cpp"
template<typename sais_index_type>
static void sufsort(sais_index_type* SA, const uint8_t* T, sais_index_type n) {
@@ -116,7 +113,7 @@ static void sufsort(sais_index_type* SA, const uint8_t* T, sais_index_type n) {
//I'd prefer to let them allocate from an array I give it, but divsuf doesn't allow that, and there
// are only half a dozen allocations per call anyways.
//This ends up in libdivsufsort if available, otherwise lite.
//This ends up in libdivsufsort if available, otherwise sais.cpp.
#include "divsufsort.h"
static void sufsort(int32_t* SA, uint8_t* T, int32_t n)
{
@@ -140,6 +137,18 @@ template<typename T> static T max(T a, T b) { return a<b ? b : a; }
namespace {
//class filecache {
// file& f;
// uint32_t crc32;
// int bytes_used;
// uint8_t bytes[65536];
//
// void append(const uint8_t * data, size_t len)
// {
//
// }
//};
struct bps_creator {
uint8_t* out;
size_t outlen;
@@ -206,7 +215,7 @@ struct bps_creator {
size_t numtargetread;
bps_creator(file* source, file* target, struct mem metadata)
bps_creator(const file& source, const file& target, struct mem metadata)
{
outlen = 0;
outbuflen = 128;
@@ -220,12 +229,10 @@ struct bps_creator {
numtargetread = 0;
append((const uint8_t*)"BPS1", 4);
appendnum(source->len);
appendnum(target->len);
appendnum(source.size());
appendnum(target.size());
appendnum(metadata.len);
append(metadata.ptr, metadata.len);
setProgress(NULL, NULL);
}
@@ -342,23 +349,7 @@ struct bps_creator {
}
bool (*prog_func)(void* userdata, size_t done, size_t total);
void* prog_dat;
static bool prog_func_null(void* userdata, size_t done, size_t total) { return true; }
void setProgress(bool (*progress)(void* userdata, size_t done, size_t total), void* userdata)
{
if (!progress) progress = prog_func_null;
prog_func=progress;
prog_dat=userdata;
}
bool progress(size_t done, size_t total)
{
return prog_func(prog_dat, done, total);
}
function<bool(size_t done, size_t total)> progress;
void finish(const uint8_t* source, const uint8_t* target)
@@ -369,9 +360,9 @@ struct bps_creator {
puts("ERROR: patch creates wrong ROM size"),abort();
#endif
appendnum32(crc32(source, sourcelen));
appendnum32(crc32(target, targetlen));
appendnum32(crc32(out, outlen));
appendnum32(crc32(arrayview<byte>(source, sourcelen)));
appendnum32(crc32(arrayview<byte>(target, targetlen)));
appendnum32(crc32(arrayview<byte>(out, outlen)));
}
struct mem getpatch()
@@ -545,10 +536,8 @@ static void create_buckets(const uint8_t* data, off_t* index, off_t len, off_t*
}
template<typename off_t>
static off_t find_index(off_t pos, const uint8_t* data, off_t datalen, const off_t* index, const off_t* reverse, off_t* buckets)
static off_t find_index(off_t pos, const uint8_t* data, off_t datalen, const off_t* index, const off_t* buckets)
{
if (reverse) return reverse[pos];
//if (datalen<2) return 0;
uint16_t bucket = read2(data+pos, datalen-pos);
//printf("p=%i b=%i\n",pos,bucket);
@@ -654,29 +643,29 @@ off_t lerp(off_t x, off_t y, float frac)
}
template<typename off_t>
static bpserror bps_create_suf_core(file* source, file* target, bool moremem, struct bps_creator * out)
static result create_suf_core(const file& source, const file& target, struct bps_creator * out)
{
#define error(which) do { err = which; goto error; } while(0)
bpserror err;
result err;
size_t realsourcelen = source->len;
size_t realtargetlen = target->len;
size_t realsourcelen = source.size();
size_t realtargetlen = target.size();
size_t overflowtest = realsourcelen + realtargetlen;
//source+target length is bigger than size_t
if (overflowtest < realsourcelen) return bps_too_big;
if (overflowtest < realsourcelen) return e_too_big;
//source+target doesn't fit in unsigned off_t
if ((size_t)(off_t)overflowtest != overflowtest) return bps_too_big;
if ((size_t)(off_t)overflowtest != overflowtest) return e_too_big;
//source+target doesn't fit in signed off_t
if ((off_t)overflowtest < 0) return bps_too_big;
if ((off_t)overflowtest < 0) return e_too_big;
//the mallocs would overflow
if (realsourcelen+realtargetlen >= SIZE_MAX/sizeof(off_t)) return bps_too_big;
if (realsourcelen+realtargetlen >= SIZE_MAX/sizeof(off_t)) return e_too_big;
if (realsourcelen+realtargetlen >= out->maxsize()) return bps_too_big;
if (realsourcelen+realtargetlen >= out->maxsize()) return e_too_big;
off_t sourcelen = realsourcelen;
@@ -686,19 +675,14 @@ static bpserror bps_create_suf_core(file* source, file* target, bool moremem, st
off_t* sorted = (off_t*)malloc(sizeof(off_t)*(realsourcelen+realtargetlen));
off_t* sorted_inverse = NULL;
if (moremem) sorted_inverse = (off_t*)malloc(sizeof(off_t)*(realsourcelen+realtargetlen));
off_t* buckets = (off_t*)malloc(sizeof(off_t)*65537);
off_t* buckets = NULL;
if (!sorted_inverse) buckets = (off_t*)malloc(sizeof(off_t)*65537);
if (!sorted || !mem_joined || (!sorted_inverse && !buckets))
if (!sorted || !mem_joined || !buckets)
{
free(mem_joined);
free(sorted);
free(sorted_inverse);
free(buckets);
return bps_out_of_mem;
return e_out_of_mem;
}
//sortedsize is how much of the target file is sorted
@@ -720,35 +704,32 @@ static bpserror bps_create_suf_core(file* source, file* target, bool moremem, st
reindex:
//this isn't an exact science
const float percSort = sorted_inverse ? 0.67 : 0.50;
const float percInv = sorted_inverse ? 0.11 : 0.10;
//const float percFind = sorted_inverse ? 0.22 : 0.40; // unused
const float percSort = 0.50;
const float percBuck = 0.10;
//const float percFind = 0.40; // unused
const size_t progPreSort = lerp(prevsortedsize, sortedsize, 0);
const size_t progPreInv = lerp(prevsortedsize, sortedsize, percSort);
const size_t progPreFind = lerp(prevsortedsize, sortedsize, percSort+percInv);
const size_t progPreBuck = lerp(prevsortedsize, sortedsize, percSort);
const size_t progPreFind = lerp(prevsortedsize, sortedsize, percSort+percBuck);
prevsortedsize = sortedsize;
if (!out->progress(progPreSort, targetlen)) error(bps_canceled);
if (out->progress(progPreSort, targetlen)) error(e_canceled);
if (target->read(mem_joined, 0, sortedsize) < (size_t)sortedsize) error(bps_io);
if (source->read(mem_joined+sortedsize, 0, sourcelen) < (size_t)sourcelen) error(bps_io);
if (target.read(arrayvieww<byte>(mem_joined, sortedsize), 0) < (size_t)sortedsize) error(e_io);
if (source.read(arrayvieww<byte>(mem_joined+sortedsize, sourcelen), 0) < (size_t)sourcelen) error(e_io);
out->move_target(mem_joined);
sufsort(sorted, mem_joined, sortedsize+sourcelen);
if (!out->progress(progPreInv, targetlen)) error(bps_canceled);
if (out->progress(progPreBuck, targetlen)) error(e_canceled);
if (sorted_inverse)
create_reverse_index(sorted, sorted_inverse, sortedsize+sourcelen);
else
create_buckets(mem_joined, sorted, sortedsize+sourcelen, buckets);
create_buckets(mem_joined, sorted, sortedsize+sourcelen, buckets);
if (!out->progress(progPreFind, targetlen)) error(bps_canceled);
if (out->progress(progPreFind, targetlen)) error(e_canceled);
}
off_t matchlen = 0;
off_t matchpos = adjust_match(find_index(outpos, mem_joined, sortedsize+sourcelen, sorted, sorted_inverse, buckets),
off_t matchpos = adjust_match(find_index(outpos, mem_joined, sortedsize+sourcelen, sorted, buckets),
mem_joined+outpos, sortedsize-outpos,
mem_joined,sortedsize+sourcelen, outpos,sortedsize,
sorted, sortedsize+sourcelen,
@@ -771,11 +752,10 @@ static bpserror bps_create_suf_core(file* source, file* target, bool moremem, st
out->finish(mem_joined+sortedsize, mem_joined);
err = bps_ok;
err = e_ok;
error:
free(buckets);
free(sorted_inverse);
free(sorted);
free(mem_joined);
@@ -783,38 +763,41 @@ error:
}
template<typename T> static bpserror bps_create_suf_pick(file* source, file* target, bool moremem, struct bps_creator * bps);
template<> bpserror bps_create_suf_pick<uint32_t>(file* source, file* target, bool moremem, struct bps_creator * bps)
template<typename T> static result create_suf_pick(const file& source, const file& target, struct bps_creator * bps);
template<> result create_suf_pick<uint32_t>(const file& source, const file& target, struct bps_creator * bps)
{
return bps_create_suf_core<int32_t>(source, target, moremem, bps);
return create_suf_core<int32_t>(source, target, bps);
}
template<> bpserror bps_create_suf_pick<uint64_t>(file* source, file* target, bool moremem, struct bps_creator * bps)
template<> result create_suf_pick<uint64_t>(const file& source, const file& target, struct bps_creator * bps)
{
bpserror err = bps_create_suf_core<int32_t>(source, target, moremem, bps);
if (err==bps_too_big) err = bps_create_suf_core<int64_t>(source, target, moremem, bps);
result err = create_suf_core<int32_t>(source, target, bps);
if (err==e_too_big) err = create_suf_core<int64_t>(source, target, bps);
return err;
}
//This one picks a function based on 32-bit integers if that fits. This halves memory use for common inputs.
//It also handles some stuff related to the BPS headers and footers.
extern "C"
bpserror bps_create_delta(file* source, file* target, struct mem metadata, struct mem * patchmem,
bool (*progress)(void* userdata, size_t done, size_t total), void* userdata, bool moremem)
result create(const file& source, const file& target, const file& metadata, file& patch,
function<bool(size_t done, size_t total)> progress)
{
bps_creator bps(source, target, metadata);
bps.setProgress(progress, userdata);
mem metamem = metadata.mmap();
bps_creator bps(source, target, metamem);
metadata.unmap(metamem.v());
bps.progress = progress;
size_t maindata = bps.outlen;
//off_t must be signed
bpserror err = bps_create_suf_pick<size_t>(source, target, moremem, &bps);
if (err!=bps_ok) return err;
result err = create_suf_pick<size_t>(source, target, &bps);
if (err!=e_ok) return err;
*patchmem = bps.getpatch();
mem patchmem = bps.getpatch();
patch.write(patchmem.v());
free(patchmem.ptr);
while ((patchmem->ptr[maindata]&0x80) == 0x00) maindata++;
if (maindata==patchmem->len-12-1) return bps_identical;
return bps_ok;
while ((patchmem.ptr[maindata]&0x80) == 0x00) maindata++;
if (maindata==patchmem.len-12-1) return e_identical;
return e_ok;
}
@@ -876,3 +859,4 @@ printf("%i/%i=%f\n",match_len_tot,match_len_n,(float)match_len_tot/match_len_n);
#endif
}
#endif
}}

1
patch/bps.cpp
View File

@@ -1,6 +1,7 @@
#include "patch.h"
namespace patch { namespace bps {
//TODO: HEAVY cleanups needed here
static uint32_t read32(uint8_t * ptr)
{
uint32_t out;

63
patch/divsufsort.h Normal file
View File

@@ -0,0 +1,63 @@
/*
* divsufsort.h for libdivsufsort-lite
* 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.
*/
#ifndef _DIVSUFSORT_H
#define _DIVSUFSORT_H 1
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
/*- Prototypes -*/
/**
* Constructs the suffix array of a given string.
* @param T[0..n-1] The input string.
* @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.
*/
int
divsufsort(const unsigned char *T, int *SA, int 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-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.
*/
int
divbwt(const unsigned char *T, unsigned char *U, int *A, int n);
#ifdef __cplusplus
} /* extern "C" */
#endif /* __cplusplus */
#endif /* _DIVSUFSORT_H */

1
patch/patch.h
View File

@@ -74,6 +74,7 @@ struct info {
//Deprecated
struct mem {
mem() : ptr(NULL), len(0) {}
mem(uint8_t* ptr, size_t len) : ptr(ptr), len(len) {}
mem(arrayview<byte> v) : ptr((byte*)v.ptr()), len(v.size()) {}
arrayvieww<byte> v() { return arrayvieww<byte>(ptr, len); }
uint8_t * ptr;

454
patch/sais.cpp Normal file
View File

@@ -0,0 +1,454 @@
/*
* sais.c for sais-lite
* Copyright (c) 2008-2010 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.
*/
//This file is heavily modified from the original <https://sites.google.com/site/yuta256/sais>;
//while the algorithm is the same, many changes were done.
//- The 'cs' parameters (1 or 4, depending on whether T is int* or uint8_t*) were replaced with C++ templates. This gave a fair speedup.
//- sais_index_type was replaced with a C++ template.
//- bwt, and various other stuff I don't use, was removed.
//- Assertions were removed, as they too showed up heavily in profiles; however, I suspect that just shifted the time taken elsewhere.
#include <stdlib.h>
#include <stdint.h>
#undef assert
#define assert(x)
#ifndef MINBUCKETSIZE
# define MINBUCKETSIZE 256
#endif
#define SAIS_LMSSORT2_LIMIT 0x3fffffff
#define SAIS_MYMALLOC(_num, _type) ((_type *)malloc((_num) * sizeof(_type)))
#define SAIS_MYFREE(_ptr, _num, _type) free((_ptr))
#define chr(_a) T[_a]
/* find the start or end of each bucket */
template<typename sais_index_type, typename TT>
static
void
getCounts(const TT *T, sais_index_type *C, sais_index_type n, sais_index_type k) {
sais_index_type i;
for(i = 0; i < k; ++i) { C[i] = 0; }
for(i = 0; i < n; ++i) { ++C[chr(i)]; }
}
template<typename sais_index_type>
static
void
getBuckets(const sais_index_type *C, sais_index_type *B, sais_index_type k, bool end) {
sais_index_type i, sum = 0;
if(end) { for(i = 0; i < k; ++i) { sum += C[i]; B[i] = sum; } }
else { for(i = 0; i < k; ++i) { sum += C[i]; B[i] = sum - C[i]; } }
}
/* sort all type LMS suffixes */
template<typename sais_index_type, typename TT>
static
void
LMSsort1(const TT *T, sais_index_type *SA,
sais_index_type *C, sais_index_type *B,
sais_index_type n, sais_index_type k) {
sais_index_type *b, i, j;
sais_index_type c0, c1;
/* compute SAl */
if(C == B) { getCounts(T, C, n, k); }
getBuckets(C, B, k, false); /* find starts of buckets */
j = n - 1;
b = SA + B[c1 = chr(j)];
--j;
*b++ = (chr(j) < c1) ? ~j : j;
for(i = 0; i < n; ++i) {
if(0 < (j = SA[i])) {
assert(chr(j) >= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert(i < (b - SA));
--j;
*b++ = (chr(j) < c1) ? ~j : j;
SA[i] = 0;
} else if(j < 0) {
SA[i] = ~j;
}
}
/* compute SAs */
if(C == B) { getCounts(T, C, n, k); }
getBuckets(C, B, k, true); /* find ends of buckets */
for(i = n - 1, b = SA + B[c1 = 0]; 0 <= i; --i) {
if(0 < (j = SA[i])) {
assert(chr(j) <= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert((b - SA) <= i);
--j;
*--b = (chr(j) > c1) ? ~(j + 1) : j;
SA[i] = 0;
}
}
}
template<typename sais_index_type, typename TT>
static
sais_index_type
LMSpostproc1(const TT *T, sais_index_type *SA,
sais_index_type n, sais_index_type m) {
sais_index_type i, j, p, q, plen, qlen, name;
sais_index_type c0, c1;
bool diff;
/* compact all the sorted substrings into the first m items of SA
2*m must be not larger than n (proveable) */
assert(0 < n);
for(i = 0; (p = SA[i]) < 0; ++i) { SA[i] = ~p; assert((i + 1) < n); }
if(i < m) {
for(j = i, ++i;; ++i) {
assert(i < n);
if((p = SA[i]) < 0) {
SA[j++] = ~p; SA[i] = 0;
if(j == m) { break; }
}
}
}
/* store the length of all substrings */
i = n - 1; j = n - 1; c0 = chr(n - 1);
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
for(; 0 <= i;) {
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) <= c1));
if(0 <= i) {
SA[m + ((i + 1) >> 1)] = j - i; j = i + 1;
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
}
}
/* find the lexicographic names of all substrings */
for(i = 0, name = 0, q = n, qlen = 0; i < m; ++i) {
p = SA[i], plen = SA[m + (p >> 1)], diff = true;
if((plen == qlen) && ((q + plen) < n)) {
for(j = 0; (j < plen) && (chr(p + j) == chr(q + j)); ++j) { }
if(j == plen) { diff = false; }
}
if(diff) { ++name, q = p, qlen = plen; }
SA[m + (p >> 1)] = name;
}
return name;
}
template<typename sais_index_type, typename TT>
static
void
LMSsort2(const TT *T, sais_index_type *SA,
sais_index_type *C, sais_index_type *B, sais_index_type *D,
sais_index_type n, sais_index_type k) {
sais_index_type *b, i, j, t, d;
sais_index_type c0, c1;
assert(C != B);
/* compute SAl */
getBuckets(C, B, k, false); /* find starts of buckets */
j = n - 1;
b = SA + B[c1 = chr(j)];
--j;
t = (chr(j) < c1);
j += n;
*b++ = (t & 1) ? ~j : j;
for(i = 0, d = 0; i < n; ++i) {
if(0 < (j = SA[i])) {
if(n <= j) { d += 1; j -= n; }
assert(chr(j) >= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert(i < (b - SA));
--j;
t = c0; t = (t << 1) | (chr(j) < c1);
if(D[t] != d) { j += n; D[t] = d; }
*b++ = (t & 1) ? ~j : j;
SA[i] = 0;
} else if(j < 0) {
SA[i] = ~j;
}
}
for(i = n - 1; 0 <= i; --i) {
if(0 < SA[i]) {
if(SA[i] < n) {
SA[i] += n;
for(j = i - 1; SA[j] < n; --j) { }
SA[j] -= n;
i = j;
}
}
}
/* compute SAs */
getBuckets(C, B, k, true); /* find ends of buckets */
for(i = n - 1, d += 1, b = SA + B[c1 = 0]; 0 <= i; --i) {
if(0 < (j = SA[i])) {
if(n <= j) { d += 1; j -= n; }
assert(chr(j) <= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert((b - SA) <= i);
--j;
t = c0; t = (t << 1) | (chr(j) > c1);
if(D[t] != d) { j += n; D[t] = d; }
*--b = (t & 1) ? ~(j + 1) : j;
SA[i] = 0;
}
}
}
template<typename sais_index_type>
static
sais_index_type
LMSpostproc2(sais_index_type *SA, sais_index_type n, sais_index_type m) {
sais_index_type i, j, d, name;
/* compact all the sorted LMS substrings into the first m items of SA */
assert(0 < n);
for(i = 0, name = 0; (j = SA[i]) < 0; ++i) {
j = ~j;
if(n <= j) { name += 1; }
SA[i] = j;
assert((i + 1) < n);
}
if(i < m) {
for(d = i, ++i;; ++i) {
assert(i < n);
if((j = SA[i]) < 0) {
j = ~j;
if(n <= j) { name += 1; }
SA[d++] = j; SA[i] = 0;
if(d == m) { break; }
}
}
}
if(name < m) {
/* store the lexicographic names */
for(i = m - 1, d = name + 1; 0 <= i; --i) {
if(n <= (j = SA[i])) { j -= n; --d; }
SA[m + (j >> 1)] = d;
}
} else {
/* unset flags */
for(i = 0; i < m; ++i) {
if(n <= (j = SA[i])) { j -= n; SA[i] = j; }
}
}
return name;
}
/* compute SA and BWT */
template<typename sais_index_type, typename TT>
static
void
induceSA(const TT *T, sais_index_type *SA,
sais_index_type *C, sais_index_type *B,
sais_index_type n, sais_index_type k) {
sais_index_type *b, i, j;
sais_index_type c0, c1;
/* compute SAl */
if(C == B) { getCounts(T, C, n, k); }
getBuckets(C, B, k, false); /* find starts of buckets */
j = n - 1;
b = SA + B[c1 = chr(j)];
*b++ = ((0 < j) && (chr(j - 1) < c1)) ? ~j : j;
for(i = 0; i < n; ++i) {
j = SA[i], SA[i] = ~j;
if(0 < j) {
--j;
assert(chr(j) >= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert(i < (b - SA));
*b++ = ((0 < j) && (chr(j - 1) < c1)) ? ~j : j;
}
}
/* compute SAs */
if(C == B) { getCounts(T, C, n, k); }
getBuckets(C, B, k, true); /* find ends of buckets */
for(i = n - 1, b = SA + B[c1 = 0]; 0 <= i; --i) {
if(0 < (j = SA[i])) {
--j;
assert(chr(j) <= chr(j + 1));
if((c0 = chr(j)) != c1) { B[c1] = b - SA; b = SA + B[c1 = c0]; }
assert((b - SA) <= i);
*--b = ((j == 0) || (chr(j - 1) > c1)) ? ~j : j;
} else {
SA[i] = ~j;
}
}
}
/* find the suffix array SA of T[0..n-1] in {0..255}^n */
template<typename sais_index_type, typename TT>
static
sais_index_type
sais_main(const TT *T, sais_index_type *SA,
sais_index_type fs, sais_index_type n, sais_index_type k) {
sais_index_type *C, *B, *D, *RA, *b;
sais_index_type i, j, m, p, q, t, name, pidx = 0, newfs;
sais_index_type c0, c1;
unsigned int flags;
assert((T != NULL) && (SA != NULL));
assert((0 <= fs) && (0 < n) && (1 <= k));
if(k <= MINBUCKETSIZE) {
if((C = SAIS_MYMALLOC(k, sais_index_type)) == NULL) { return -2; }
if(k <= fs) {
B = SA + (n + fs - k);
flags = 1;
} else {
if((B = SAIS_MYMALLOC(k, sais_index_type)) == NULL) { SAIS_MYFREE(C, k, sais_index_type); return -2; }
flags = 3;
}
} else if(k <= fs) {
C = SA + (n + fs - k);
if(k <= (fs - k)) {
B = C - k;
flags = 0;
} else if(k <= (MINBUCKETSIZE * 4)) {
if((B = SAIS_MYMALLOC(k, sais_index_type)) == NULL) { return -2; }
flags = 2;
} else {
B = C;
flags = 8;
}
} else {
if((C = B = SAIS_MYMALLOC(k, sais_index_type)) == NULL) { return -2; }
flags = 4 | 8;
}
if((n <= SAIS_LMSSORT2_LIMIT) && (2 <= (n / k))) {
if(flags & 1) { flags |= ((k * 2) <= (fs - k)) ? 32 : 16; }
else if((flags == 0) && ((k * 2) <= (fs - k * 2))) { flags |= 32; }
}
/* stage 1: reduce the problem by at least 1/2
sort all the LMS-substrings */
getCounts(T, C, n, k); getBuckets(C, B, k, true); /* find ends of buckets */
for(i = 0; i < n; ++i) { SA[i] = 0; }
b = &t; i = n - 1; j = n; m = 0; c0 = chr(n - 1);
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
for(; 0 <= i;) {
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) <= c1));
if(0 <= i) {
*b = j; b = SA + --B[c1]; j = i; ++m;
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
}
}
if(1 < m) {
if(flags & (16 | 32)) {
if(flags & 16) {
if((D = SAIS_MYMALLOC(k * 2, sais_index_type)) == NULL) {
if(flags & (1 | 4)) { SAIS_MYFREE(C, k, sais_index_type); }
if(flags & 2) { SAIS_MYFREE(B, k, sais_index_type); }
return -2;
}
} else {
D = B - k * 2;
}
assert((j + 1) < n);
++B[chr(j + 1)];
for(i = 0, j = 0; i < k; ++i) {
j += C[i];
if(B[i] != j) { assert(SA[B[i]] != 0); SA[B[i]] += n; }
D[i] = D[i + k] = 0;
}
LMSsort2(T, SA, C, B, D, n, k);
name = LMSpostproc2(SA, n, m);
if(flags & 16) { SAIS_MYFREE(D, k * 2, sais_index_type); }
} else {
LMSsort1(T, SA, C, B, n, k);
name = LMSpostproc1(T, SA, n, m);
}
} else if(m == 1) {
*b = j + 1;
name = 1;
} else {
name = 0;
}
/* stage 2: solve the reduced problem
recurse if names are not yet unique */
if(name < m) {
if(flags & 4) { SAIS_MYFREE(C, k, sais_index_type); }
if(flags & 2) { SAIS_MYFREE(B, k, sais_index_type); }
newfs = (n + fs) - (m * 2);
if((flags & (1 | 4 | 8)) == 0) {
if((k + name) <= newfs) { newfs -= k; }
else { flags |= 8; }
}
assert((n >> 1) <= (newfs + m));
RA = SA + m + newfs;
for(i = m + (n >> 1) - 1, j = m - 1; m <= i; --i) {
if(SA[i] != 0) {
RA[j--] = SA[i] - 1;
}
}
if(sais_main(RA, SA, newfs, m, name) != 0) {
if(flags & 1) { SAIS_MYFREE(C, k, sais_index_type); }
return -2;
}
i = n - 1; j = m - 1; c0 = chr(n - 1);
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
for(; 0 <= i;) {
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) <= c1));
if(0 <= i) {
RA[j--] = i + 1;
do { c1 = c0; } while((0 <= --i) && ((c0 = chr(i)) >= c1));
}
}
for(i = 0; i < m; ++i) { SA[i] = RA[SA[i]]; }
if(flags & 4) {
if((C = B = SAIS_MYMALLOC(k, sais_index_type)) == NULL) { return -2; }
}
if(flags & 2) {
if((B = SAIS_MYMALLOC(k, sais_index_type)) == NULL) {
if(flags & 1) { SAIS_MYFREE(C, k, sais_index_type); }
return -2;
}
}
}
/* stage 3: induce the result for the original problem */
if(flags & 8) { getCounts(T, C, n, k); }
/* put all left-most S characters into their buckets */
if(1 < m) {
getBuckets(C, B, k, true); /* find ends of buckets */
i = m - 1, j = n, p = SA[m - 1], c1 = chr(p);
do {
q = B[c0 = c1];
while(q < j) { SA[--j] = 0; }
do {
SA[--j] = p;
if(--i < 0) { break; }
p = SA[i];
} while((c1 = chr(p)) == c0);
} while(0 <= i);
while(0 < j) { SA[--j] = 0; }
}
induceSA(T, SA, C, B, n, k);
if(flags & 2) { SAIS_MYFREE(B, k, sais_index_type); }
if(flags & (1 | 4)) { SAIS_MYFREE(C, k, sais_index_type); }
return pidx;
}