/****   Copyright (c) 1997-2003  Constantionos J. Siniolakis and
 ****                            Alexandros V. Gerbessiotis
 ****
 ****   Permission to use, copy, modify, and distribute this software,
 ****   and to incorporate it, in whole or in part, into other software,
 ****   is hereby granted without fee, provided that
 ****     (1) the above copyright notice and this permission notice appear in
 ****         all copies of the source code, and the above copyright notice
 ****         appear in clearly visible form on all supporting documentation
 ****         and distribution media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects; use this code on your own personal risk
 ****/

#include <stdio.h>
#include        <string.h>
#include	<memory.h>
#include 	<limits.h>
#include        "radixs.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))
#define MIN(X,Y) (((X)< (Y))?(X):(Y))

/*** FLog2 */
int FLog2 (int x)
{
  register int l;
  if (x <= 1) {
      return(0);
  } else {
      l = 0;
      while (x = Div2(x)) l++;
      return(l);
  }
}

/*** CLog2 */
int CLog2 (int x)
{
  register int l;
  if (x <= 1) {
     return(0);
  } else {
      x--; l = 1;
      while (x = Div2(x)) l++;
      return(l);
  }
}

/*** Exp2 */
int Exp2 (int x)
{
  return (int)((unsigned int)(1) << x);
}



void radixsorti(int n,unsigned int mx_int,int *sbits,int *nbits)
{
  register int tmp;
  *sbits = CLog2(mx_int); *nbits = FLog2(n);
  tmp = CDiv(*sbits, *nbits); *nbits = CDiv(*sbits, tmp);
  if (*nbits >= (*sbits - 1)) *nbits = *sbits;
}

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

  trg= (int *) malloc(n*sizeof(int));
/* check_if_null((void *)trg,"sortint_mrg","temp"); */

  if (n > mx_int) {
    r=mx_int+1;
    cnt=(int *) malloc(r*sizeof(int));
/* check_if_null((void *)cnt,"sortint_cntparameters","cnt"); */
    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((void *)cnt);
  } 
  else {
    radixsorti(n, mx_int, &sbits, &nbits);
    rbits = nbits; r = Exp2(rbits) + 1;
    cnt= (int *) malloc(r*sizeof(int));
/*  check_if_null((void *)cnt,"sortint_cntparameters","cnt"); */
    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((void *)cnt);
  }

  free((void *)trg);
}

void int32sort ( unsigned int *src, int n )
{
  radixsort(src, n, 2147483647);
}

void cintsort(unsigned int *src,int n)
{
  unsigned int *trg;
  int   cnt[256];
  register int          i,r;
  register int          j,idx;
  
  if (n <= 1) return;
  trg= (unsigned int *) malloc(n*sizeof(unsigned int));
#ifdef SANITY_CHECK
  check_if_null((void *)trg,"sortint_cntmx","trg");
#endif
  r = 256;
  memset(cnt,0,r*sizeof(int)); 
  for (j=0;j<n;j++) cnt[src[j] & 0x000000FF]++; 
  for (i=1;i<r;i++) cnt[i] += cnt[i-1];
  for (j=n-1;j>=0;j--) {
      idx= --cnt[src[j] & 0x000000FF]; 
      trg[idx]=src[j];
  }

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

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

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

  free((void *)trg);
}



void sintsort(unsigned int *src,int n)
{
  unsigned int *trg;
  static   int   cnt[65536];
  register int          i,r;
  register int          j,idx;
  
  if (n <= 1) return;
  trg= (unsigned int *) malloc(n*sizeof(unsigned int));
#ifdef SANITY_CHECK
  check_if_null((void *)trg,"sortint_cntmx","trg");
#endif
  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 = 65536;
  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);
}

void rdxinit(int n,unsigned int mx_int,int *sbits,int *nbits)
{
  register int tmp ;

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

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

  trg= (int *) malloc(n*sizeof(int));

  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 {
      rdxinit              (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 rdxint ( unsigned int *src, int n )
{
  rdx(src, n, 2147483647) ;
}

void sortint2_cntmx(unsigned int *src,int n )
{
  unsigned 	int 	*trg;
  static 	int	cnt[65536];
  register 	int	i,j,r,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 = 65536 ;
  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) ;
}



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

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

  r=65536 ;
  for(i=0;i<r;i++) cnt[i]=0;

  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 = 65536 ;
  for(i=0;i<r;i++) cnt[i]=0;
  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) ;
}



