#include <iostream>
#include <string>
#include <cmath>
#include <ctime>

//
// Larsson-Sadakane Algorithm
//
void radixSort_LS(int* S, int Slen, int* SA, int* ISA){	
		// Treat each element as one digit, execute radix sort, and update ISA.
		int n = Slen;
		int m=0;
		for(int i=0; i<n; i++) 	// Calculate the maximum.
			if(m < S[i]) m = S[i];
		m++; // Total number of elements
		int* count = new 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; }  // update ISA
		for(int i = n-1; 0 <= i; i--){ SA[count[S[i]]-1] = i; count[S[i]]--; }
	}	

void ternary_split_quicksort_LS
					(int* SA, int l, int r, int* ISA, int h) {
		if(l > r) return;	 // Exit when the empty range is input.
		if(l == r){ // Set ISA[SA[l]] to l when the input is singleton.
			ISA[SA[l]] = l; return; }
		// Choose a pivot from ISA[SA[]+h].		
		int v = ISA[SA[l + ((rand()*RAND_MAX+rand())%(r-l+1)) ]+h];  
			// rand() generates a 2-byte integer at random. RAND_MAX = 32767 = 2^15 - 1.
		int i = l; int mi = l; int j = r; int mj = r; int tmp;
		for (;;) {  // Compare values according to ISA[SA[]+h]
			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;  // update ISA
		ternary_split_quicksort_LS(SA, j, r, ISA, h);
	}

int* compute_SA_LS(int* S, int Slen){
		// Assume that all elements in S are positive except that the last is zero.
		int n = Slen;
		int* SA  = new int[n];
		int* ISA = new int[n];
		// Basic step.
		radixSort_LS(S, Slen, SA, ISA);
		//Iterative step.
		bool iteration; int h=1;
		do{		
			iteration = false;
			int next_i;
			for(int i=0; i < n; i = next_i){
				next_i = ISA[SA[i]]+1;		// SA と ISA が下記の ternary_split で更新されるので現在の値を取っておく
				if(i < ISA[SA[i]]){			// If some blocks are divided.
					ternary_split_quicksort_LS(SA, i, ISA[SA[i]], ISA, h);
					// Repeat the while-loop. 
					iteration = true;
				}
			}
			h = h*2;
		}while(iteration);
		return(SA);
	}

void generateRandomNumbers(int* target, int targetLen){
		// Generate an array in which elements are selected from {1,2,3,4} at random.
		for(int i=0; i<targetLen-1; i++) 	target[i] = rand()%4 + 1;
		// Set the last element to 0 !
		target[targetLen-1] = 0;
}


void main(int argc, char** argv){

		int slen = 100000;
		int* s = new int[slen]; 	
		generateRandomNumbers(s, slen);		// Assume that the last element is 0 !

		clock_t start_time, end_time, etime;
		start_time = clock();
		int* SA_LS = compute_SA_LS(s, slen);
		end_time = clock();
		etime = (end_time - start_time);  // CLOCKS_PER_SEC;
		cout << "Time for Larsson-Sadakane\t" << etime << " msec" << endl;

}