#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "nbtri.h"

/* 
++ 
  Module	Tri de tableaux (C)
  Lib	LibsUtil
  include	nbtri.h

	Tri de tableaux et indexation.
--
*/

/*=========================================================================*/
/*
++
  void HeapSort(int n,double *ra_int)
	On reordonne par ordre numerique croissant
	le tableau ra_int[n]: Numerical Recipes mode C.
--
*/
void HeapSort(int n,double *ra_int)
{
   int l,j,ir,i;
   double rra,*ra;

   /* attention, Numerical recipes prend des tableaux de 1 a n on remet
      de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
   ra = ra_int-1;

   l=(n >> 1)+1;
   ir=n;
   for (;;) {
      if (l > 1)
         rra=ra[--l];
      else {
         rra=ra[ir];
         ra[ir]=ra[1];
         if (--ir == 1) {
            ra[1]=rra;
            return;
         }
      }
      i=l;
      j=l << 1;
      while (j <= ir) {
         if (j < ir && ra[j] < ra[j+1]) ++j;
         if (rra < ra[j]) {
            ra[i]=ra[j];
            j += (i=j);
         }
         else j=ir+1;
      }
      ra[i]=rra;
   }
}

/*=========================================================================*/
/*
++
  void HeapSortF(int n,float *ra_int)
	On reordonne par ordre numerique croissant
	le tableau ra_int[n]: Numerical Recipes mode C.
--
*/
void HeapSortF(int n,float *ra_int)
{
   int l,j,ir,i;
   float rra,*ra;

   /* attention, Numerical reciepes prend des tableaux de 1 a n on remet
      de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
   ra = ra_int-1;

   l=(n >> 1)+1;
   ir=n;
   for (;;) {
      if (l > 1)
         rra=ra[--l];
      else {
         rra=ra[ir];
         ra[ir]=ra[1];
         if (--ir == 1) {
            ra[1]=rra;
            return;
         }
      }
      i=l;
      j=l << 1;
      while (j <= ir) {
         if (j < ir && ra[j] < ra[j+1]) ++j;
         if (rra < ra[j]) {
            ra[i]=ra[j];
            j += (i=j);
         }
         else j=ir+1;
      }
      ra[i]=rra;
   }
}

/*=========================================================================*/
/*
++
  void HeapSortF2(int n,float *ra_int,float *ra2_int)
	On reordonne par ordre numerique croissant
	le tableau ra_int[n] et ra2_int en parralelle.
	Numerical Recipes mode C.
--
*/
void HeapSortF2(int n,float *ra_int,float *ra2_int)
{
   int l,j,ir,i;
   float rra,rra2,*ra,*ra2;

   /* attention, Numerical reciepes prend des tableaux de 1 a n on remet
      de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
   ra = ra_int-1;
   ra2 = ra2_int-1;

   l=(n >> 1)+1;
   ir=n;
   for (;;) {
      if (l > 1)
        {rra=ra[--l];   rra2=ra2[l];}
      else {
         rra=ra[ir];    rra2=ra2[ir];
         ra[ir]=ra[1];   ra2[ir]=ra2[1];
         if (--ir == 1) {
            ra[1]=rra;  ra2[1]=rra2;
            return;
         }
      }
      i=l;
      j=l << 1;
      while (j <= ir) {
         if (j < ir && ra[j] < ra[j+1]) ++j;
         if (rra < ra[j]) {
            ra[i]=ra[j];  ra2[i]=ra2[j];
            j += (i=j);
         }
         else j=ir+1;
      }
      ra[i]=rra;   ra2[i]=rra2;
   }
}

/*=========================================================================*/
/*
++
  int_4 tri_double ( double *tab, int_4 *indx,int_4 N)
	Methode de tri sans finesse (double boucles).
| entree : tab -> tableau de double de longueur N (0 a N-1 )
| sortie : indx -> tableau d'entiers de longueur N : la case i   
|          contient le classement du ieme element dans tab  
| Retourne:  0 si tri impossible
|            1 si tri reussi
--
*/
int_4 tri_double ( double *tab, int_4 *indx,int_4 N)
{
int_4 i,j,k;

if (N<=0) return (-1);

for (i=0; i<N ; i++) indx[i]=i;
if (N==1) return (N);

/* classement dans l'ordre croissant */
for (i=1; i< N ; i++)
  { 

  if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
    { 
     k = indx[i-1];
     indx[i-1] = indx[i];
     indx[i] = k;

     if ( i > 1 )
        for (j=i ; j>=1 ; j--)
          { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
            { 
            k = indx[j-1];
            indx[j-1] = indx[j];
            indx[j] = k;
    } } } }

return (N) ;
}

/*=========================================================================*/
/*
++
  int_4 tri_float ( float *tab, int_4 *indx,int_4 N)
	Methode de tri sans finesse (double boucles).
| entree : tab -> tableau de flottant de longueur N (0 a N-1 )
| sortie : indx -> tableau d'entiers de longueur N : la case i   
|          contient le classement du ieme element dans tab  
| Retourne:  0 si tri impossible
|            1 si tri reussi
--
*/
int_4 tri_float ( float *tab, int_4 *indx,int_4 N)
{
int_4 i,j,k;

if (N<=0) return (-1);

for (i=0; i<N ; i++) indx[i]=i;
if (N==1) return (N);

/* classement dans l'ordre croissant */
for (i=1; i< N ; i++)
  { 

  if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
    { 
     k = indx[i-1];
     indx[i-1] = indx[i];
     indx[i] = k;

     if ( i > 1 )
        for (j=i ; j>=1 ; j--)
          { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
            { 
            k = indx[j-1];
            indx[j-1] = indx[j];
            indx[j] = k;
    } } } }

return (N) ;
}

/*=========================================================================*/
/*
++
  int_4 tri_entier ( int_4 *tab,int_4 *indx,int_4 N)
	Methode de tri sans finesse (double boucles).
| entree : tab -> tableau d'entiers de longueur N (0 a N-1 )
| sortie : indx -> tableau d'entiers de longueur N : la case i   
|          contient le classement du ieme element dans tab  
| Retourne:  0 si tri impossible
|            1 si tri reussi
--
*/
int_4 tri_entier ( int_4 *tab,int_4 *indx,int_4 N)
{
int_4 i,j,k;

if (N<=0) return (-1);

for (i=0; i<N ; i++) indx[i]=i;
if (N==1) return (N);

/* classement dans l'ordre croissant */
for (i=1; i< N ; i++)
  { 

  if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
    { 
     k = indx[i-1];
     indx[i-1] = indx[i];
     indx[i] = k;

     if ( i > 1 )
        for (j=i ; j>=1 ; j--)
          { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
            { 
            k = indx[j-1];
            indx[j-1] = indx[j];
            indx[j] = k;
    } } } }

return (N) ;
}

/*===================================================================================*/
/*
++
  int qSort_Float(const void *a1,const void *a2)
	Fonction de tri de `float' a utiliser dans qsort.
--
*/
int qSort_Float(const void *a1,const void *a2)
{
if( *((float *) a1) < *((float *) a2) ) return(-1);
if( *((float *) a1) > *((float *) a2) ) return( 1);
return(0);
}

/*===================================================================================*/
/*
++
  int qSort_Dble(const void *a1,const void *a2)
	Fonction de tri de `double' a utiliser dans qsort.
--
*/
int qSort_Dble(const void *a1,const void *a2)
{
if( *((double *) a1) < *((double *) a2) ) return(-1);
if( *((double *) a1) > *((double *) a2) ) return( 1);
return(0);
}

/*===================================================================================*/
/*
++
  int qSort_Int(const void  *a1,const void  *a2)
	Fonction de tri de `int' a utiliser dans qsort.
--
*/
int qSort_Int(const void  *a1,const void  *a2)
{
if( *((int *) a1) < *((int *) a2) ) return(-1);
if( *((int *) a1) > *((int *) a2) ) return( 1);
return(0);
}

/*===================================================================================*/
/*
++
  int qSort_Ushort(const void  *a1,const void  *a2)
	Fonction de tri de `unsigned short' a utiliser dans qsort.
--
*/
int qSort_Ushort(const void  *a1,const void  *a2)
{
if( *((unsigned short *) a1) < *((unsigned short *) a2) ) return(-1);
if( *((unsigned short *) a1) > *((unsigned short *) a2) ) return( 1);
return(0);
}

/*===================================================================================*/
/*
++
  int qSort_Short(const void  *a1,const void  *a2)
	Fonction de tri de `short' a utiliser dans qsort.
--
*/
int qSort_Short(const void  *a1,const void  *a2)
{
if( *((short *) a1) < *((short *) a2) ) return(-1);
if( *((short *) a1) > *((short *) a2) ) return( 1);
return(0);
}

/*===================================================================================*/
/*
++
  void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
	Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
	such that arr[indx[j]] is in ascending order for j=1,2,,N.
	The input quantities n and arr are not changed.
--
*/
#define SWAP_INDEXR4(a,b) itemp=(a);(a)=(b);(b)=itemp;
#define M 7
#define NSTACK 50
void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
/* encore du Num.Rec. avec tableaux commencant a 1. */
{
        float *arr; int_4 *indx;
	int_4 i,indxt,ir=n,itemp,j,k,l=1;
	int_4 jstack=0;
	float a;
        int_4 istack[NSTACK+1];

        arr = arr_c-1;
        indx = indx_c-1;

	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=1;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP_INDEXR4(indx[k],indx[l+1]);
			if (arr[indx[l+1]] > arr[indx[ir]]) {
				SWAP_INDEXR4(indx[l+1],indx[ir])
			}
			if (arr[indx[l]] > arr[indx[ir]]) {
				SWAP_INDEXR4(indx[l],indx[ir])
			}
			if (arr[indx[l+1]] > arr[indx[l]]) {
				SWAP_INDEXR4(indx[l+1],indx[l])
			}
			i=l+1;
			j=ir;
			indxt=indx[l];
			a=arr[indxt];
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;
				SWAP_INDEXR4(indx[i],indx[j])
			}
			indx[l]=indx[j];
			indx[j]=indxt;
			jstack += 2;
			if (jstack > NSTACK) {
                          printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
                          exit(-1);
			}
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
        /* pour avoir un retour d'indice entre 0 et n-1 */
	for (j=0;j<n;j++) indx_c[j]=indx_c[j]-1;
}
#undef SWAP_INDEXR4
#undef M
#undef NSTACK

/*===================================================================================*/
/*
++
  void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
	Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
	such that arr[indx[j]] is in ascending order for j=1,2,,N.
	The input quantities n and arr are not changed.
--
*/
#define SWAP_INDEXR8(a,b) itemp=(a);(a)=(b);(b)=itemp;
#define M 7
#define NSTACK 50
void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
/* encore du Num.Rec. avec tableaux commencant a 1. */
{
        double *arr; int_4 *indx;
	int_4 i,indxt,ir=n,itemp,j,k,l=1;
	int_4 jstack=0;
	double a;
        int_4 istack[NSTACK+1];

        arr = arr_c-1;
        indx = indx_c-1;

	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=1;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			SWAP_INDEXR8(indx[k],indx[l+1]);
			if (arr[indx[l+1]] > arr[indx[ir]]) {
				SWAP_INDEXR8(indx[l+1],indx[ir])
			}
			if (arr[indx[l]] > arr[indx[ir]]) {
				SWAP_INDEXR8(indx[l],indx[ir])
			}
			if (arr[indx[l+1]] > arr[indx[l]]) {
				SWAP_INDEXR8(indx[l+1],indx[l])
			}
			i=l+1;
			j=ir;
			indxt=indx[l];
			a=arr[indxt];
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;
				SWAP_INDEXR8(indx[i],indx[j])
			}
			indx[l]=indx[j];
			indx[j]=indxt;
			jstack += 2;
			if (jstack > NSTACK) {
                          printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
                          exit(-1);
			}
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
        /* pour avoir un retour d'indice entre 0 et n-1 */
	for (j=0;j<n;j++) indx_c[j]=indx_c[j]-1;
}
#undef SWAP_INDEXR8
#undef M
#undef NSTACK