source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/nbtri.c

#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_4 tri_rapide_I (int_4 *datum,int_4 *index,int_4 N)
        algorythme de tri rapide
        ( p114-119 THE ART OF COMPUTER PROGRAMMING, vol. 3 SORTING AND
        SEARCHING de D.E. Knuth)
| datum  ( entree/sortie ) -> vecteur de dimension N contenant des ENTIERS
|        desordonnes en entree. Ces elements ressortent ordonnes.
| index  ( sortie ) -> vecteur de dimension N contenant des entiers.
|        Le contenu de la case i de index nous indique la place
|        du ieme element du datum originel.
| tri_rapide = -1 si echec, 1 si reussite de l'operation
--
*/
int_4 tri_rapide_I (int_4 *datum,int_4 *index,int_4 N)
{
/* 14 = nb max d'entrees que la pile (stack) peut contenir. Ce 14 limite la longueur 
             maximale possible pour les vecteurs a 32 768 ( = 2**15 ) */

int_4 stklo[14], stkhi[14], hi, nstak, i, limlo, limhi, lo, ikey;
int_4 dkey;

/* initialisations */

for (i=0; i<=N-1 ;i++)  index[i]=i;

nstak=0;
limlo=0; limhi=N-1;

grande_boucle:
  dkey = *(datum+limlo);
  ikey = *(index+limlo);

/* compare tous les elements d'un ss-vecteur entre limlo et limhi avec la donnee-cle courante */

  lo=limlo; hi=limhi;
  
sous_boucle_1:
  if ( lo == hi ) goto lo_egal_hi;

  if ( *(datum+hi) <= dkey ) goto remplacement_lo;
  hi = hi - 1;

/* le pointeur  hi doit pointer une donnee plus petite que la cle qui va etre remplacer */

  goto sous_boucle_1;

remplacement_lo:
  *(datum+lo) = *(datum+hi);
  *(index+lo) = *(index+hi);
  lo = lo + 1;
  
sous_boucle_2:
  if ( lo == hi ) goto lo_egal_hi;

  if ( *(datum+lo) >= dkey ) goto remplacement_hi;

  lo = lo + 1;
  goto sous_boucle_2;

remplacement_hi:
  *(datum+hi) = *(datum+lo);
  *(index+hi) = *(index+lo);
  hi = hi - 1;

/* le pointeur lo doit pointer une donnee plus grande que la cle qui va etre remplacer */

  goto sous_boucle_1;

lo_egal_hi:

/* lo et hi sont egaux, et pointent sur une valeur qui va etre remplacer. Tant que toutes les valeurs sous ce point sont inferieures a la cle et que toutes les valeurs apres ce point sont superieures a la cle, c'est la que nous remettrons la cle dans le vecteur */

  *(datum+lo) = dkey;
  *(index+lo) = ikey;

/* A ce point du ss-prog. toutes les donnees entre limlo et lo-1 inclus sont inferieurs a datum(lo), et toutes les donnes entre lo+1 et limhi sont superieures a datum(lo)  
   Si les deux ss-tableaux ne contiennent pas plus d'un element, on prend l'intervale le plus recent de la pile ( si le stack est vide c'est fini ). Si le plus grand des deux ss-tableaux contient plus d'un element et si le plus petit contient au plus un element, alors on oublie le plus petit et on reduit l'autre. Si le plus petit ss-tableau contient au moins 2 elements, alors on place le plus grand ss-tableau dans la pile et on processe le ss-tableau.  */
  
  if ( (limhi-lo) > (lo-limlo) ) goto cas_1;

/* cas 1 : le ss-tableau inferieur est plus long. Si il contient un element au plus, alors on prend l'intervalle du stack le plus recent, on le ramene et on travaille dessus */

  if ( (lo-limlo) <= 1) goto test_fin;

/* si le ss-tableau superieur ( le + court ss-tableau ) contient  un element au plus alors on processe le  ss-tableau inferieur ( le + long), mais si le ss-tableau superieur contient plus d'un element, alors on place le  ss-tableau inferieur ( le + long) sur la pile et on processe le ss-tableau superieur. */

  if ( (limhi-lo) >= 2) goto cas_1b;

/* cas 1a : si le ss-tableau superieur ( le + court ss-tableau ) contient  un element au plus alors on revient en arriere et on agit sur le  ss-tableau inferieur ( le + long) */

  limhi=lo-1;
  goto grande_boucle;

/* cas 1b :  si le ss-tableau superieur ( le + court ss-tableau ) contient plus d'un element, alors on place le  ss-tableau inferieur ( le + long) sur la pile et on processe le ss-tableau superieur. */

cas_1b:
  nstak=nstak+1;
  *(stklo+nstak)=limlo;
  *(stkhi+nstak)=lo-1;
  limlo=lo+1;
  goto grande_boucle;
 
/* cas 2 : le ss-tableau superieur est plus long, si il contient  un element au plus alors on agit sur l'intervalle le plus recent de la pile. */

cas_1:
  if ( (limhi-lo) <= 1) goto test_fin;

/* si le ss-tableau inferieur ( le + court ss-tableau ) contient  un element au plus alors on processe le  ss-tableau superieur ( le + long), mais si le ss-tableau inferieur contient plus d'un element, alors on place le  ss-tableau superieur ( le + long) sur la pile et on processe le ss-tableau inferieur. */

   if ( (lo-limlo) >= 2) goto cas_2b;

/*  cas 2a : si le ss-tableau inferieur ( le + court ss-tableau ) contient  un element au plus alors on revient en arriere et on agit sur le  ss-tableau superieur ( le + long) */

  limlo=lo+1;
  goto grande_boucle;
  
/* cas 2b :  si le ss-tableau inferieur ( le + court ss-tableau ) contient plus d'un element, alors on place le  ss-tableau superieur ( le + long) sur la pile et on processe le ss-tableau inferieur. */

cas_2b:
  nstak=nstak+1;
  *(stklo+nstak)=lo+1;
  *(stkhi+nstak)=limhi;
  limlo=lo-1;
  goto grande_boucle;
   
/* on prend l'intervalle le plus recent de la pile. Si le stack est vide, c'est fini !!! */
        
test_fin:
  if (nstak <= 0) return (1);
  limlo = *(stklo+nstak);
  limhi = *(stkhi+nstak);
  nstak=nstak-1;
  goto grande_boucle;

}

/*=========================================================================*/
/*
++
  int_4 tri_rapide_F (float *datum,int_4 *index,int_4 N)
        Idem tri_rapide_I mais pour un tableau de flottants
         REMPLI avec des ENTIERS.
--
*/
int_4 tri_rapide_F (float *datum,int_4 *index,int_4 N)
{
/* ATTENTION, Ce programme tri un tableau de flottants REMPLI avec des ENTIERS!!!!!
  14 = nb max d'entrees que la pile (stack) peut contenir. Ce 14 limite la longueur 
             maximale possible pour les vecteurs a 32 768 ( = 2**15 ) */

int_4 stklo[14], stkhi[14], hi, nstak, i, limlo, limhi, lo, ikey;
float dkey;

/* initialisations */

for (i=0; i<=N-1 ; i++) index[i]=i;
nstak=0;
limlo=0; limhi=N-1;

grande_boucle:
  dkey = *(datum+limlo);
  ikey = *(index+limlo);

/* compare tous les elements d'un ss-vecteur entre limlo et limhi avec la donnee-cle courante */

  lo=limlo; hi=limhi;
  
sous_boucle_1:
  if ( lo == hi ) goto lo_egal_hi;

  if ( *(datum+hi) <= dkey ) goto remplacement_lo;
  hi = hi - 1;

/* le pointeur  hi doit pointer une donnee plus petite que la cle qui va etre remplacer */

  goto sous_boucle_1;

remplacement_lo:
  *(datum+lo) = *(datum+hi);
  *(index+lo) = *(index+hi);
  lo = lo + 1;
  
sous_boucle_2:
  if ( lo == hi ) goto lo_egal_hi;

  if ( *(datum+lo) >= dkey ) goto remplacement_hi;

  lo = lo + 1;
  goto sous_boucle_2;

remplacement_hi:
  *(datum+hi) = *(datum+lo);
  *(index+hi) = *(index+lo);
  hi = hi - 1;

/* le pointeur lo doit pointer une donnee plus grande que la cle qui va etre remplacer */

  goto sous_boucle_1;

lo_egal_hi:

/* lo et hi sont egaux, et pointent sur une valeur qui va etre remplacer. Tant que toutes les valeurs sous ce point sont inferieures a la cle et que toutes les valeurs apres ce point sont superieures a la cle, c'est la que nous remettrons la cle dans le vecteur */

  *(datum+lo) = dkey;
  *(index+lo) = ikey;

/* A ce point du ss-prog. toutes les donnees entre limlo et lo-1 inclus sont inferieurs a datum(lo), et toutes les donnes entre lo+1 et limhi sont superieures a datum(lo)  
   Si les deux ss-tableaux ne contiennent pas plus d'un element, on prend l'intervale le plus recent de la pile ( si le stack est vide c'est fini ). Si le plus grand des deux ss-tableaux contient plus d'un element et si le plus petit contient au plus un element, alors on oublie le plus petit et on reduit l'autre. Si le plus petit ss-tableau contient au moins 2 elements, alors on place le plus grand ss-tableau dans la pile et on processe le ss-tableau.  */
  
  if ( (limhi-lo) > (lo-limlo) ) goto cas_1;

/* cas 1 : le ss-tableau inferieur est plus long. Si il contient un element au plus, alors on prend l'intervalle du stack le plus recent, on le ramene et on travaille dessus */

  if ( (lo-limlo) <= 1) goto test_fin;

/* si le ss-tableau superieur ( le + court ss-tableau ) contient  un element au plus alors on processe le  ss-tableau inferieur ( le + long), mais si le ss-tableau superieur contient plus d'un element, alors on place le  ss-tableau inferieur ( le + long) sur la pile et on processe le ss-tableau superieur. */

  if ( (limhi-lo) >= 2) goto cas_1b;

/* cas 1a : si le ss-tableau superieur ( le + court ss-tableau ) contient  un element au plus alors on revient en arriere et on agit sur le  ss-tableau inferieur ( le + long) */

  limhi=lo-1;
  goto grande_boucle;

/* cas 1b :  si le ss-tableau superieur ( le + court ss-tableau ) contient plus d'un element, alors on place le  ss-tableau inferieur ( le + long) sur la pile et on processe le ss-tableau superieur. */

cas_1b:
  nstak=nstak+1;
  *(stklo+nstak)=limlo;
  *(stkhi+nstak)=lo-1;
  limlo=lo+1;
  goto grande_boucle;
 
/* cas 2 : le ss-tableau superieur est plus long, si il contient  un element au plus alors on agit sur l'intervalle le plus recent de la pile. */

cas_1:
  if ( (limhi-lo) <= 1) goto test_fin;

/* si le ss-tableau inferieur ( le + court ss-tableau ) contient  un element au plus alors on processe le  ss-tableau superieur ( le + long), mais si le ss-tableau inferieur contient plus d'un element, alors on place le  ss-tableau superieur ( le + long) sur la pile et on processe le ss-tableau inferieur. */

   if ( (lo-limlo) >= 2) goto cas_2b;

/*  cas 2a : si le ss-tableau inferieur ( le + court ss-tableau ) contient  un element au plus alors on revient en arriere et on agit sur le  ss-tableau superieur ( le + long) */

  limlo=lo+1;
  goto grande_boucle;
  
/* cas 2b :  si le ss-tableau inferieur ( le + court ss-tableau ) contient plus d'un element, alors on place le  ss-tableau superieur ( le + long) sur la pile et on processe le ss-tableau inferieur. */

cas_2b:
  nstak=nstak+1;
  *(stklo+nstak)=lo+1;
  *(stkhi+nstak)=limhi;
  limlo=lo-1;
  goto grande_boucle;
   
/* on prend l'intervalle le plus recent de la pile. Si le stack est vide, c'est fini !!! */
        
test_fin:
  if (nstak <= 0) return (1);
  limlo = *(stklo+nstak);
  limhi = *(stkhi+nstak);
  nstak=nstak-1;
  goto grande_boucle;

}

/*===================================================================================*/
/*
++
  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;
                        }
                }
        }
}
#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;
                        }
                }
        }
}
#undef SWAP_INDEXR8
#undef M
#undef NSTACK