#include #include #include #include #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 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 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 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