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

/* **************************************************************************
 *        can be compiled in 32 or 64 bits
 * ************************************************************************** */

//#define _debug
//unblock if interested in displaying of the interim data for all steps


/* **************************************************************************
 * The following are support routines for producing the test data
 * ************************************************************************** */


// returns the length of the longest Lyndon prefix of x
// function maxLyndon --------------------------------------------------------
int maxLyndon(char* x,int len) {
   int i;
   int p=1;
   int max=1;

   for(i = 1; i<len; i++) {
      if (x[i-p] != x[i]) {
         if (x[i] < x[i-p]) 
            return max;
         else{
            max = i+1; 
            p = i+1;
         }
      }
   }
   return max;
}//end maxLyndon


// compute Lyndon array
// function LyndonArr ------------------------------------------------------
void LyndonArr(char* x,int len,int la[]) {
   int i;
   for(i = 0; i < len; i++)
      la[i] = maxLyndon(&x[i],len-i);
  
   /*printf("Lyndon array:\n");
   for(i=0; i<len; i++)
     if (la[i]<10)
       printf(" %d ",la[i]);
     else
       printf("%2d ",la[i]);
   putchar('\n');*/
}//end LyndonArr


// compare two substrings x[from1..to1], x[from2..to2]
// returns < 0 if x[from1..to1] is smaller than x[from2..to2] in lexorder
// returns > 0 if x[from1..to1] is bigger than x[from2..to2] in lexorder
// returns = 0 if x[from1..to1] equals x[from2..to2]
// function ss_comp ----------------------------------------------------------
int ss_comp(char* x,int from1,int to1,int from2,int to2) {
  int t1, t2;

  for(t1=from1, t2=from2;     ; t1++, t2++) {
    if (x[t1]<x[t2]) return -1;
    if (x[t1]>x[t2]) return 1;
    if (t1==to1 || t2==to2) break;
  }
  if (t1==to1 && t2==to2) return 0;
  if (t1==to1) return -1;
  return 1;
}//end ss_comp



// compute suffix array of x[0..len-1] by brute force
// function SA -------------------------------------------------------------
void SA(char *x,int len,int sa[]) {
  int i, j, t;
  for(i=0; i<len; i++) sa[i]=i;
  for(i=len-1; i>0; i--) {
    for(j=0; j<i; j++) {
      if (ss_comp(x,sa[j],len-1,sa[j+1],len-1)>0) {
        t=sa[j];
        sa[j]=sa[j+1];
        sa[j+1]=t;
      }
    }
  }
#ifdef _debug
  printf("suffix array:\n");
  for(int i=0; i<len; i++) 
    printf("%d ",sa[i]);
  putchar('\n');
  for(int i=0; i<len; i++)
    printf("[%2d] %s\n",sa[i],&x[sa[i]]);
  putchar('\n');
  //getchar();
#endif
}//end SA


// display nicely an integer array
// function show -----------------------------------------------------------
void show(int *a,int len) {
  int i;
  for(i=0; i<len; i++) 
    if (a[i]<10)
      printf(" %d ",a[i]);
    else
      printf("%2d ",a[i]);
  putchar('\n');
}//end show



// display nicely a partially ordered Lyndon array
// function showPSLA ---------------------------------------------------------
void showPSLA(char* x,int len,int** psla,int indent=0) {
  int i;
  for(i=0; i<indent; i++) putchar(' ');
  show(psla[0],len);
  for(i=0; i<indent; i++) putchar(' ');
  show(psla[1],len);

  int k;
  for(i=0; i<len; i++) {
    if (psla[1][i]!=0) k=psla[1][i];
    printf("%2d: ",psla[0][i]);
    for(int j=0; j<len; j++) {
      if (j<psla[0][i])
        putchar('.');
      else
        putchar(x[j]);
    }
    printf("   ");
    for(int j=0; j<k; j++) printf("%c",x[psla[0][i]+j]);
    putchar('\n');
  }
}//end showPSLA


// reverse order of entries in an integer subarray 
// function reverse ---------------------------------------------------------
void reverse(int* a,int from,int to) {
  int j, k, t;
  j = (to-from+1)/2;
  for(k=0; k<j; k++) {
    t=a[from+k];
    a[from+k]=a[to-k];
    a[to-k]=t;
  }
}//end reverse


// given an intput string x[0..len-1], compute a partially ordered
// Lyndon array -- first it computes by brute force the sorted
// Lyndon array, then it reverses the order in the groups to produce
// just a partially ordered Lyndon array
// function PSLA -----------------------------------------------------------
void PSLA(char* x,int len,int **sla) {
  int i, j, t, temp;
  int* la=new int[len];

  LyndonArr(x,len,la);

  for(i=0;i<len;i++) sla[0][i]=i;

  // bubble sort of suffixes
  for(i=len-1; i>=0; i--) {
    for(j=0; j<i; j++) {
      t=ss_comp(x,sla[0][j],len-1,sla[0][j+1],len-1);
      if (t>0) {
        temp=sla[0][j];
        sla[0][j]=sla[0][j+1];
        sla[0][j+1]=temp;
      }
    }
  }
  for(i=0; i<len; i++)
    sla[1][i]=la[sla[0][i]];

  for(i=len-1; i>0;i--) {
    if (sla[1][i-1]==sla[1][i]) {
      if (ss_comp(x,sla[0][i-1],sla[0][i-1]+sla[1][i-1]-1,sla[0][i],sla[0][i]+sla[1][i]-1)==0)
        sla[1][i]=0;
    }
  }

#ifdef _debug
  printf("sla:\n");
  show(sla[0],len);
  show(sla[1],len);
  //getchar();
#endif

  // reverse order in the groups to get just a partially sorted SLA
  int start, end;
  for(i=0; i<len; i=end+1) {
    if (sla[1][i]!=0) start=i;
    for(j=i+1; j<len && sla[1][j]==0; j++);
    end=j-1;
    // group is from start..end
    if (start==end) continue;
    reverse(sla[0],start,end);
  }
}

/* **************************************************************************
 * this is the main function --- given in the input a partially ordered 
 * Lyndon array, it will compute and 'returns' via the paramater int **psla 
 * the sorted version
 * ************************************************************************** */
// function PSLA2SLA ---------------------------------------------------------
void PSLA2SLA(char* x,int len,int **psla) {
  // allocate working memory as a bunch of segments
  int *G=new int[6+9*len]; 
  int *I=G+len;
  int *SG=I+len;
  int *_G=SG+len;
  int *_SG=_G+len;
  int *_I=_SG+len;
  int *BucketStart=_I+len;
  int *BucketEnd=BucketStart+len+2;
  int *Buckets=BucketEnd+len+2;

  // in many situation we want to indicate 'empty' item, we will use len+2 for it and call it _null
  int _null=len+2;
  int len2=len+2;

  // aliases -- must be pointed to a segment before they can be used
  int *GMem;
  int *DeltaStart;
  int *DeltaEnd;
  int *Deltas;
  int *sgPtr;
  int *sgStart;
  int *inSerted;

  int i, j, g, k, delta;

  // meaning of G[i]:   I[i] belongs to group G[i]
  // meaning of SG[i]:  I[i] belongs to subgroup SG[i]
  // meaning of _G[i]:  _I[i] belongs to group _G[i]
  // meaning of _SG[i]:  _I[i] belongs to subgroup _SG[i]
  // meaning of GMem[i]:   i belongs to group GMem[i]

  // we start with (G,SG,psla[0]), so we use _G as GMem and psla[0] stands for I
  GMem=_G;
  // compute G and GMem
  for(i=g=0; i<len; i++) {
    if (psla[1][i]==0) {
      G[i]=G[i-1];
      GMem[psla[0][i]]=GMem[psla[0][i-1]];
    }else{
      G[i]=GMem[psla[0][i]]=++g;
    }
  }
  // compute SG
  for(i=0; i<len; i++) {
    if (psla[1][i]!=0) k=psla[1][i];
    delta=psla[0][i]+k;
    if (delta==len)
      SG[i]=0;
    else
      SG[i]=GMem[delta];
  }
#ifdef _debug
  printf("triples:\n");
  show(psla[0],len);
  show(G,len);
  show(SG,len);
  //getchar();
#endif

  // we might have groups 1..len, and hence subgroups 0..len, so we need arrays of length len+1
  // now we are going to radix sort triples (psla[0][i],G[i],SG[i]) -- ascending in G, ascending in SG

  // first pass, sorting (psla[0],G,SG) ascending by by SG, the result of the sort is in (_I,_G,_SG)
  for(i=0; i<len2; i++) BucketStart[i]=_null;  // initialize BucketStarts
  for(i=0; i<len; i++) {
    // put i representing (psla[0][i],G[i],SG[i]) in bucket SG[i]
    if (BucketStart[SG[i]]==_null) {
      BucketEnd[SG[i]]=BucketStart[SG[i]]=i;
      Buckets[i]=_null;
    }else{
      Buckets[BucketEnd[SG[i]]]=i;
      BucketEnd[SG[i]]=i;
      Buckets[i]=_null;
    }
  }
  // now we traverse the buckets and put the result in (_I,_G,_SG) while initializing
  // BucketStarts
  k=-1;
  for(i=0; i<len2; i++) {
    if (BucketStart[i]==_null) continue;
    k++;
    _I[k]=psla[0][BucketStart[i]];
    _G[k]=G[BucketStart[i]];
    _SG[k]=SG[BucketStart[i]];
    if (Buckets[BucketStart[i]]==_null) {
      BucketStart[i]=_null;
    }else{
      for(j=Buckets[BucketStart[i]]; j!=_null; j=Buckets[j]) {
        k++;
        _I[k]=psla[0][j];
        _G[k]=G[j];
        _SG[k]=SG[j];
      }
      BucketStart[i]=_null;
    }
  }
#ifdef _debug
  printf("triples after first pass (sorted by SG):\n");
  show(_I,len);
  show(_G,len);
  show(_SG,len);
  //getchar();
#endif

  // second pass, sorting (_I,_G,_SG) ascending by by _G, the result of the sort is in (I,G,SG)
  for(i=0; i<len; i++) {
    // put i representing (_I[i],_G[i],_SG[i]) in bucket _G[i]
    if (BucketStart[_G[i]]==_null) {
      BucketEnd[_G[i]]=BucketStart[_G[i]]=i;
      Buckets[i]=_null;
    }else{
      Buckets[BucketEnd[_G[i]]]=i;
      BucketEnd[_G[i]]=i;
      Buckets[i]=_null;
    }
  }
  // now we traverse the buckets and put the result in (I,G,SG) while initializing BucketStart
  k=-1;
  for(i=0; i<len2; i++) {
    if (BucketStart[i]==_null) continue;
    k++;
    I[k]=_I[BucketStart[i]];
    G[k]=_G[BucketStart[i]];
    SG[k]=_SG[BucketStart[i]];
    if (Buckets[BucketStart[i]]==_null) {
      BucketStart[i]=_null;
    }else{
      for(j=Buckets[BucketStart[i]]; j!=_null; j=Buckets[j]) {
        k++;
        I[k]=_I[j];
        G[k]=_G[j];
        SG[k]=_SG[j];
      }
      BucketStart[i]=_null;
    }
  }

#ifdef _debug
  printf("triples after second pass (sorted by G):\n");
  show(I,len);
  show(G,len);
  show(SG,len);
  //getchar();
#endif

  // so we have the triples properly sorted in (I,G,SG), we need to preserve psla[0] and psla[1]
  // since we do not need BucketStart,BucketEnd,Buckets any longer
  DeltaStart=BucketStart;
  DeltaEnd=BucketEnd;
  Deltas=Buckets;
  // compute inverse Delta relation
  // we do not have to initialize DeltaStart[] as it is BucketStart[] and it is already initialized
  for(i=0; i<len; i++) {
    if (psla[1][i]!=0) k=psla[1][i];
    delta=psla[0][i]+k;
    if (DeltaStart[delta]==_null) {
      DeltaStart[delta]=DeltaEnd[delta]=psla[0][i];
      Deltas[psla[0][i]]=_null;
    }else{
      Deltas[DeltaEnd[delta]]=psla[0][i];
      DeltaEnd[delta]=psla[0][i];
      Deltas[psla[0][i]]=_null;
    }
  }

#ifdef _debug
  printf("inverse delta:\n");
  for(i=0; i<len; i++) {
    if ((k=DeltaStart[I[i]])==_null) continue;
    printf("%d->%d",I[i],k);
    while(true) {
      k=Deltas[k];
      if (k==_null) break;
      printf(",%d->%d",I[i],k);
    }
    putchar('\n');
  }
  if ((k=DeltaStart[len])!=_null) {
    printf("X->%d",k);
    while(true) {
      k=Deltas[k];
      if (k==_null) break;
      printf(",X->%d",I[i],k);
    }
  }
  //getchar();
#endif

  // so we have triples sorted in (I,G,SG)
  // we have inverse delta relation in DeltaStart, Deltas
  // we need to preserve psla[1]

  // our target is psla[0]
  // we have available: _I,_G,_SG, DeltaEnd

  // we use _I for sgPtr, _G for sgStart, _SG for inSerted
  // in this loop while computing sgPtr and sgStart, we also initialize inSerted[]
  sgPtr=_I;
  sgStart=_G;
  inSerted=_SG;
  int last=0;
  sgPtr[I[0]]=0;
  sgStart[0]=0;
  inSerted[0]=0;
  for(i=1; i<len; i++) {
    inSerted[i]=0;
    if (G[i-1]==G[i]) {
      if (SG[i-1]!=SG[i]) { // starting new subgroup
        sgPtr[I[i]]=++last;
        sgStart[last]=i;
      }else{
        sgPtr[I[i]]=last;
      }
    }else{ // starting new group
      sgPtr[I[i]]=++last;
      sgStart[last]=i;
    }
  }

//-----------------------------------------------------------------------
// a few macros for copying from I[] to psla[0][]
#ifdef _debug
  #define insert(t) if (inSerted[t]==0) {psla[0][sgStart[sgPtr[t]]]=t;\
                                     sgStart[sgPtr[t]]++; inSerted[t]=1;\
                     printf("inserted %d\n",t);}
#else
  #define insert(t) if (inSerted[t]==0) {psla[0][sgStart[sgPtr[t]]]=t;\
                                     sgStart[sgPtr[t]]++; inSerted[t]=1;}
#endif
#define insert_chain(t)   if ((k=DeltaStart[t])!=_null) {\
                insert(k)\
                while(true) {\
                  k=Deltas[k];\
                  if (k==_null) break;\
                  insert(k)\
                }\
              }
//-----------------------------------------------------------------------


  // we are now copying I to psla[0] according to invers delta relat.
  insert_chain(len)
  for(i=0; i<len; i++)
    insert_chain(psla[0][i])
}//end PSLA2SLA




// function main -----------------------------------------
int main() {
  // input string
  //char x[]="abbabbabbaaabbabbaa";
  //char x[]="raindrainage";
  //char x[]="abbaabababababcdabscdbab";
  //char x[]="bscdbab";
  //char x[]="aaaa";
  //char x[]="defghabc";
  //char x[]="aabbccddeeff";
  //char x[]="abababababababababababab";
  char x[]="abbabbaaabbabbabb";

  printf("input string: %s\n",x);

  int len = strlen(x);
  
  int **psla=new int*[2];
  psla[0]=new int[len];
  psla[1]=new int[len];

  int *sa=new int[len];
  
  SA(x,len,sa);
  printf("independently computed suffix array:\n");
  show(sa,len);

  PSLA(x,len,psla); // create a test psla
  printf("input psla:\n");
  show(psla[0],len);
  show(psla[1],len);

  // sort it
  PSLA2SLA(x,len,psla);

  printf("sorted sla:\n");
  show(psla[0],len);
  show(psla[1],len);
  putchar('\n');

  // check if correct, i.e. compare the saved suffix array sa[] with psla[0] 
  int *errors=new int[len];
  bool err=false;

  for(int i=0; i<len; i++) {
    if (sa[i]==psla[0][i])
      errors[i]=1;
    else{
      err=true;
      errors[i]=0;
    }
  }
  if (err) {
    printf("error: sa[] and psla[0][] disagree:\n");
    show(sa,len);
    show(psla[0],len);
    show(errors,len);
  }else{
    printf("sa[] and psla[0][] agree\n");
  }
  //getchar();
  return 0;
}//end main