#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <stdio.h>
#include <queue>
#include <string>

using namespace std;

typedef unsigned char byte;

byte *comparefunc_buf = NULL;
int comparefunc_bufsize = 0;

void bwt_transform(byte *input, byte *output, int size, int *orig_index);
void bwt_reverse(byte *input, byte *output, int size, int orig_index);

#define BLOCK_SIZE 200000

int main(int argc, char **argv) {
	
	if (argc !=2 ) {
		printf ("Usage: %s inputfilename\n", argv[0]);
		exit(0);
	}
	FILE * in_file;
	
	in_file = fopen(argv[1], "r");
	if (in_file == 	NULL) {
		printf("Cannot input input file.");
		exit(0);
	}
	
	byte orig_seq[BLOCK_SIZE];
	
	while (1) {
		int size = fread( (char *) orig_seq, 1, BLOCK_SIZE, in_file );
        if ( size == 0 ) break; 
        
		printf("BWT processing chunk of %d bytes\n", size );
        fwrite( (char *) &orig_seq, 1, size, stdout );
		
		int orig_index;
		byte bwt_transf_seq[size];
		byte test[size];
    
		bwt_transform (orig_seq, bwt_transf_seq, size, &orig_index);
		bwt_reverse (bwt_transf_seq, test, size, orig_index);

		assert (!memcmp (orig_seq, test, size));
		printf ("Result is the same as input, that is: <%.*s>\n", size, test);
		// Print out encode/decode results:
		printf ("Input : <%.*s>\n", size, orig_seq);
		printf ("Output: <%.*s>\n", size, bwt_transf_seq);

	}
	
	
    return 0;
}


int comparefunc(const void *l, const void *r)
{
    int li = *(const int*)l ;
	int ri = *(const int*)r ;
	int ac = comparefunc_bufsize;
	
    while (comparefunc_buf[li] == comparefunc_buf[ri])
    {
        if (++li == comparefunc_bufsize)
            li = 0;
        if (++ri == comparefunc_bufsize)
            ri = 0;
        if (!--ac)
            return 0;
    }
    if (comparefunc_buf[li] > comparefunc_buf[ri])
        return 1;
    else
        return -1;
}

void bwt_transform(byte *input, byte *output, int size, int *orig_index)
{
    int indices[size];
    int i;

    for(i=0; i<size; i++)
        indices[i] = i;
    comparefunc_buf = input;
    comparefunc_bufsize = size;
    qsort (indices, size, sizeof(int), comparefunc );

    for (i=0; i<size; i++)
        output[i] = input[(indices[i]+size-1)%size];
		
    for (i=0; i<size; i++)
    {
        if (indices[i] == 1) {
            *orig_index = i;
            return;
        }
    }
}

void bwt_reverse(byte *input, byte *output, int size, int orig_index)
{
    byte F[size];
    int buckets[256];
    int i,j,k;
    int indices[size];

    for (i=0; i<256; i++)
        buckets[i] = 0;
		
    for (i=0; i<size; i++)
        buckets[input[i]] ++;
		
    for (i=0,k=0; i<256; i++)
        for (j=0; j<buckets[i]; j++)
            F[k++] = i;

    for (i=0,j=0; i<256; i++)
    {
        while (i>F[j] && j<size)
            j++;
        buckets[i] = j; 
    }
	
    for(i=0; i<size; i++)
        indices[buckets[input[i]]++] = i;
		
    for(i=0,j=orig_index; i<size; i++)
    {
        output[i] = input[j];
        j=indices[j];
    }
}