MergeSort based blocksorting algorithm

[Piotr Tarsa]

Here's my newly invented (or not, probably someone did such thing before, as it's really simple) blocksorting algorithm, which I'll probably use as building block in my complex QRotSort algorithm.

I've tried stuffing suffix.lcp'th byte from each suffix into suffix.lcp, so then sometimes only comparing those bytes could suffice to figure out the suffix order and lcps, but in the end the gain was really small (and complexity is higher) and probably some other technique will bring some improvement.

Currently I'm thinking about multi-way MergeSort, but that will be tricky. Also the number of ways is yet to be decided.

/* 
 * File:   main.cpp
 * Author: Piotr Tarsa
 *
 */

#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <ctime>

#include <algorithm>

#define COMPUTE_STATISTICS

struct _index_t {
    int32_t index;
    int32_t lcp;
};

typedef _index_t index_t;

size_t length;
index_t *SAin;
index_t *SAout;
u_int8_t *sourceBuffer;
u_int8_t *outputBuffer;

int32_t min(int32_t former, int32_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 = min(8000000, length);
    printf("File length: %ld\n", length);
    fseek(sourceFile, 0, SEEK_SET);

    SAin = new index_t[length];
    SAout = new index_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);

    for (int32_t i = 0; i < length; i++) {
        SAin[i].index = i;
        SAin[i].lcp = 0;
    }


#ifdef COMPUTE_STATISTICS
    int64_t comparisons = 0;
    int64_t quickChecks = 0;
#endif

    int32_t sortedLength = 1;
    while (sortedLength < length) {

        int32_t leftStart = 0;
        int32_t leftEnd = sortedLength;
        int32_t rightStart = sortedLength;
        int32_t rightEnd = min(length, sortedLength * 2);
        int32_t outputIndex = 0;

        do {
            int32_t toLeftLcp = 0;
            int32_t toRightLcp = 0;
            int32_t lastIndex = 0;

            {
                int32_t firstPtr = SAin[leftStart].index;
                int32_t secondPtr = SAin[rightStart].index;
                int32_t equal = 0;
                for (; equal < DepthLimit; equal++) {
                    if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                        break;
                    }
                    firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                    secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                }
#ifdef COMPUTE_STATISTICS
                comparisons++;
#endif
                if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
                    toRightLcp = equal;
                } else {
                    lastIndex = SAin[rightStart].index;
                    toRightLcp = SAin[rightStart].lcp;
                    rightStart++;
                    toLeftLcp = equal;
                }
            }

            while ((leftStart < leftEnd) && (rightStart < rightEnd)) {

                SAout[outputIndex].index = lastIndex;

                if (toLeftLcp > toRightLcp) {
                    SAout[outputIndex++].lcp = toLeftLcp;
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else if (toLeftLcp < toRightLcp) {
                    SAout[outputIndex++].lcp = toRightLcp;
                    lastIndex = SAin[rightStart].index;
                    toRightLcp = SAin[rightStart].lcp;
                    rightStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else { // toLeftLcp == toRightLcp
                    int32_t firstPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                    int32_t secondPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                    int32_t equal = toLeftLcp;
                    for (; equal < DepthLimit; equal++) {
                        if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                            break;
                        }
                        firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                        secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                    }
#ifdef COMPUTE_STATISTICS
                    comparisons++;
#endif
                    if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                        SAout[outputIndex++].lcp = toLeftLcp;
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
                        toRightLcp = equal;
                    } else {
                        SAout[outputIndex++].lcp = toRightLcp;
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
                        toLeftLcp = equal;
                    }
                }
            }

            if (leftStart == leftEnd) {
                SAout[outputIndex].index = lastIndex;
                SAout[outputIndex].lcp = toRightLcp;
                outputIndex++;
#ifdef COMPUTE_STATISTICS
                quickChecks++;
#endif
                while (rightStart < rightEnd) {
                    SAout[outputIndex++] = SAin[rightStart++];
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                }
            } else { // rightStart == rightEnd
                SAout[outputIndex].index = lastIndex;
                SAout[outputIndex].lcp = toLeftLcp;
                outputIndex++;
#ifdef COMPUTE_STATISTICS
                quickChecks++;
#endif
                while (leftStart < leftEnd) {
                    SAout[outputIndex++] = SAin[leftStart++];
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                }
            }

            leftStart = rightEnd;
            leftEnd = min(leftStart + sortedLength, length);
            rightStart = leftEnd;
            rightEnd = min(rightStart + sortedLength, length);
        } while (rightStart < rightEnd);

        while (leftStart < leftEnd) {
            SAout[outputIndex++] = SAin[leftStart++];
#ifdef COMPUTE_STATISTICS
            quickChecks++;
#endif
        }

        {
            index_t *temp = SAin;
            SAin = SAout;
            SAout = temp;
        }
        sortedLength *= 2;
    }

#ifdef COMPUTE_STATISTICS
    printf("Comparisons:           %11ld\n", comparisons);
    printf("Quick checks:          %11ld\n", quickChecks);
#endif


    for (int32_t i = 0; i < length; i++) {
        outputBuffer[i] = sourceBuffer[clampOnce(SAin[i].index + length - 1, length)];
    }

    {
        clock_t time = clock();
        for (int32_t i = 1; i < length; i++) {
            int32_t firstPtr = SAin[i - 1].index;
            int32_t secondPtr = SAin[i].index;
            int32_t equal = 0;
            for (; equal < DepthLimit; equal++) {
                if (sourceBuffer[firstPtr] > sourceBuffer[secondPtr]) {
                    printf("Wrong order on position: %d, equal prefix length: %d\n", i - 1, equal);
                    break;
                } else if (sourceBuffer[firstPtr] < sourceBuffer[secondPtr]) {
                    break;
                }
                firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 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;
}

[Piotr Tarsa]

OK. I've managed to make three-way version. Unfortunately, it's only little faster, like 10 % on enwik8. But it's always something.

/* 
 * File:   main.cpp
 * Author: Piotr Tarsa
 *
 */

#include <cassert>
#include <cstdio>
#include <cstdlib>
#include <ctime>

#include <algorithm>

#define THREE_WAY_MERGE
#define INITIAL_SORTED_LENGTH 5

#define COMPUTE_STATISTICS

struct _index_t {
    u_int32_t index;
    u_int32_t lcp;
};

typedef _index_t index_t;

size_t length;
index_t *SAin;
index_t *SAout;
u_int8_t *sourceBuffer;
u_int8_t *outputBuffer;

int32_t min(int32_t former, int32_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);

    SAin = new index_t[length];
    SAout = new index_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);

    for (u_int32_t i = 0; i < length; i++) {
        SAin[i].index = i;
        SAin[i].lcp = 0;
    }
    
#ifdef COMPUTE_STATISTICS
    int64_t comparisons = 0;
    int64_t quickChecks = 0;
#endif

#if INITIAL_SORTED_LENGTH > 1
    for (u_int32_t currentStart = 0; currentStart < length; currentStart += INITIAL_SORTED_LENGTH) {
        u_int32_t currentEnd = currentStart + INITIAL_SORTED_LENGTH - 1;
        if (currentEnd >= length) {
            currentEnd = length - 1;
        }
        for (u_int32_t current = currentStart + 1; current <= currentEnd; current++) {
            int32_t currentIndex = SAin[current].index;
            u_int32_t currentLcp = 0;
            int32_t firstPtr = currentIndex;
            u_int8_t firstValue = sourceBuffer[firstPtr];
            u_int32_t j = current;
            do {
#ifdef COMPUTE_STATISTICS
                comparisons++;
#endif
                u_int32_t newLcp = currentLcp;
                int32_t secondPtr = clampOnce(SAin[j - 1].index + currentLcp, length);
                u_int8_t secondValue = sourceBuffer[secondPtr];
                while ((firstValue == secondValue) && (newLcp < length)) {
                    newLcp++;
                    firstPtr = (firstPtr + 1) == length ? 0 : firstPtr + 1;
                    firstValue = sourceBuffer[firstPtr];
                    secondPtr = (secondPtr + 1) == length ? 0 : secondPtr + 1;
                    secondValue = sourceBuffer[secondPtr];
                }
                if (secondValue <= firstValue) {
                    SAin[j - 1].lcp = newLcp;
                    break;
                }
                currentLcp = newLcp;
                do {
                    SAin[j] = SAin[j - 1];
                    j--;
                } while ((j > currentStart) && (currentLcp < SAin[j - 1].lcp));
            } while ((j > currentStart) && (currentLcp <= SAin[j - 1].lcp));
            SAin[j].index = currentIndex;
            SAin[j].lcp = currentLcp;
        }
    }
#endif

    u_int32_t sortedLength = INITIAL_SORTED_LENGTH;
    while (sortedLength < length) {

        clock_t time = clock();
#ifdef COMPUTE_STATISTICS
        int64_t lastComparisons = comparisons;
        int64_t lastQuickChecks = quickChecks;
#endif

        u_int32_t first = 0;
        int32_t outputIndex = 0;

        do {
            int32_t leftStart = first;
            int32_t leftEnd = min(length, leftStart + sortedLength);
            int32_t middleStart = leftEnd;
            int32_t middleEnd = min(length, middleStart + sortedLength);
            int32_t rightStart = middleEnd;
#ifdef THREE_WAY_MERGE
            int32_t rightEnd = min(length, rightStart + sortedLength);
#else
            int32_t rightEnd = rightStart;
#endif

            int32_t toLeftLcp = 0;
            int32_t toMiddleLcp = 0;
            int32_t toRightLcp = 0;
            int32_t lastIndex = 0;

            if ((leftStart < leftEnd) && (middleStart < middleEnd) && (rightStart < rightEnd)) {
                int32_t leftPtr = SAin[leftStart].index;
                int32_t middlePtr = SAin[middleStart].index;
                int32_t rightPtr = SAin[rightStart].index;
                u_int8_t leftValue = sourceBuffer[leftPtr];
                u_int8_t middleValue = sourceBuffer[middlePtr];
                u_int8_t rightValue = sourceBuffer[rightPtr];
                int32_t equal = 0;
                while (true) {
                    if (leftValue != middleValue) {
                        if (leftValue < middleValue) {
                            if (leftValue < rightValue) {
                                lastIndex = SAin[leftStart].index;
                                toLeftLcp = SAin[leftStart].lcp;
                                leftStart++;
                                toMiddleLcp = equal;
                                toRightLcp = equal;
                            } else if (leftValue > rightValue) {
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toLeftLcp = equal;
                                toMiddleLcp = equal;
                            } else {
                                toMiddleLcp = equal;
                                if (equal < DepthLimit) {
                                    do {
                                        equal++;
                                        leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                        leftValue = sourceBuffer[leftPtr];
                                        rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                        rightValue = sourceBuffer[rightPtr];
                                    } while ((leftValue == rightValue) && (equal < DepthLimit));
                                }
                                if (leftValue <= rightValue) {
                                    lastIndex = SAin[leftStart].index;
                                    toLeftLcp = SAin[leftStart].lcp;
                                    leftStart++;
                                    toRightLcp = equal;
                                } else {
                                    lastIndex = SAin[rightStart].index;
                                    toRightLcp = SAin[rightStart].lcp;
                                    rightStart++;
                                    toLeftLcp = equal;
                                }
                            }
                        } else if (middleValue > rightValue) {
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toLeftLcp = equal;
                            toMiddleLcp = equal;
                        } else if (middleValue < rightValue) {
                            lastIndex = SAin[middleStart].index;
                            toMiddleLcp = SAin[middleStart].lcp;
                            middleStart++;
                            toLeftLcp = equal;
                            toRightLcp = equal;
                        } else {
                            toLeftLcp = equal;
                            if (equal < DepthLimit) {
                                do {
                                    equal++;
                                    middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                    middleValue = sourceBuffer[middlePtr];
                                    rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                    rightValue = sourceBuffer[rightPtr];
                                } while ((middleValue == rightValue) && (equal < DepthLimit));
                            }
                            if (middleValue <= rightValue) {
                                lastIndex = SAin[middleStart].index;
                                toMiddleLcp = SAin[middleStart].lcp;
                                middleStart++;
                                toRightLcp = equal;
                            } else {
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toMiddleLcp = equal;
                            }
                        }
                        break;
                    } else if (leftValue != rightValue) {
                        if (leftValue > rightValue) {
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toLeftLcp = equal;
                            toMiddleLcp = equal;
                        } else {
                            toRightLcp = equal;
                            if (equal < DepthLimit) {
                                do {
                                    equal++;
                                    leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                    leftValue = sourceBuffer[leftPtr];
                                    middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                    middleValue = sourceBuffer[middlePtr];
                                } while ((leftValue == middleValue) && (equal < DepthLimit));
                            }
                            if (leftValue <= middleValue) {
                                lastIndex = SAin[leftStart].index;
                                toLeftLcp = SAin[leftStart].lcp;
                                leftStart++;
                                toMiddleLcp = equal;
                            } else {
                                lastIndex = SAin[middleStart].index;
                                toMiddleLcp = SAin[middleStart].lcp;
                                middleStart++;
                                toLeftLcp = equal;
                            }
                        }
                        break;
                    }
                    if (equal < DepthLimit) {
                        leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                        leftValue = sourceBuffer[leftPtr];
                        middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                        middleValue = sourceBuffer[middlePtr];
                        rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                        rightValue = sourceBuffer[rightPtr];
                        equal++;
                    }
                    if (equal == DepthLimit) {
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
                        toMiddleLcp = equal;
                        toRightLcp = equal;
                        break;
                    }
                }
#ifdef COMPUTE_STATISTICS
                comparisons++;
#endif
            } else if ((leftStart < leftEnd) && (middleStart < middleEnd)) {
                int32_t leftPtr = SAin[leftStart].index;
                int32_t middlePtr = SAin[middleStart].index;
                int32_t equal = 0;
                for (; equal < DepthLimit; equal++) {
                    if (sourceBuffer[leftPtr] != sourceBuffer[middlePtr]) {
                        break;
                    }
                    leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                    middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                }
#ifdef COMPUTE_STATISTICS
                comparisons++;
#endif
                if (sourceBuffer[leftPtr] <= sourceBuffer[middlePtr]) {
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
                    toMiddleLcp = equal;
                } else {
                    lastIndex = SAin[middleStart].index;
                    toMiddleLcp = SAin[middleStart].lcp;
                    middleStart++;
                    toLeftLcp = equal;
                }
            } else if (leftStart < leftEnd) {
                lastIndex = SAin[leftStart].index;
                toLeftLcp = SAin[leftStart].lcp;
                leftStart++;
#ifdef COMPUTE_STATISTICS
                quickChecks++;
#endif
            }

            while ((leftStart < leftEnd) && (middleStart < middleEnd) && (rightStart < rightEnd)) {

                SAout[outputIndex].index = lastIndex;

                if ((toLeftLcp > toMiddleLcp) && (toLeftLcp > toRightLcp)) {
                    SAout[outputIndex++].lcp = toLeftLcp;
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else if ((toLeftLcp < toMiddleLcp) && (toMiddleLcp > toRightLcp)) {
                    SAout[outputIndex++].lcp = toMiddleLcp;
                    lastIndex = SAin[middleStart].index;
                    toMiddleLcp = SAin[middleStart].lcp;
                    middleStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else if ((toLeftLcp < toRightLcp) && (toMiddleLcp < toRightLcp)) {
                    SAout[outputIndex++].lcp = toRightLcp;
                    lastIndex = SAin[rightStart].index;
                    toRightLcp = SAin[rightStart].lcp;
                    rightStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else {
#ifdef COMPUTE_STATISTICS
                    comparisons++;
#endif
                    if (toLeftLcp < toMiddleLcp) { // && toMiddleLcp == toRightLcp
                        int32_t middlePtr = clampOnce(SAin[middleStart].index + toMiddleLcp, length);
                        int32_t rightPtr = clampOnce(SAin[rightStart].index + toMiddleLcp, length);
                        int32_t equal = toMiddleLcp;
                        for (; equal < DepthLimit; equal++) {
                            if (sourceBuffer[middlePtr] != sourceBuffer[rightPtr]) {
                                break;
                            }
                            middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                            rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                        }
                        if (sourceBuffer[middlePtr] <= sourceBuffer[rightPtr]) {
                            SAout[outputIndex++].lcp = toMiddleLcp;
                            lastIndex = SAin[middleStart].index;
                            toMiddleLcp = SAin[middleStart].lcp;
                            middleStart++;
                            toRightLcp = equal;
                        } else {
                            SAout[outputIndex++].lcp = toRightLcp;
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toMiddleLcp = equal;
                        }
                    } else if (toLeftLcp > toMiddleLcp) { // toLeftLcp == toRightLcp
                        int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                        int32_t rightPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                        int32_t equal = toLeftLcp;
                        for (; equal < DepthLimit; equal++) {
                            if (sourceBuffer[leftPtr] != sourceBuffer[rightPtr]) {
                                break;
                            }
                            leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                            rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                        }
                        if (sourceBuffer[leftPtr] <= sourceBuffer[rightPtr]) {
                            SAout[outputIndex++].lcp = toLeftLcp;
                            lastIndex = SAin[leftStart].index;
                            toLeftLcp = SAin[leftStart].lcp;
                            leftStart++;
                            toRightLcp = equal;
                        } else {
                            SAout[outputIndex++].lcp = toRightLcp;
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toLeftLcp = equal;
                        }
                    } else if (toLeftLcp > toRightLcp) { // && toLeftLcp == toMiddleLcp
                        int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                        int32_t middlePtr = clampOnce(SAin[middleStart].index + toLeftLcp, length);
                        int32_t equal = toLeftLcp;
                        for (; equal < DepthLimit; equal++) {
                            if (sourceBuffer[leftPtr] != sourceBuffer[middlePtr]) {
                                break;
                            }
                            leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                            middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                        }
                        if (sourceBuffer[leftPtr] <= sourceBuffer[middlePtr]) {
                            SAout[outputIndex++].lcp = toLeftLcp;
                            lastIndex = SAin[leftStart].index;
                            toLeftLcp = SAin[leftStart].lcp;
                            leftStart++;
                            toMiddleLcp = equal;
                        } else {
                            SAout[outputIndex++].lcp = toMiddleLcp;
                            lastIndex = SAin[middleStart].index;
                            toMiddleLcp = SAin[middleStart].lcp;
                            middleStart++;
                            toLeftLcp = equal;
                        }
                    } else { // toLeftLcp == toMiddleLcp == toRightLcp
                        int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                        int32_t middlePtr = clampOnce(SAin[middleStart].index + toLeftLcp, length);
                        int32_t rightPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                        u_int8_t leftValue = sourceBuffer[leftPtr];
                        u_int8_t middleValue = sourceBuffer[middlePtr];
                        u_int8_t rightValue = sourceBuffer[rightPtr];
                        int32_t equal = toLeftLcp;
                        while (true) {
                            if (leftValue != middleValue) {
                                if (leftValue < middleValue) {
                                    if (leftValue < rightValue) {
                                        SAout[outputIndex++].lcp = toLeftLcp;
                                        lastIndex = SAin[leftStart].index;
                                        toLeftLcp = SAin[leftStart].lcp;
                                        leftStart++;
                                        toMiddleLcp = equal;
                                        toRightLcp = equal;
                                    } else if (leftValue > rightValue) {
                                        SAout[outputIndex++].lcp = toRightLcp;
                                        lastIndex = SAin[rightStart].index;
                                        toRightLcp = SAin[rightStart].lcp;
                                        rightStart++;
                                        toLeftLcp = equal;
                                        toMiddleLcp = equal;
                                    } else {
                                        toMiddleLcp = equal;
                                        if (equal < DepthLimit) {
                                            do {
                                                equal++;
                                                leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                                leftValue = sourceBuffer[leftPtr];
                                                rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                                rightValue = sourceBuffer[rightPtr];
                                            } while ((leftValue == rightValue) && (equal < DepthLimit));
                                        }
                                        if (leftValue <= rightValue) {
                                            SAout[outputIndex++].lcp = toLeftLcp;
                                            lastIndex = SAin[leftStart].index;
                                            toLeftLcp = SAin[leftStart].lcp;
                                            leftStart++;
                                            toRightLcp = equal;
                                        } else {
                                            SAout[outputIndex++].lcp = toRightLcp;
                                            lastIndex = SAin[rightStart].index;
                                            toRightLcp = SAin[rightStart].lcp;
                                            rightStart++;
                                            toLeftLcp = equal;
                                        }
                                    }
                                } else if (middleValue > rightValue) {
                                    SAout[outputIndex++].lcp = toRightLcp;
                                    lastIndex = SAin[rightStart].index;
                                    toRightLcp = SAin[rightStart].lcp;
                                    rightStart++;
                                    toLeftLcp = equal;
                                    toMiddleLcp = equal;
                                } else if (middleValue < rightValue) {
                                    SAout[outputIndex++].lcp = toMiddleLcp;
                                    lastIndex = SAin[middleStart].index;
                                    toMiddleLcp = SAin[middleStart].lcp;
                                    middleStart++;
                                    toLeftLcp = equal;
                                    toRightLcp = equal;
                                } else {
                                    toLeftLcp = equal;
                                    if (equal < DepthLimit) {
                                        do {
                                            equal++;
                                            middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                            middleValue = sourceBuffer[middlePtr];
                                            rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                            rightValue = sourceBuffer[rightPtr];
                                        } while ((middleValue == rightValue) && (equal < DepthLimit));
                                    }
                                    if (middleValue <= rightValue) {
                                        SAout[outputIndex++].lcp = toMiddleLcp;
                                        lastIndex = SAin[middleStart].index;
                                        toMiddleLcp = SAin[middleStart].lcp;
                                        middleStart++;
                                        toRightLcp = equal;
                                    } else {
                                        SAout[outputIndex++].lcp = toRightLcp;
                                        lastIndex = SAin[rightStart].index;
                                        toRightLcp = SAin[rightStart].lcp;
                                        rightStart++;
                                        toMiddleLcp = equal;
                                    }
                                }
                                break;
                            } else if (leftValue != rightValue) {
                                if (leftValue > rightValue) {
                                    SAout[outputIndex++].lcp = toRightLcp;
                                    lastIndex = SAin[rightStart].index;
                                    toRightLcp = SAin[rightStart].lcp;
                                    rightStart++;
                                    toLeftLcp = equal;
                                    toMiddleLcp = equal;
                                } else {
                                    toRightLcp = equal;
                                    if (equal < DepthLimit) {
                                        do {
                                            equal++;
                                            leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                            leftValue = sourceBuffer[leftPtr];
                                            middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                            middleValue = sourceBuffer[middlePtr];
                                        } while ((leftValue == middleValue) && (equal < DepthLimit));
                                    }
                                    if (leftValue <= middleValue) {
                                        SAout[outputIndex++].lcp = toLeftLcp;
                                        lastIndex = SAin[leftStart].index;
                                        toLeftLcp = SAin[leftStart].lcp;
                                        leftStart++;
                                        toMiddleLcp = equal;
                                    } else {
                                        SAout[outputIndex++].lcp = toMiddleLcp;
                                        lastIndex = SAin[middleStart].index;
                                        toMiddleLcp = SAin[middleStart].lcp;
                                        middleStart++;
                                        toLeftLcp = equal;
                                    }
                                }
                                break;
                            }
                            if (equal < DepthLimit) {
                                leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                leftValue = sourceBuffer[leftPtr];
                                middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                middleValue = sourceBuffer[middlePtr];
                                rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                rightValue = sourceBuffer[rightPtr];
                                equal++;
                            }
                            if (equal == DepthLimit) {
                                SAout[outputIndex++].lcp = toLeftLcp;
                                lastIndex = SAin[leftStart].index;
                                toLeftLcp = SAin[leftStart].lcp;
                                leftStart++;
                                toMiddleLcp = equal;
                                toRightLcp = equal;
                                break;
                            }
                        }
                    }
                }
            }

            if ((leftStart == leftEnd) && (middleStart < middleEnd)) {
                leftStart = middleStart;
                leftEnd = middleEnd;
                toLeftLcp = toMiddleLcp;
            } else if ((middleStart < middleEnd) && (rightStart == rightEnd)) {
                rightStart = middleStart;
                rightEnd = middleEnd;
                toRightLcp = toMiddleLcp;
            }

            while ((leftStart < leftEnd) && (rightStart < rightEnd)) {

                SAout[outputIndex].index = lastIndex;

                if (toLeftLcp > toRightLcp) {
                    SAout[outputIndex++].lcp = toLeftLcp;
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else if (toLeftLcp < toRightLcp) {
                    SAout[outputIndex++].lcp = toRightLcp;
                    lastIndex = SAin[rightStart].index;
                    toRightLcp = SAin[rightStart].lcp;
                    rightStart++;
#ifdef COMPUTE_STATISTICS
                    quickChecks++;
#endif
                } else { // toLeftLcp == toRightLcp
                    int32_t firstPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                    int32_t secondPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                    int32_t equal = toLeftLcp;
                    for (; equal < DepthLimit; equal++) {
                        if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                            break;
                        }
                        firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                        secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                    }
#ifdef COMPUTE_STATISTICS
                    comparisons++;
#endif
                    if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                        SAout[outputIndex++].lcp = toLeftLcp;
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
                        toRightLcp = equal;
                    } else {
                        SAout[outputIndex++].lcp = toRightLcp;
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
                        toLeftLcp = equal;
                    }
                }
            }

            if (rightStart < rightEnd) {
                leftStart = rightStart;
                leftEnd = rightEnd;
                toLeftLcp = toRightLcp;
            }

            SAout[outputIndex].index = lastIndex;
            SAout[outputIndex].lcp = toLeftLcp;
            outputIndex++;

            while (leftStart < leftEnd) {
                SAout[outputIndex++] = SAin[leftStart++];
#ifdef COMPUTE_STATISTICS
                quickChecks++;
#endif
            }

#ifdef THREE_WAY_MERGE
            first += sortedLength * 3;
#else
            first += sortedLength * 2;
#endif
        } while (first < length);

        if (outputIndex != length) {
            puts("dsflj");
        }

        {
            index_t *temp = SAin;
            SAin = SAout;
            SAout = temp;
        }

#ifdef COMPUTE_STATISTICS
        printf("One pass merge time: %9ld, stride: %11u, comparisons: %11ld, quick checks: %11ld\n",
                clock() - time, sortedLength, comparisons - lastComparisons,
                quickChecks - lastQuickChecks);
#else 
        printf("One pass merge time: %ld, stride: %u\n", clock() - time, sortedLength);
#endif

#ifdef THREE_WAY_MERGE
        sortedLength *= 3;
#else
        sortedLength *= 2;
#endif

    }

#ifdef COMPUTE_STATISTICS
    printf("Comparisons:           %11ld\n", comparisons);
    printf("Quick checks:          %11ld\n", quickChecks);
#endif

    //return 0;


    for (u_int32_t i = 0; i < length; i++) {
        outputBuffer[i] = sourceBuffer[clampOnce(SAin[i].index + length - 1, length)];
    }

    if (false) {
        int64_t lcpSum = 0;
        clock_t time = clock();
        for (u_int32_t i = 1; i < length; i++) {
            lcpSum += SAin[i - 1].lcp;
            int32_t leftPtr = SAin[i - 1].index;
            int32_t middlePtr = SAin[i].index;
            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;
            }
            if (equal != SAin[i - 1].lcp) {
                puts("błąd");
            }
        }
        printf("Rotation order checking time: %ld\n", clock() - time);
        printf("Average LCP: %lf\n", (double) lcpSum / length);
    }

    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;
}

[wheezil]

I'm not sure what is the complete context of your code, but would a multi-core parallelism implementation be more helpful than low-level code optimization? Given that almost everyone has two cores now, and servers have 8-64, this would seem to be an obvious thing. Have you looked at using e.g. CILK++ or Intel TBB or OpenMP?

[Piotr Tarsa]

It's supposed to be a part of a multithreaded concurrent block-sorter. It's pointless to make a threaded algorithm which will be equally performant on four cores than best non-threaded algorithm on one core. I also care about the performance on low number of cores. I am focusing on eg performance/ watt, so both scaling and single-core performance are equally importatnt.