Code:

/*
* File: main.cpp
* Author: Piotr Tarsa
*
*/
#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <queue>
#include <stack>
#include <utility>
#include <vector>
//#define USE_CONTEXT_0_NOT_1
struct _index_t {
int32_t index;
int32_t info;
};
typedef struct _index_t index_t;
const int32_t decodingLookupBits = 10;
struct _context_t {
u_int8_t index[256];
u_int8_t value[256];
u_int8_t logRange[256];
u_int32_t cumulativeRange[256];
u_int8_t decodedIndex[1 << decodingLookupBits];
int32_t numSymbols;
};
typedef struct _context_t context_t;
struct _huffman_node_t {
int32_t weight;
int32_t index;
int32_t leftChild;
int32_t rightChild;
int operator()(const _huffman_node_t &_this, const _huffman_node_t & _that) const {
return _this.weight > _that.weight;
}
};
typedef struct _huffman_node_t huffman_node_t;
size_t length;
int32_t *SA;
u_int8_t *sourceBuffer;
u_int8_t *outputBuffer;
template<typename T>
T min(T former, T latter) {
return former < latter ? former : latter;
}
template<typename T>
T max(T former, T latter) {
return former > latter ? former : latter;
}
int32_t clampOnce(int32_t value, int32_t limit) {
return (value >= limit) ? value - limit : value;
}
/*
*
*/
int main(int argc, char** argv) {
freopen("report.txt", "w", stderr);
//FILE *sourceFile = fopen("/home/piotrek/Dokumenty/BWT/qsufsort.c", "rb");
FILE *sourceFile = fopen("/home/piotrek/Pobrane/corpora/enwik/enwik8", "rb");
fseek(sourceFile, 0, SEEK_END);
length = ftell(sourceFile);
const int32_t DepthLimit = length;
printf("File length: %ld\n", length);
fseek(sourceFile, 0, SEEK_SET);
SA = new int32_t[length];
sourceBuffer = new u_int8_t[length];
outputBuffer = new u_int8_t[length];
size_t readBytes = fread(sourceBuffer, sizeof (u_int8_t), length, sourceFile);
assert(readBytes == length);
fclose(sourceFile);
clock_t time = clock();
int32_t counts[256][256];
for (int32_t i = 0; i < 256; i++) {
for (int32_t j = 0; j < 256; j++) {
counts[i][j] = 0;
}
}
{
u_int8_t lastChar = sourceBuffer[length - 1];
for (u_int32_t i = 0; i < length; i++) {
#ifdef USE_CONTEXT_0_NOT_1
lastChar = 0;
#endif
const u_int8_t currentChar = sourceBuffer[i];
counts[lastChar][currentChar]++;
lastChar = currentChar;
}
}
context_t contexts[256];
{
huffman_node_t nodes[512];
std::priority_queue<huffman_node_t, std::vector<huffman_node_t>, huffman_node_t> heap;
for (int32_t context = 0; context < 256; context++) {
bool needRescaling = false;
do {
int32_t numSymbols = 0;
for (int32_t i = 0; i < 256; i++) {
if (counts[context][i] > 0) {
huffman_node_t huffmanNode;
huffmanNode.weight = counts[context][i];
contexts[context].value[numSymbols] = i;
contexts[context].index[i] = numSymbols;
huffmanNode.index = numSymbols++;
heap.push(huffmanNode);
}
}
contexts[context].numSymbols = numSymbols;
if (numSymbols == 0) {
continue;
}
int32_t firstFree = 0;
int32_t rootNode = 0;
while (true) {
huffman_node_t firstNode = heap.top();
heap.pop();
if (heap.empty()) {
rootNode = firstFree;
nodes[firstFree++] = firstNode;
break;
}
huffman_node_t secondNode = heap.top();
heap.pop();
huffman_node_t newNode;
newNode.weight = firstNode.weight + secondNode.weight;
newNode.index = -1;
newNode.leftChild = firstFree;
nodes[firstFree++] = firstNode;
newNode.rightChild = firstFree;
nodes[firstFree++] = secondNode;
heap.push(newNode);
}
int32_t currentLogRange = 32;
std::stack<std::pair<int32_t, int32_t> > stack;
int32_t currentNode = rootNode;
while (true) {
if (nodes[currentNode].index == -1) {
currentLogRange--;
needRescaling |= (currentLogRange < 0);
stack.push(std::make_pair(nodes[currentNode].rightChild, currentLogRange));
currentNode = nodes[currentNode].leftChild;
} else {
contexts[context].logRange[nodes[currentNode].index] = currentLogRange;
if (stack.empty()) {
break;
} else {
currentNode = stack.top().first;
currentLogRange = stack.top().second;
stack.pop();
}
}
}
u_int32_t currentLow = 0;
for (int32_t i = 0; i < numSymbols; i++) {
contexts[context].cumulativeRange[i] = currentLow;
currentLow += 1 << contexts[context].logRange[i];
}
u_int8_t decodedIndex = numSymbols - 1;
for (int32_t i = (1 << decodingLookupBits) - 1; i >= 0; i--) {
u_int32_t currentLow = i << (32 - decodingLookupBits);
while (contexts[context].cumulativeRange[decodedIndex] > currentLow) {
decodedIndex--;
}
contexts[context].decodedIndex[i] = decodedIndex;
}
if (needRescaling) {
for (u_int32_t symbol = 0; symbol < 256; symbol++) {
counts[context][symbol] += 1;
counts[context][symbol] /= 2;
}
puts("rescaling");
}
} while (needRescaling);
}
}
bool useCompressedInput = true;
{
int64_t totalCost = 0;
for (int32_t context = 0; context < 256; context++) {
for (int32_t symbol = 0; symbol < 256; symbol++) {
if (counts[context][symbol] > 0) {
totalCost += counts[context][symbol] *
(32 - contexts[context].logRange[contexts[context].index[symbol]]);
}
}
}
double ratio = (double) totalCost / 8. / length;
printf("Compression ratio: %lf\n", ratio);
if ((length < 1000) || (ratio > 0.9)) {
useCompressedInput = false;
}
}
printf("Process time: %ld\n", clock() - time);
// compute compressed input
if (useCompressedInput) {
time = clock();
for (u_int32_t i = 0; i < length; i++) {
outputBuffer[i] = 0;
}
int32_t outputIndex = 0;
int32_t numFFdwords = 0;
int32_t numHeadBits = 0;
int32_t numTailBits = 0;
u_int32_t headDword = 0;
u_int32_t tailDword = 0;
u_int32_t currentLow = 0;
int32_t lastByte = sourceBuffer[length - 1];
for (u_int32_t i = 0; i <= length; i++) {
#ifdef USE_CONTEXT_0_NOT_1
lastByte = 0;
#endif
int32_t currentByte = (i == length) ? 0 : sourceBuffer[i];
int32_t currentIndex = (i == length) ? 0 : contexts[lastByte].index[currentByte];
int32_t currentLogRange = (i == length) ? 0 : contexts[lastByte].logRange[currentIndex];
int32_t bitsToShift = (i == length) ? 32 : (32 - currentLogRange);
u_int32_t newLow = currentLow + ((i == length) ? *(u_int32_t*) (sourceBuffer) :
contexts[lastByte].cumulativeRange[currentIndex]);
if (newLow < currentLow) { // overflow
if (numHeadBits == 32) {
if (numTailBits > 0) {
tailDword += 1 << (32 - numTailBits);
}
if (tailDword == 0) {
*(u_int32_t*) (outputBuffer + outputIndex) = headDword + 1;
outputIndex += 4;
while (numFFdwords > 0) {
numFFdwords--;
*(u_int32_t*) (outputBuffer + outputIndex) = 0;
outputIndex += 4;
}
headDword = 0;
numHeadBits = numTailBits;
numTailBits = 0;
} else if (tailDword == 0xffffffff) {
numFFdwords++;
tailDword = 0;
numTailBits = 0;
}
} else {
if (numHeadBits > 0) {
headDword += 1 << (32 - numHeadBits);
}
}
}
if (numHeadBits == 32) {
if (numTailBits + bitsToShift < 32) {
if (bitsToShift > 0) {
tailDword += (newLow >> currentLogRange) << (currentLogRange - numTailBits);
numTailBits += bitsToShift;
}
} else {
int32_t tailBitsToShift = 32 - numTailBits;
tailDword += newLow >> numTailBits;
if (tailDword == 0xffffffff) {
numFFdwords++;
} else {
*(u_int32_t*) (outputBuffer + outputIndex) = headDword;
outputIndex += 4;
while (numFFdwords > 0) {
numFFdwords--;
*(u_int32_t*) (outputBuffer + outputIndex) = 0xffffffff;
outputIndex += 4;
}
headDword = tailDword;
numHeadBits = 32;
}
tailDword = (currentLogRange + tailBitsToShift < 32) ?
(newLow >> currentLogRange) << (currentLogRange + tailBitsToShift) : 0;
numTailBits = bitsToShift - tailBitsToShift;
}
} else {
if (numHeadBits + bitsToShift <= 32) {
headDword += (newLow >> currentLogRange) << (currentLogRange - numHeadBits);
numHeadBits += bitsToShift;
} else {
int32_t headBitsToShift = 32 - numHeadBits;
headDword += newLow >> numHeadBits;
tailDword = (newLow >> currentLogRange) << (currentLogRange + headBitsToShift);
numHeadBits = 32;
numTailBits = bitsToShift - headBitsToShift;
}
}
currentLow = newLow << bitsToShift;
lastByte = currentByte;
}
if (numHeadBits > 0) {
*(u_int32_t*) (outputBuffer + outputIndex) = headDword;
outputIndex += 4;
while (numFFdwords > 0) {
numFFdwords--;
*(u_int32_t*) (outputBuffer + outputIndex) = 0xffffffff;
outputIndex += 4;
}
if (numTailBits > 0) {
*(u_int32_t*) (outputBuffer + outputIndex) = tailDword;
outputIndex += 4;
}
}
printf("Output index: %d\n", outputIndex);
printf("Process time: %ld\n", clock() - time);
}
if (true) {
// test print context
int32_t testedContext = sourceBuffer[length - 1];
context_t c = contexts[testedContext];
for (int32_t index = 0; index < c.numSymbols; index++) {
printf("%3d %3d %2d %8x\n", index, c.value[index], c.logRange[index],
c.cumulativeRange[index]);
}
puts("");
for (int32_t symbol = 0; symbol < 256; symbol++) {
if (counts[testedContext][symbol] > 0) {
printf("%3d %3d\n", symbol, c.index[symbol]);
}
}
for (int32_t low = 0; low < 1024; low++) {
printf("%3x %3d\n", low << 2, c.decodedIndex[low]);
}
}
// verify compressed input
if (useCompressedInput) {
bool foundBrokenIndex = false;
time = clock();
u_int32_t currentLow = *(u_int32_t*) (outputBuffer + 0);
u_int32_t bufferDword = *(u_int32_t*) (outputBuffer + 4);
int32_t bufferBits = 32;
int32_t outputIndex = 8;
int32_t lastByte = sourceBuffer[length - 1];
for (u_int32_t i = 0; i < length; i++) {
#ifdef USE_CONTEXT_0_NOT_1
lastByte = 0;
#endif
context_t ¤tContext = contexts[lastByte];
int32_t currentIndex =
currentContext.decodedIndex[currentLow >> (32 - decodingLookupBits)];
while (((currentIndex + 1) < currentContext.numSymbols) &&
(currentLow >= currentContext.cumulativeRange[currentIndex + 1])) {
currentIndex++;
}
lastByte = currentContext.value[currentIndex];
if ((lastByte != sourceBuffer[i]) && (!foundBrokenIndex)) {
foundBrokenIndex = true;
printf("First broken index: %d\n", i);
}
int32_t bitsToShift = 32 - currentContext.logRange[currentIndex];
if (bitsToShift != 0) {
currentLow -= currentContext.cumulativeRange[currentIndex];
currentLow <<= bitsToShift;
currentLow += bufferDword >> (32 - bitsToShift);
if (bitsToShift < bufferBits) {
bufferDword <<= bitsToShift;
bufferBits -= bitsToShift;
} else if (bitsToShift > bufferBits) {
bufferDword = *(u_int32_t*) (outputBuffer + outputIndex);
outputIndex += 4;
bitsToShift -= bufferBits;
currentLow += bufferDword >> (32 - bitsToShift);
bufferDword <<= bitsToShift;
bufferBits = 32 - bitsToShift;
} else {
bufferDword = *(u_int32_t*) (outputBuffer + outputIndex);
outputIndex += 4;
bitsToShift -= bufferBits;
bufferBits = 32;
}
}
}
printf("Process time: %ld\n", clock() - time);
length = outputIndex;
}
fflush(stdout);
return 0;
for (u_int32_t i = 0; i < length; i++) {
SA[i] = i;
}
for (u_int32_t i = 0; i < length; i++) {
outputBuffer[i] = sourceBuffer[clampOnce(SA[i] + length - 1, length)];
}
{
clock_t time = clock();
for (u_int32_t i = 1; i < length; i++) {
int32_t leftPtr = SA[i - 1];
int32_t middlePtr = SA[i];
int32_t equal = 0;
for (; equal < DepthLimit; equal++) {
if (sourceBuffer[leftPtr] > sourceBuffer[middlePtr]) {
printf("Wrong order on position: %d, equal prefix length: %d\n", i - 1, equal);
break;
} else if (sourceBuffer[leftPtr] < sourceBuffer[middlePtr]) {
break;
}
leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
}
}
printf("Rotation order checking time: %ld\n", clock() - time);
}
FILE *outputFile = fopen("/home/piotrek/Pulpit/bwtoutput.mergerotsort", "wb");
size_t writtenBytes = fwrite(outputBuffer, sizeof (u_int8_t), length, outputFile);
assert(writtenBytes == length);
fclose(outputFile);
return 0;
}