#include        <stdio.h>
#include        <string.h>
#include	<memory.h>
#include 	<limits.h>
#include        "bsp.h"
#include        "my_bsp.h"
#define Div2(x)               ((x) >> 1)
#define Times2(x)             ((x) << 1)
#define CDiv(x, y)            (((x) != 0) ? ((((x) - 1)/(y)) + 1) : 0)
#define Mask(x, sbit, nbits)  ((x >> sbit) & ~(UINT_MAX << nbits))


void 		add(result,left,right,size)
int     	*result,*left,*right;
int     	size;

{
        register int    i;

        for (i=0;i<size;i++)
                result[i]=left[i]+right[i];
}


int      compare(left,right)
data            *left,*right;
{
        return(*left - *right);
/*  for double and strings use this
        if (*left > *right)
                return(1);
        else
                if (*left < *right)
                        return(-1);
                else
                        return(0);
*/
}

/*  for double and strings use this */
int      compare_d(left,right)
data            *left,*right;
{
/*
        if (*left > *right)
                return(1);
        else
                return(-1);
*/
        if (*left > *right)
                return(1);
        else
                if (*left < *right)
                        return(-1);
                else
                        return(0);
}

/* for integral data */
int      ncompare(data *left,data *right,int prl,int prr, int inl, int inr)
{
        return( (*left != *right)?(*left-*right):
                                  ((prl != prr)? (prl - prr): (inl -inr)) );
}

/* for double etc */
int      ncompare_d(data *left,data *right,int prl,int prr, int inl, int inr)
{
	register int	comp_result;
/*
	if (*left > *right ) comp_result = 1;
	 else
	if (*left < *right ) comp_result = -1;
	  else
			    comp_result = 0;

*/
	if (*left > *right ) comp_result = 1;
	 else
			      comp_result = -1;
        return( (*left != *right)?(comp_result):
                                  ((prl != prr)? (prl - prr): (inl -inr)) );
}



int FLog2 (int x)
{
  register int l ;

  if (x <= 1) return(0);
  else
    {
      l = 0 ;
      while (x = Div2(x)) ++l ;
      return(l);
    }
}

int CLog2 (int x)

{
  register int l ;

  if (x <= 1) return(0);
  else
    {
      --x ; l = 1 ;
      while (x = Div2(x)) ++l ;
      return(l);
    }
}

int Exp2 (int x)

{
  return (int)((unsigned int)(1) << x) ;
}


void sort_ins ( char    *src,
                int      n,
                int      sbytes,
                int    (*compare)() )

{
  char   *tmp ;
  register int     i, idx ;

  if ((n <= 1)  || (sbytes <= 0)) return ;

  tmp= (char *) malloc(sbytes*sizeof(char));
  if (tmp == NULL) {
		printf("error: malloc\n");
		exit(123);
	}


  for (i = 1 ; i < n ; ++i)
    {
      memcpy(tmp, &src[i*sbytes], sbytes) ;
      idx = i ;
      while ((idx > 0) && compare(&src[(idx-1)*sbytes], tmp) > 0)
        {
          memcpy(&src[idx*sbytes], &src[(idx-1)*sbytes], sbytes) ;
          --idx ;
        }
      memcpy(&src[idx*sbytes], tmp, sbytes) ;
    }
 
  free(tmp) ;
}



void sortint_ins ( int   *src, 
                   int    n )

{
  register int tmp ;
  register int i, idx ;

  if (n <= 1) return ;

  for (i = 1 ; i < n ; ++i)
    {
      tmp = src[i] ;
      idx = i ; 
      while ((idx > 0) && (src[idx-1] > tmp))
        {
          src[idx] = src[idx-1] ;
          --idx ; 
        }
      src[idx] = tmp ;
    }
}

void sortdbl_ins ( double  *src,
                   int      n )

{
  double tmp ;
  register int    i, idx ;

  if (n <= 1) return ;

  for (i = 1 ; i < n ; ++i)
    {
      tmp = src[i] ;
      idx = i ;
      while ((idx > 0) && (src[idx-1] > tmp))
        {
          src[idx] = src[idx-1] ;
          --idx ;
        }
      src[idx] = tmp ;
    }
}



void qs ( int  *src, int   n )
{
  register int  mid, tmp ;
  register int  l, r ;
  
  if (n <= 1) return ;

  l = 0 ; r = n - 1 ;

  mid = src[get_ir_myrandom(n)] ;

  while (l <= r)
    {
      while (src[l] < mid) ++l ;
      while (src[r] > mid) --r ;
      if (l <= r)
        {
          tmp = src[l] ;
          src[l] = src[r] ;
          src[r] = tmp ;
          ++l ; --r ;
        }
    }
  
  if (r > 0) 
    {
      if ((r + 1) < 16) sortint_ins(src, r + 1) ;
      else qs(src, r + 1) ;
    }
  if (l < n) 
    {
      if ((n - l) < 16) sortint_ins(&src[l], n - l) ;
      else qs(&src[l], n - l) ;
    }
}
           
void sortint_q ( int   *src, 
                 int   n )
{
  if (n <= 1) return ;
  init_myrandom(1000000);
 
  qs(src, n) ;
}

void qs_dbl ( double  *src, int   n )
{
  double  mid, tmp ;
  register int  l, r ;
  register int x;

  if (n <= 1) return ;

  l = 0 ; r = n - 1 ;

  mid = src[get_ir_myrandom(n)] ;

  while (l <= r)
    {
      while (src[l] < mid) ++l ;
      while (src[r] > mid) --r ;
      if (l <= r)
        {
          tmp = src[l] ;
          src[l] = src[r] ;
          src[r] = tmp ;
          ++l ; --r ;
        }
    }
 
  if (r > 0)
    {
      if ((r + 1) < 16) sortdbl_ins(src, r + 1) ;
      else qs_dbl(src, r + 1) ;
    }
  if (l < n)
    {
      if ((n - l) < 16) sortdbl_ins(&src[l], n - l) ;
      else qs_dbl(&src[l], n - l) ;
    }
}



void sortdbl_q ( double   *src,
                 int   n )
{
  if (n <= 1) return ;

  init_myrandom(1000000);
  qs_dbl(src, n) ;
}

void sort_mrg ( char    *src, 
                int   n, 
                int   sbytes, 
                int    (*compare)() )

{
  char   *temp ;
  char   *tmp, *s, *t ;
  int   phase ;
  register int idx;
  register int  i, j, l ;
  int  q, r ;


  if ((n <= 1)  || (sbytes <= 0)) return ;

  temp= (char *) malloc(n*sbytes*sizeof(char));
  if (temp == NULL) {
		printf("error: malloc\n");
		exit(123);
	}

  for (idx = 0 ; idx < n ; idx += 16    )
    {
      sort_ins(&src[idx*sbytes], MIN(16    , (n - idx)), sbytes, compare) ;
    }

  s = src ; t = temp ;
  for (phase = 16     ; phase < n ; phase = Times2(phase))
    { 
      for (idx = 0 ; idx < n ; idx += Times2(phase))
        {
          q = MIN(idx + Times2(phase), n) ; r = MIN(idx + phase, n) ;
          i = r ; j = idx ; l = idx ;  
          while ((i < q) && (j < r))
            {
              if (compare(&s[i*sbytes], &s[j*sbytes]) < 0)
                {
                  memcpy(&t[l*sbytes], &s[i*sbytes], sbytes) ;
                  ++i ; ++l ; 
                }
              else
                {
                  memcpy(&t[l*sbytes], &s[j*sbytes], sbytes) ;
                  ++j ; ++l ;
                }
             }

           if (i < q) memcpy(&t[l*sbytes], &s[i*sbytes], (q - i) * sbytes) ;    
           if (j < r) memcpy(&t[l*sbytes], &s[j*sbytes], (r - j) * sbytes) ;
        }
 
      tmp = s ; s = t ; t = tmp ;
    }

  if (s != src) memcpy(src, s, n * sbytes) ; 

  free(temp);
}

void sortint_mrg ( int  *src, 
                   int  n )
{
  int  *temp ;
  int  *tmp, *s, *t ;
  int   phase ;
  register int idx;
  register int  i, j, l ;
  int  q, r ;


  if (n <= 1) return ;

  temp= (int *) malloc(n*sizeof(int));
  if (temp == NULL) {
		printf("error: malloc\n");
		exit(123);
	}

  for (idx = 0 ; idx < n ; idx += 16    )
    {
      sortint_ins(&src[idx], MIN(16    , (n - idx))) ;
    }

  s = src ; t = temp ;
  for (phase = 16     ; phase < n ; phase = Times2(phase))
    { 
      for (idx = 0 ; idx < n ; idx += Times2(phase))
        {
          q = MIN(idx + Times2(phase), n) ; r = MIN(idx + phase, n) ;
          i = r ; j = idx ; l = idx ;  
          while ((i < q) && (j < r))
            {
              if (s[i] < s[j])
                {
                  t[l] = s[i] ;
                  ++i ; ++l ; 
                }
              else
                {
                  t[l] = s[j] ;
                  ++j ; ++l ;
                }
             }

           if (i < q) memcpy(&t[l], &s[i], (q - i) * sizeof(int)) ;    
           if (j < r) memcpy(&t[l], &s[j], (r - j) * sizeof(int)) ;
        }
 
      tmp = s ; s = t ; t = tmp ;
    }

  if (s != src) memcpy(src, s, n * sizeof(int)) ; 

  free(temp) ;  
}

void sortdbl_mrg ( double  *src, 
                   int     n )
{
  double  *temp ;
  double  *tmp, *s, *t ;
  int     phase ;
  register int idx;
  register int  i, j, l ;
  int     q, r ;


  if (n <= 1) return ;

  temp= (double *) malloc(n*sizeof(double));
  if (temp == NULL) {
		printf("error: malloc\n");
		exit(123);
	}


  for (idx = 0 ; idx < n ; idx += 16    )
    {
      sortdbl_ins(&src[idx], MIN(16    , (n - idx))) ;
    }

  s = src ; t = temp ;
  for (phase = 16 ; phase < n ; phase = Times2(phase))
    { 
      for (idx = 0 ; idx < n ; idx += Times2(phase))
        {
          q = MIN(idx + Times2(phase), n) ; r = MIN(idx + phase, n) ;
          i = r ; j = idx ; l = idx ;  
          while ((i < q) && (j < r))
            {
              if (s[i] < s[j])
                {
                  t[l] = s[i] ;
                  ++i ; ++l ; 
                }
              else
                {
                  t[l] = s[j] ;
                  ++j ; ++l ;
                }
             }

           if (i < q) memcpy(&t[l], &s[i], (q - i) * sizeof(double)) ;    
           if (j < r) memcpy(&t[l], &s[j], (r - j) * sizeof(double)) ;
        }
 
      tmp = s ; s = t ; t = tmp ;
    }

  if (s != src) memcpy(src, s, n * sizeof(double)) ; 

  free(temp) ;  
}


void sortint_cntparameters ( int           n,
                             unsigned int  mx_int,
                             int          *sbits,
                             int          *nbits )

{
  int tmp ;

  *sbits = CLog2(mx_int) ; *nbits = FLog2(n) ;
  tmp = CDiv(*sbits, *nbits) ; *nbits = CDiv(*sbits, tmp) ;
  if (*nbits >= (*sbits - 1)) *nbits = *sbits ;
}

void sortint_cnt ( unsigned int *src,
                   int             n,
                   unsigned int  mx_int )

{
  int *trg, *cnt ;
  int  lbit, sbits, nbits, rbits ;
  register int i,j, idx,r ;
 
  if ((n <= 1) || (mx_int <= 0)) return ;

  trg= (int *) malloc(n*sizeof(int));
  if (trg == NULL) {
		printf("error: malloc\n");
		exit(123);
	}


  if (n > mx_int)
    {
      r = mx_int + 1 ;
  cnt= (int *) malloc(r*sizeof(int));
  if (cnt == NULL) {
		printf("error: malloc\n");
		exit(123);
	}

      memset(cnt, 0, r * sizeof(int)) ;
      for (j = 0 ; j < n ; ++j)  cnt[src[j]]++ ;
      for (i = 1 ; i < r ; ++i) cnt[i] +=  cnt[i-1] ;
      for (j = n - 1 ; j >= 0 ; --j)
        {
          idx = --cnt[src[j]] ;
          trg[idx] = src[j] ;
        }
      memcpy(src, trg, n * sizeof(int)) ;

      free(cnt) ;
    }
  else
    {
      sortint_cntparameters(n, mx_int, &sbits, &nbits) ;
      rbits = nbits ; r = Exp2(rbits) + 1 ;
  cnt= (int *) malloc(r*sizeof(int));
  if (cnt == NULL) {
		printf("error: malloc\n");
		exit(123);
	}

      for (lbit = 0 ; lbit < sbits ; lbit += rbits)
        {
          rbits = MIN(nbits, sbits - lbit) ; r = Exp2(rbits) + 1 ;

          memset(cnt, 0, r * sizeof(int)) ;
          for (j = 0 ; j < n ; ++j)  cnt[Mask(src[j], lbit, rbits)]++ ;
          for (i = 1 ; i < r ; ++i) cnt[i] +=  cnt[i-1] ;
          for (j = n - 1 ; j >= 0 ; --j)
            {
              idx = --cnt[Mask(src[j], lbit, rbits)] ;
              trg[idx] = src[j] ;
            }

          memcpy(src, trg, n * sizeof(int)) ;
        }

      free(cnt) ;
    }

  free(trg) ;
}

/*
void sortint_cntmx ( unsigned int *src, int n )
{
  sortint_cnt(src, n, 2147483647) ;
}
*/

void sortint_cntmx ( unsigned int   *src, 
                     int           n )


{
  unsigned int *trg ;
  static int    cnt[65536] ;
  int           i, r ;
  int          j, idx ;
  
  if (n <= 1) return ;

  trg= (unsigned int *) malloc(n*sizeof(unsigned int));
  if (trg == NULL) {
		printf("error: malloc\n");
		exit(123);
	}

  r = 65536 ;
  memset(cnt, 0, r * sizeof(int)) ; 
  for (j = 0 ; j < n ; ++j)  cnt[src[j] & 0x0000FFFF]++ ; 
  for (i = 1 ; i < r ; ++i) cnt[i] += cnt[i-1] ;
  for (j = n - 1 ; j >= 0 ; --j)
    {
      idx = --cnt[src[j] & 0x0000FFFF] ; 
      trg[idx] = src[j] ;
    }

  r = 32768 ;
  memset(cnt, 0, r * sizeof(int)) ; 
  for (j = 0 ; j < n ; ++j)  cnt[trg[j] >> 16]++ ;
  for (i = 1 ; i < r ; ++i) cnt[i] += cnt[i-1] ;
  for (j = n - 1 ; j >= 0 ; --j)
    {
      idx = --cnt[trg[j] >> 16] ;
      src[idx] = trg[j] ;
    }

  free( (void *)trg) ;
}



