// this is a simple C code for Algorithm 3 of the paper 
// "Reconstructing a suffix array"

// note that for the suffix array and the input string x we use C indexing 
// starting from 0 rather than the more common indexing used in the paper
// starting from 1
// Therefore we cannot use 0 in the usual meaning, and hence introduce a
// special symbol null
// We also introduce a sspecial symbol empty for use with the stack
// Note that the stack is created with size N which may not be needed,
// see the paper

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int N; // length of the input string

#define null -1
#define empty -2

void algorithm3(int* LSA,int* LCP,char* x);
void sort(int* Next,int* from,int* to,int type,char* x);
void modLCP(int* Next,int* LCP, int from,int to,int type,int last);
int* compute_LSA(char* x);


#define _debug
// uncomment the above line and recompile for debugging

//The following macors were just used for debugging
#define ShowInd() {\
	printf("     ");\
	for(int _i=0;_i<N;_i++)\
	 printf("%2d ",_i);\
	putchar('\n');fflush(stdout);}
#define ShowLSA() {\
	printf("LSA :");\
	for(int _i=0;_i<N;_i++)\
	  if (LSA[_i]==null)\
            printf(" . ");\
	  else \
           printf("%2d ",LSA[_i]);putchar('\n');fflush(stdout);}
#define ShowLCP() {\
	printf("LCP :");\
	for(int _i=0;_i<N;_i++)\
          if (LCP[_i]==null)\
	    printf(" . ");\
	  else \
            printf("%2d ",LCP[_i]);putchar('\n');fflush(stdout);}
#define ShowLCP1() {\
	printf("LCP1:");\
	for(int _i=0;_i<N;_i++)\
       	  if (LCP1[_i]==null)\
	    printf(" . "); \
	  else \
            printf("%2d ",LCP1[_i]);putchar('\n');fflush(stdout);}
#define ShowNext() {\
	printf("Next:");\
	for(int _i=0;_i<N;_i++)\
       	  if (Next[_i]==null)\
	    printf(" . ");\
	  else \
            printf("%2d ",Next[_i]);putchar('\n');fflush(stdout);}
#define ShowTail() {\
	printf("Tail:"); \
        for(int _i=0;_i<N;_i++)\
  	  if (Tail[_i]==null)\
	    printf(" . ");\
	  else \
            printf("%2d ",Tail[_i]);putchar('\n');fflush(stdout);}
#define ShowAll() {\
      	ShowInd();ShowLCP();ShowNext();ShowTail();}
#define ShowString(){\
	printf("x   :");\
       	for(int _i=0;i<N;i++)\
          printf("%c",x[_i]);fflush(stdout);}
#define TraverseNext() {\
	printf("Trav:");\
        int _i=start;\
	while(_i!=null){\
          printf("%2d ",_i);_i=Next[_i];}\
	putchar('\n');fflush(stdout);}


// create and compute the suffix array of the input string
// function compute_LSA ----------------------------------------------------
int* compute_LSA(char* x)
{  
  int* LSA=(int*) malloc(N*sizeof(int)); // allocate memory
  // the following code must be replaced by the actual computation of 
  // the suffix array of the input string rather than the following
  // suffix array for the sample input string "abcaabcdabccabdabea"
  
  int i=0;
  LSA[i++]=18; LSA[i++]=3; LSA[i++]=0; LSA[i++]=8; LSA[i++]=4;
  LSA[i++]=12; LSA[i++]=15; LSA[i++]=1; LSA[i++]=9; LSA[i++]=5;
  LSA[i++]=13; LSA[i++]=16; LSA[i++]=2; LSA[i++]=11; LSA[i++]=10;
  LSA[i++]=6; LSA[i++]=7; LSA[i++]=14; LSA[i++]=17;
  return LSA;
}// end compute_LSA


// create and compuet the LCP array of the input string
// function compute_LCP -----------------------------------------------------
int* compute_LCP(char* x)
{
  int* LCP=(int*) malloc(N*sizeof(int)); // allocate memory
  // the following code must be replaced by the actual computation of 
  // the LCP array of the input string rather than the following
  // LCP array for the sample input string "abcaabcdabccabdabea"
  int i=0;
  LCP[i++]=null; LCP[i++]=1; LCP[i++]=1; LCP[i++]=3; LCP[i++]=3; LCP[i++]=2;
  LCP[i++]=2; LCP[i++]=0; LCP[i++]=2; LCP[i++]=2; LCP[i++]=1; LCP[i++]=1;
  LCP[i++]=0; LCP[i++]=2; LCP[i++]=1; LCP[i++]=1; LCP[i++]=0; LCP[i++]=3;
  LCP[i++]=0;
  return LCP;
}//end compute_LCP


// function main ---------------------------------------------------------
int main()
{
  char x[]="abcaabcdabccabdabea";  // input string
  N=strlen(x);                     // set the global length of x
  int* LSA=compute_LSA(x);         // create & compute LSA for x  
  int* LCP=compute_LCP(x);         // create & compute LCP for x

#ifdef _debug
  printf("Input configuration:\n");
  ShowInd();
  ShowLSA();
  ShowLCP();
#endif
  algorithm3(LSA,LCP,x);          // reconstruct LSA,LCP according to the 
                                  // permutation of the alphabet
                                  // the permutation is defined in the 
				  // function sort()
#ifdef _debug
  printf("Output configuration:\n");
  ShowInd();
  ShowLSA();
  ShowLCP();
#endif

  return 0;
}//end main



// the entries in Next between 'from' and 'to' inclusive are sorted according 
// to the permutation with respect to the letter at position 'type' as all 
// suffixes in the family between 'from' and 'to' inclusive have LCP=type
// sort() is assumed to modify 'from' and 'to' accordingly (if they change 
// order)
// function sort -----------------------------------------------------------
void sort(int* Next,int* from,int* to,int type,char* x)
{
  // we are sorting according to the letter at possition type for each suffix
  // p[a]=a
  // p[b]=c
  // p[c]=d
  // p[d]=b
  // p[e]=e
  // d moves behind a or to the front

  // the following code must be replaced by a real sort

  int k, cur, prev;

  //show the family before sort, only used for debugging
#ifdef _debug
  printf("%d-family before sort:\n",type);
  for(cur=*from;cur!=null;cur=Next[cur])
    printf("%2d ",cur);
  putchar('\n');
#endif

  // traverse the family, pluck out 'd' and put it behind 'from'
  prev=*from;
  cur=Next[prev];
  while(cur!=null) {
    if (x[cur+type]=='d') {
      // cur is moving, so pluck it out
      Next[prev]=Next[cur];
      // update to if need be
      if (cur==*to)
        *to=prev;
      if (x[*from+type]=='a') {  // cur goes behind 'from', 'from' unchanged
        k=Next[*from];
        Next[*from]=cur;
        Next[cur]=k;
      }else{                    // cur becomes 'from'
        Next[cur]=*from;
        *from=cur;
      }
      break;
    }else{
      prev=cur;
      cur=Next[cur];
      continue;
    }
  }//endwhile
  // so the fmaily is now sorted properly, 'from' and 'to' are updated properly

  // show the family after sort, only used for debugging
#ifdef _debug
  printf("%d-family after sort:\n",type);
  for(cur=*from;cur!=null;cur=Next[cur])
    printf("%2d ",cur);
  putchar('\n');
#endif
}//end sort



// the LCP values in the freshly sorted family between 'from' and 'to' 
// inclusive are updated properly, 'last' is the LCP value of the very last 
// member of the family before the sort of the family
// function modLCP ----------------------------------------------------------
void modLCP(int* Next,int* LCP, int from,int to,int type,int last) 
{
  int i;

  // show LCP of the family before modification, only used for debugging
#ifdef _debug
  printf("LCP before:\n");
  for(i=from;i!=null;i=Next[i])
    printf("%2d ",LCP[i]);
  putchar('\n');
#endif
  
  for(i=from; i!=null; i=Next[i])
    if (LCP[i] < type)
      LCP[i]=type;
  LCP[to]=last;

  // show LCP of the family after modification, only used for debugging
#ifdef _debug
  printf("LCP after:\n");
  for(i=from;i!=null;i=Next[i])
    printf("%2d ",LCP[i]);
  putchar('\n');
#endif
}//end modLCP


// these are Stack handling macros
#define Pop()  ((top<0)?(empty):(Stack[top--]))
#define Top()  ((top<0)?(empty):(Stack[top]))   // show top of the stack
#define Top_1() ((top<=0)?(empty):(Stack[top-1]))   // show one before top 
#define Push(a)  Stack[++top]=a
#define ShowStack()\
        for(int _i=top;_i>=0;_i--)\
          if (_i==top)\
            printf("Stack:\n%d\n",Stack[_i]);\
         else\
            printf("%d\n",Stack[_i])
    
#define lcp(cur) ((Tail[cur]==null)?(LCP[cur]):(LCP[Tail[Tail[cur]]]))
#define nextTop() ((Top()==null)?(start):(Next[Top()]))

// input LSA,LCP of the input string x
// output modified LSA,LCP according to the permutation as defined in sort()
// function algorithm3 -----------------------------------------------------
void algorithm3(int* LSA,int* LCP1,char* x)
{
  int* Next=new int[N];   // create Next[]
  int* LCP=new int[N];    // create LCP[]
  int* Stack=new int[N];  // create Stack - the size of stack need not bee N (see the paper), 
                          // however for simplicity here we are suing N
  int top=null;           // initialize Stack to empty
  int* Tail;              // will become Tail[]

  int start;              // "pointer" to beginning of the list Next[]
  int cur;                // "pointer" to current "node" of the list Next[]
  int prev;               // "pointer" to previous "node" of the list Next[]
  int from;               // "pointer" to beginning of the family to be sorted
  int LE;                 // "pointer" to node immediately preceeding 'from'
  int RE;                 // "pointer" to node immediately following the end 
                          // of the family
  int type;               // the size of the common prefix of all suffixes in 
                          // the family
  int modify_start;       // a flag to indicate if 'start' needs to be modified
                          // after sorting the family
  int small_family;       // a flag indicating if the family has two (small) 
                          // or more members
  int last;               // the LCP of the very last member of the family 
                          // before sort
  int famend;             // "pointer" to the last member of the family 
                          // after sort

  int i, j, type1;        // auxiliary variables
  
  
  // transform LSA to Next[], and modify LCP1 into LCP
  start=LSA[0];
  LCP[start]=LCP1[1];
  for(i=0;i<N-1;i++) {
    Next[LSA[i]]=LSA[i+1];
    LCP[LSA[i]]=LCP1[i+1];
  }
  Next[LSA[N-1]]=null;
  LCP[LSA[N-1]]=null;

  // use the memory of LCP1 for Tail and set it up
  Tail=LCP1;
  for(i=0;i<N;i++)
    Tail[i]=null;

  // initialize prev,cur and Stack for the first family and compute its type
  prev=null;
  cur=start;
  while(lcp(cur)==0) {
    prev=cur;
    cur=Next[cur];
  }
  Push(prev);
  type=lcp(cur);

  while(1) {

  // IDENTIFY AND EXTRACT FAMILY
    if (type==lcp(cur)) {
      prev=cur;
      cur=Next[cur];
      continue;
    }else if (type<lcp(cur)) {  // lcp went up, so we start a new family
      type=lcp(cur);
      Push(prev);
      prev=cur;
      cur=Next[cur];
      continue;
    }else{         // lcp went down, so we have identified a family
        // the family starts at nextTop() and ends at cur
      from=nextTop();
      // the left end is Top(), the right end is Next[cur]
      LE=Top();
      RE=Next[cur];
#ifdef _debug
      printf("%d-family from %d to %d (LE=%d,RE=%d)\n",type,from,cur,LE,RE);
#endif
      // extract it 
      Next[cur]=null;
    }

#ifdef _debug
    printf("Before sort:\n");
    ShowAll();
#endif
    
  // SORT THE FAMILY
    // remember if the family begins with 'start' as it may move and will 
    // need to be updated
    modify_start=(from==start);
    // remember if we are sorting small family or not, will be needed for 
    // update of 'prev' and 'cur'
    small_family=(Next[from]==cur);
    // rember last LCP
    last=LCP[cur];

    // sort the family
    sort(Next,&from,&cur,type,x);
    // modify the LCP of the sorted family
    modLCP(Next,LCP,from,cur,type,last);

    // 'from' and 'cur' were properly updated by sort()
    // update 'start', if need be
    if (modify_start) start=from;
    // update 'prev'
    if (small_family)
      prev=from;
    else
      for(prev=from;Next[prev]!=cur;prev=Next[prev]);

    // update the attachment of the family to LE
    if (LE!=null)
      Next[LE]=from;

#ifdef _debug
    printf("Before flattening:\n");
    ShowAll();
#endif

  // FLATTEN THE FAMILY 
  // and compute the family end
    famend=cur;
    for(i=from; i!=null; i=Next[i]) {
      if (Tail[i]!=null) {
        if (i==cur)
          famend=Tail[Tail[i]];
        Next[Tail[Tail[i]]]=Next[i];
        Next[i]=Tail[i];
        Tail[Tail[i]]=null;
        Tail[i]=null;
      }
    }

#ifdef _debug
    printf("before verticalization:\n");
    ShowAll();
#endif

  // VERTICALIZE THE FAMILY
    Tail[from]=Next[from];
    Tail[Tail[from]]=famend;
    Next[from]=RE;
    cur=from;
    prev=LE;

#ifdef _debug
    printf("After verticalization:\n");
    ShowAll();
#endif

    // multipop()
    // is implemented her differently than as discribed in the paper
    // rather than a function, it is simply handled through flow control
    // 'break' corresponds to multipop() returning 'false', 'continue'
    // corresponds to multipop() returning 'true'
    // The cases A1-A4, B1-B3 refer to the cases as described in the paper
 
  if (Top_1()==empty) {
    if (LCP[famend]==0) { // cases A1, A2, A3
      Pop();
      prev=cur;
      cur=Next[cur];
      if (cur==null) break; // break the while loop to perform the final 
                            // sort of the 0-family
                            // 'prev' points to the very last item in Next
      type=LCP[cur];
      while(type==0) {
        prev=cur;
        cur=Next[cur];
        type=LCP[cur];
      }
      if (type==null) {   // we had one or more stand-alone suffixes at 
	                  // the end of the arry
      prev=cur;           // since 'cur' points to the very last iten in Next, 
                          // we make 'prev' point there
      break;              // break the while loop to perform the final sort 
                          // of the 0-family    
      }
      // so we have not reached the end, we are starting a new family
      Push(prev);
      continue;
    }else{ // case A4
      type=LCP[famend];    // decrease type
      prev=cur;
      cur=Next[cur];
      continue;
    }
  }else{  // hence Top_1()!=empty
    if (Top_1()==null)
      type1=LCP[start];
    else
      type1=LCP[Next[Top_1()]];

    if (LCP[famend]>type1) { // case B1 
      type=LCP[famend];    // decrease type
      prev=cur;
      cur=Next[cur];
      continue;
    }else if (LCP[famend]==type1) { // case B2
      Pop();
      type=type1;
      prev=cur;
      cur=Next[cur];
      if (cur==null) break; // break the while loop to perform the final 
                            // sort of the 0-family  
                            // 'prev' points to the very last item in Next
      continue;
    }else{ // hence LCP[famend]<type1  --- case B3
      Pop();
      type=type1;
      continue;
    }
  }
  }//end while


  // final sort of the 0-family
#ifdef _debug
  printf("0-family from %d to %d (LE=%d,RE=%d)\n",start,prev,null,null);
  printf("Before sort:\n");
  ShowAll();
#endif

  sort(Next,&start,&prev,0,x); 
  
#ifdef _debug
  printf("before flattening:\n");
  ShowAll();
#endif

  for(i=start; i!=null; i=Next[i]) {
    if (Tail[i]!=null) {
        Next[Tail[Tail[i]]]=Next[i];
        Next[i]=Tail[i];
        Tail[Tail[i]]=null;
        Tail[i]=null;
    }
  }

#ifdef _debug
  printf("done:\n");
  ShowAll();
#endif

  // transform Next to LSA
  for(j=0,i=start; i!=null; j++,i=Next[i])
    LSA[j]=i;
  // transform LCP to LCP1
  LCP1[0]=null;
  for(j=1,i=start; Next[i]!=null; j++,i=Next[i])
    LCP1[j]=LCP[i];
  
}//end algorithm3