#ifdef MYDEBUG
#include  <stdio.h>
#include  <stdlib.h>
#include  <limits.h>
#include  <string.h>
#endif
#include  <math.h>
#include  "ssrt.h"
#include  "sids.h"
#include  "../avgmain.h"
#define   TWO31 ((unsigned int)1<<31)

#define CHARD char
#define STANDARD
#define INLINE_BSEARCH

double log(double);
#undef  MEMCPY
#define MEMCPY(XX,YY,SS) {DES=(XX);SRC=(YY);M=(SS); do {*DES++=*SRC++;} while (--M); }

#undef Malloc
#define Malloc(x,y)    malloc((x))
#undef Free
#define Free(x,y,z)  free((void*)(x))

#ifdef MYDEBUG
void prvector(int *ll,int dim, char *s) 
{ 
  int l,i;
fprintf(stdout,"%10s\n",s);
for(l=0;l<dim     ;l++) {
  fprintf(stdout, "[%2d]=%d ",l,ll[l]);
 }
  fprintf(stdout, "\n"); fflush(stdout);
}
void prarray(int *la,int dim, char *s) 
{
  int l,i;
fprintf(stdout,"%10s\n",s);
for(l=0;l<dim     ;l++) {
  fprintf(stdout, "[%2d] ",l);
 for(i=0;i<dim     ;i++) {
  fprintf(stdout, "%3d ,",la[i*dim     +l]);
 }
  fprintf(stdout, "\n");
  fflush(stdout);
}
 fflush(stdout);
}
#endif

int dsort
  (void  *akeys,
   int n,
   int splitters,
   char *tempres,
   int nbytes,
   void (*GENSORT)(void *, int , int , int (*)(const void*,const void*)),
   int  (*compare)(const void *,const void*),
   int  percent,
   int  alexg
  )
{
  register int oversample,skperbuck,sktotal,ksperbuck,size;
  register int i,j,k,l,s,left,right;
  register int boundary,bl,br,boff;
  register int	*bucket;   /* local bucket size */
  register int	*buck;   /* local bucket size */
  register int	*lbuck;  /* auxiliary local */

  register CHARD *sample;  /* key sample array */
  register CHARD *split, *sp;   /* splitter array */
  register CHARD *keys,  *kp;

  register int LL,RR,MM,RES;
  register char  *SRC,*DES,*PP;
  register int M=0;
  register int NBUCKETS;

  register CHARD *SRC1,*SRC2,*FIRST1,*LAST1,*FIRST2,*LAST2,*TEMP,*TMP;
  register int M1,M2;

  keys = (char *)akeys;
  NBUCKETS=splitters+1;

  /**** PART 1: Initializations */
  /* percent:percent of theor. value*/
#ifdef OLD
  oversample=(int)ceil(((percent/100.0)*(double)log(log((double)n)))+1.0);
  if (splitters <  256 )
        oversample =  oversample *4;
    else {
       if ( (splitters > 256) && (splitters <  750 ))
           oversample =  oversample * 1.5;
       else
           oversample =  MAX(oversample/4,1);
    }
#endif
  oversample = (percent/100.0)*alexg;

  skperbuck=NBUCKETS*oversample; /* sample size per bucket */
  sktotal   =skperbuck*NBUCKETS;           /* total sample size */
  ksperbuck = n/(NBUCKETS); /* bucket size in keys  */
  size=ksperbuck*nbytes;
#ifdef DEBUG_CHECK
  fprintf(stdout,"splitters=%d oversample=%d SamplePerB=%d, SampleTotal=%d n=%d KeysPerB=%d\n",
splitters,oversample,skperbuck,sktotal,n,ksperbuck); fflush(stdout);
#endif

  if ((splitters <1) || (sktotal > n/5)){
    MEMCPY(tempres,keys,n*nbytes);
    (*GENSORT)((void *)tempres,n,nbytes,compare);
    return(n);
  }
  sample = (char *) Malloc(sktotal*nbytes,"sample");       /* SampleKeys   Storage */
/*
  memset(sample,0,sktotal*nbytes);
*/
  split  = (char *) Malloc(splitters*nbytes,"split");     /* SplitterKeys Storage */
/*
  memset(split,0,splitters*nbytes);
*/
  bucket = (int *)  Malloc(NBUCKETS*sizeof(int ),"bucket");/* Bucket Size Info */
/*
  memset(bucket,0,NBUCKETS*sizeof(int));
*/

  buck  = (int *) Malloc((NBUCKETS)*(NBUCKETS)*sizeof(int ),"buck"); /* bucket size */
  lbuck = (int *) Malloc((NBUCKETS)*(NBUCKETS)*sizeof(int ),"lbuck"); /* bucket size */
/*
  memset(buck, 0,NBUCKETS*NBUCKETS*sizeof(int));
  memset(lbuck,0,NBUCKETS*NBUCKETS*sizeof(int));
*/

  /**** PART 5: Sort subsequences: All except the last have ksperbuck keys*/
  for(i=0,j=0;i<splitters;i++,j+=size){
    GENSORT((void *)(keys+j),ksperbuck,nbytes,compare);
  }
  GENSORT((void *)(keys+size*splitters),(n-ksperbuck*splitters),nbytes,compare);

/* Example: Say n= 11 S=ns/skperbuck =3
   o  o  o  o  o  o  o  o  o  o  o
            ^           ^        ^
   splitter 1           2        3 maximum
   CEIL(...) =4 FLOOR(...) =3 i.e. interspace is respectively
              3             2 
   omega_n = lgn
   oversample = r = ceil(omega_n);
   samplekeys per bucket = skperbuck = oversample * NBUCKETS
   samplekeys total      = sktotal   = skperbuck  * NBUCKETS = r NBUCKETS**2
   keys per bucket       = ksperbuck = n/NBUCKETS , last = n - ksperbuck*splitters
   Sort sktotal and start picking every skperbuck-th key for splitter.
 */
 /* Collect sample */
  for(j=0,k=0,i=0;j<n;j+=ksperbuck,k+=skperbuck,i++){
    l=((i<splitters)?ksperbuck:(n-ksperbuck*splitters)); /* # keysin i-th buck 0..splitters */
/* Determine where to change spacing from CEIL to FLOOR    */
    boundary=LEFT(l,skperbuck); 
/* CEIL  spacing for the sample on left  of boundary */
    bl      =CEIL(l,skperbuck); 
/* FLOOR spacing for the sample on right of boundary */
    br      =FLOOR(l,skperbuck);
/* j : base address -1 because 0... indexing not 1...indexing */
    boff    =j-1;   
    for(s=0;s<boundary;s++) {
      MEMCPY( (sample +(k+s)*nbytes), (keys+(boff+(s+1)*bl)*nbytes), nbytes);
    }
/* j-1 as before, boundary*bl from previous loop*/
    boff    =j-1+boundary*bl;  
/* Save one sample for the maximum */
    for(s=boundary;s<(skperbuck-1);s++) { 
     MEMCPY((sample+(k+s)*nbytes),(keys+(boff+(s+1-boundary)*br)*nbytes),nbytes);
    }
     /* Fill-in the maximum as the last sample key */
    MEMCPY( (sample+(k+skperbuck-1)*nbytes), (keys+(j+l-1)*nbytes),nbytes);
  }

 TMP = (char *) Malloc(n*nbytes,"TMP");
#ifdef SORT_THE_SAMPLE_KEYS_ALTERNATIVE
 (*GENSORT)((void *)sample,sktotal,nbytes,compare);
 /* Sort sample keys : sktotal in number */
/*  */
 /* Move here the GENSORT line above*/
  Merge
  for(i=0,j=skperbuck,s=skperbuck*nbytes;i<splitters;i++,j+=skperbuck,s+=skperbuck*nbytes) {
    MERGEFUNCTION( sample,(sample+s),j,skperbuck,nbytes,compare);
  }
 if (sktotal < (n/200)) {
  (*GENSORT)((void *)sample,sktotal,nbytes,compare);
 }
 else {
 }
*/
#endif
#ifdef SORT_THE_SAMPLE_KEYS_ALT
#endif
   for(i=1,k=2;i<NBUCKETS;i=(i<<1),k=(k<<1)) {
    for(left=0,right=i;((left<NBUCKETS)&&(right<NBUCKETS));left +=k,right+=k){
        boundary =MIN((right+i),(NBUCKETS));
        bl       = left*skperbuck*nbytes;
        br       = right*skperbuck*nbytes;
        s        = i*skperbuck;
        if (boundary == NBUCKETS) s=sktotal-right*skperbuck;
/*
        MERGEFUNCTION(
         (sample+bl),(sample+br),i*skperbuck,s,nbytes,compare);
*/
         SRC1= sample +bl;
         SRC2= sample +br;
         M1=i*skperbuck;
         M2=s;
 FIRST1=SRC1;FIRST2=SRC2;
 LAST1 =SRC1+M1*nbytes;
 LAST2 =SRC2+M2*nbytes;
 TEMP=TMP;

 for (;((FIRST1 != LAST1) && (FIRST2 != LAST2));TEMP+=nbytes)
    if (compare(FIRST2,FIRST1)<0){
      MEMCPY(TEMP,FIRST2,nbytes);
      FIRST2+=nbytes;
    }
    else {
      MEMCPY(TEMP,FIRST1,nbytes);
      FIRST1+=nbytes;
    }
 if ( FIRST1 != LAST1)
   MEMCPY(TEMP,FIRST1,(LAST1-FIRST1));
 if ( FIRST2 != LAST2)
   MEMCPY(TEMP,FIRST2,(LAST2-FIRST2));

  MEMCPY(SRC1,TMP,M1*nbytes);
  MEMCPY(SRC2,&TMP[M1*nbytes],M2*nbytes);
     }
   }


 /* Pick splitters */
  boff=(skperbuck-1)*nbytes; bl=skperbuck*nbytes;
  for(s=0;s<splitters;s++) {
/* MEMCPY( (split+s*nbytes), (sample+((s+1)*skperbuck-1)*nbytes), nbytes); */
   MEMCPY( (split+s*nbytes), (sample+boff+s*bl), nbytes); 
    bucket[s]=0;
  }
  bucket[splitters]=0;

  /* Binary search splitters into keys */
  for(k=0,kp=keys;k<=splitters;k++,kp+=size) {
    l=((k<splitters)?ksperbuck:(n-ksperbuck*splitters)); 
    sp=split+(splitters-1)*nbytes;
    for(i=splitters-1,bucket[splitters]=l;i>=0; i--,sp-=nbytes) {
#ifdef INLINE_BSEARCH
      LL=0; RR=l-1; MM=RR/2 ;PP=kp;
      while (LL<=RR) {
        RES=  compare((PP+MM*nbytes),sp);
        if (RES < 0) LL=MM+1;
          else       RR=MM-1;
        MM=(LL+RR)/2 ;
      }
      bucket[i]=LL;
#else
      bucket[i] = searchbin(sp,kp,l,compare,nbytes);
#endif
      bucket[i+1] -= bucket[i]; /* number of keys greater than or equal to splitter*/
      lbuck[(i+1)*NBUCKETS+k] = buck[(i+1)*NBUCKETS+k]=bucket[i+1];
    } 
    lbuck[k]=buck[k] = bucket[0];
  }


  /* We have sizes of buckets based on splitters from step above */
  /* Compute below first position to be occupied by a key after bucket rearrangement */
  /* Add vertically(columnwise) to determine for the l-th bucket the i-th set */
  for(i=0;i<NBUCKETS;i++) {
    lbuck[i*NBUCKETS]=0;
  }
  for(i=0;i<NBUCKETS;i++) {
    for(l=1;l<NBUCKETS;l++) {
      lbuck[i*NBUCKETS+l] = lbuck[i*NBUCKETS+l-1] + buck[i*NBUCKETS+l-1];
      bucket[i] += buck[i*NBUCKETS+l];
    }
      bucket[i] += buck[i*NBUCKETS+0];
  }
  for(i=0;i<=splitters;i++) {
    bucket[i] = lbuck[i*NBUCKETS+splitters]+buck[i*NBUCKETS+splitters];
  }
  s=bucket[0];bucket[0]=0;  /* Compute base addresses for each bucket */
  for(i=1;i<=splitters;i++) {
    l=bucket[i];
    bucket[i] = bucket[i-1]+s; 
    s=l;
  }
  /* Update lbuck by adding base address (bucket) to the  i-th column */
 for(i=0;i<NBUCKETS;i++) {
   for(l=0;l<NBUCKETS;l++) {
      lbuck[i*NBUCKETS+l] +=  bucket[i];
    }
  }

  /* Copy elements around: buck[l][i] elements to the address starting at lbuck[l][i] */
  /* Determine from indirectly. */
  for(l=0;l<=splitters;l++) {
    s=((l<splitters)?ksperbuck:(n-ksperbuck*splitters)); 
    j=s*l;
    for(i=0;i<=splitters;i++) {
     MEMCPY((tempres+lbuck[i*NBUCKETS+l]*nbytes), (keys+j*nbytes), (buck[i*NBUCKETS+l]*nbytes));
     j += buck[i*NBUCKETS+l];
    }
  }
#ifdef MYDEBUG
prarray( buck,NBUCKETS," buck/Update6");
prarray(lbuck,NBUCKETS,"lbuck/Update6");
prvector(bucket,NBUCKETS,"bucket/Update6");
#endif

#ifdef ALEXG_AJD
#endif
 for(l=0;l<NBUCKETS;l++) {
   for(i=1,k=2;i<NBUCKETS;i=(i<<1),k=(k<<1)) {
    for(left=0,right=i;((left<NBUCKETS)&&(right<NBUCKETS));left +=k,right+=k){
        /* left     is index of left  sequence to be Merged */
        /* right    is index of right sequence to be Merged */
        boff     = bucket[l]*nbytes;
        boundary =MIN((right+i),(NBUCKETS));
        bl       = lbuck[l*NBUCKETS+left];
        br       = lbuck[l*NBUCKETS+right];
        s        = lbuck[l*NBUCKETS+boundary];
        if (boundary == (NBUCKETS))
          s=((l<splitters)?(bucket[l+1]):(n));
/*
        MERGEFUNCTION(
         (tempres+bl*nbytes),(tempres+br*nbytes),(br-bl),(s-br),nbytes,compare);
*/
         SRC1= tempres+bl*nbytes;
         SRC2= tempres+br*nbytes;
         M1=br-bl;
         M2=s-br;
 FIRST1=SRC1;FIRST2=SRC2;
 LAST1 =SRC1+M1*nbytes;
 LAST2 =SRC2+M2*nbytes;
 TEMP=TMP;

 for (;((FIRST1 != LAST1) && (FIRST2 != LAST2));TEMP+=nbytes)
    if (compare(FIRST2,FIRST1)<0){
      MEMCPY(TEMP,FIRST2,nbytes);
      FIRST2+=nbytes;
    }
    else {
      MEMCPY(TEMP,FIRST1,nbytes);
      FIRST1+=nbytes;
    }
 if ( FIRST1 != LAST1)
   MEMCPY(TEMP,FIRST1,(LAST1-FIRST1));
 if ( FIRST2 != LAST2)
   MEMCPY(TEMP,FIRST2,(LAST2-FIRST2));

  MEMCPY(SRC1,TMP,M1*nbytes);
  MEMCPY(SRC2,&TMP[M1*nbytes],M2*nbytes);
     }
   }
 }
 Free((void *)TMP,n*nbytes,"TMP");
#ifdef NOT_STANDARD
 for(l=0;l<splitters;l++) {
    GENSORT((void *)(tempres+bucket[l]*nbytes),bucket[l+1]-bucket[l],nbytes,compare);
 }
    GENSORT((void *)(tempres+bucket[splitters]*nbytes),n-bucket[splitters],nbytes,compare);
#endif



/********* EXIT and FREE SPACE BELOW ***/
  Free((void *) lbuck,NBUCKETS*NBUCKETS*sizeof(int),"lbuck");
  Free((void *) buck,NBUCKETS*NBUCKETS*sizeof(int),"buck");
  Free((void *) bucket,NBUCKETS*sizeof(int),"bucket");
  Free((void *) split,splitters*nbytes,"split");
  Free((void *) sample,sktotal*nbytes,"sample");
return(n);

}




