//-------------------------------------------------------------
// Larsson-Sadakane's algorithm
// To compile, perform: gcc -std=c99 SALS.s -o SALS
// To run, SALS　<length of the string>. e.g., SALS 100
//-------------------------------------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <string.h>
#include "MT.h"     // Use the Mersenne Twister.

#define DEBUG
// これを定義すると #ifdef <code> #endif の中の <code> が実行されるので、デバッグの時だけ動かしたいコードを入れるのに便利です. デバッグしないときはコメントではずしましょう.

void dump_SA_ISA(int *SA, int *ISA, int h, int n){
    printf("h = %i, n = %i\n", h, n);
    printf("SA = ");  for(int i=0; i<n; i++) printf("%i ", SA[i]);  printf("\n");
    printf("ISA = "); for(int i=0; i<n; i++) printf("%i ", ISA[i]); printf("\n");
}

void radixSort_LS(int *S, int *SA, int *ISA, int n){
    int m = 0;
    for(int i = 0; i < n; i++){ //
        if(m < S[i]) m = S[i];
    }
    m++;    // m=5 for {$,A,C,G,T}
    int *count = (int *)malloc(sizeof(int) * m);
    for(int i=0; i<m; i++)      count[i]=0;              //
    for(int i=0; i<n; i++)      count[S[i]]++;           //
    for(int i=1; i<m; i++)      count[i] += count[i-1];  //
    for(int i=n-1; 0<=i; i--)   ISA[i] = count[S[i]]-1;  //
    for(int i=n-1; 0<=i; i--)   SA[--count[S[i]]] = i;   //
    free(count);
}

void ternary_split_quicksort_LS(int *SA, int l, int r, int *ISA, int h){
    if(l > r)   return;
    if(l == r){
        ISA[SA[l]] = l; return;
    }
    int pivot = l + genrand_int32() % (r-l+1);
    int v = ISA[SA[pivot]+h];
    int i=l; int mi=l; int j=r; int mj=r; int tmp;
    for(;;){
        for( ; i <= j && ISA[SA[i]+h] <= v; i++)
            if(ISA[SA[i]+h] == v){
                tmp = SA[i]; SA[i] = SA[mi]; SA[mi] = tmp; mi++;
            }
        for( ; i <= j && v <= ISA[SA[j]+h]; j--)
            if(ISA[SA[j]+h] == v){
                tmp = SA[j]; SA[j] = SA[mj]; SA[mj] = tmp; mj--;
            }
        if(i > j) break;
        tmp = SA[i]; SA[i] = SA[j]; SA[j] = tmp; i++; j--;
    }
    for(mi--, i--; l <= mi; mi--, i--){
        tmp = SA[i]; SA[i] = SA[mi]; SA[mi] = tmp;
    }
    for(mj++, j++; mj <= r; mj++, j++){
        tmp = SA[j]; SA[j] = SA[mj]; SA[mj] = tmp;
    }
    ternary_split_quicksort_LS(SA, l, i, ISA, h); //
    for(int k=i+1; k<j; k++)  ISA[SA[k]] = j-1;   //
    ternary_split_quicksort_LS(SA, j, r, ISA, h); //
}

void SALS(int *S, int *SA, int n){
    int *ISA = (int *)malloc(sizeof(int) * n);
    // Basic step
    radixSort_LS(S, SA, ISA, n);
    // Iterative step
    int iteration;
    int h=1;
    do{
        #ifdef DEBUG
        dump_SA_ISA(SA, ISA, h, n);
        #endif
        iteration = 0;
        int next_i;
        for(int i=0; i<n; i=next_i){
            next_i = ISA[SA[i]]+1;  //
            if(i < ISA[SA[i]]){     //
                ternary_split_quicksort_LS(SA, i, ISA[SA[i]], ISA, h);  //
                iteration = 1;      //
            }
        }
        h = h*2;    //
    }while(iteration == 1);
    free(ISA);
}

int less(int *S, int aI, int aJ, int n){
    int i, j;
    for(i=aI, j=aJ; i<n && j<n; i++, j++){
        if(S[i] < S[j]) return(1); else if(S[i] > S[j]) return(0);
    }
    //
    if(aI < aJ) return(0); else return(1);
}


int main(int argc, char *argv[])
{
    if(argc != 2){
        fprintf(stderr, "Usage: SALS <Length of an input string>\n");
        exit(EXIT_FAILURE);
    }
    int n = atoi(argv[1]);
    if(n <= 0){
        fprintf(stderr, "The argument must be a positive integer\n");
        exit(EXIT_FAILURE);
    }
    
    int *S  = (int *)malloc(sizeof(int) * n);  // Use the heap for the input string.
    if(S == NULL)  exit(EXIT_FAILURE);
    int *SA = (int *)malloc(sizeof(int)*n);
    if(SA == NULL){  free(S); exit(EXIT_FAILURE); }
    
    // We encode $,A,C,G,T by 0,1,2,3,4 so that  $=0 < A=1 < C=2 < G=3 < T=4
    for(int i=0; i<n-1; i++)
        #ifdef DEBUG
        S[i] = 1;   //
        #else
        S[i] = genrand_int32() % 4;     //
        #endif
    S[n-1] = 0;     //
    #ifdef DEBUG
    printf("S = "); for(int i=0; i<n; i++) printf("%i ", S[i]); printf("\n");
    #endif
    
    struct timeval s, e;    //
    gettimeofday(&s, NULL); //
    SALS(S, SA, n);         //
    gettimeofday(&e, NULL); //
    float comp_time = (e.tv_sec - s.tv_sec) + (e.tv_usec - s.tv_usec)*1.0E-6;
    
    for(int i=0; i<n-1; i++){ //
        if( less(S, SA[i], SA[i+1], n) == 0 ){   //
            fprintf(stderr, "Failure to sort the target array.\n");
            free(S); free(SA); exit(EXIT_FAILURE);
        }
    }
    printf("The input of %i in size was sorted successfully in %f sec.\n", n, comp_time);
    
    free(S);
    free(SA);
    return EXIT_SUCCESS;
}