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

Last change on this file since 3223 was 658, checked in by ansari, 26 years ago

no message
File size: 21.7 KB
RevLine 
[658]1#include <unistd.h>
2#include <stdlib.h>
3#include <stdio.h>
4#include <math.h>
5#include "nbtri.h"
6
7/*
8++
9 Module Tri de tableaux (C)
10 Lib LibsUtil
11 include nbtri.h
12
13 Tri de tableaux et indexation.
14--
15*/
16
17/*=========================================================================*/
18/*
19++
20 void HeapSort(int n,double *ra_int)
21 On reordonne par ordre numerique croissant
22 le tableau ra_int[n]: Numerical Recipes mode C.
23--
24*/
25void HeapSort(int n,double *ra_int)
26{
27 int l,j,ir,i;
28 double rra,*ra;
29
30 /* attention, Numerical recipes prend des tableaux de 1 a n on remet
31 de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
32 ra = ra_int-1;
33
34 l=(n >> 1)+1;
35 ir=n;
36 for (;;) {
37 if (l > 1)
38 rra=ra[--l];
39 else {
40 rra=ra[ir];
41 ra[ir]=ra[1];
42 if (--ir == 1) {
43 ra[1]=rra;
44 return;
45 }
46 }
47 i=l;
48 j=l << 1;
49 while (j <= ir) {
50 if (j < ir && ra[j] < ra[j+1]) ++j;
51 if (rra < ra[j]) {
52 ra[i]=ra[j];
53 j += (i=j);
54 }
55 else j=ir+1;
56 }
57 ra[i]=rra;
58 }
59}
60
61/*=========================================================================*/
62/*
63++
64 void HeapSortF(int n,float *ra_int)
65 On reordonne par ordre numerique croissant
66 le tableau ra_int[n]: Numerical Recipes mode C.
67--
68*/
69void HeapSortF(int n,float *ra_int)
70{
71 int l,j,ir,i;
72 float rra,*ra;
73
74 /* attention, Numerical reciepes prend des tableaux de 1 a n on remet
75 de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
76 ra = ra_int-1;
77
78 l=(n >> 1)+1;
79 ir=n;
80 for (;;) {
81 if (l > 1)
82 rra=ra[--l];
83 else {
84 rra=ra[ir];
85 ra[ir]=ra[1];
86 if (--ir == 1) {
87 ra[1]=rra;
88 return;
89 }
90 }
91 i=l;
92 j=l << 1;
93 while (j <= ir) {
94 if (j < ir && ra[j] < ra[j+1]) ++j;
95 if (rra < ra[j]) {
96 ra[i]=ra[j];
97 j += (i=j);
98 }
99 else j=ir+1;
100 }
101 ra[i]=rra;
102 }
103}
104
105/*=========================================================================*/
106/*
107++
108 void HeapSortF2(int n,float *ra_int,float *ra2_int)
109 On reordonne par ordre numerique croissant
110 le tableau ra_int[n] et ra2_int en parralelle.
111 Numerical Recipes mode C.
112--
113*/
114void HeapSortF2(int n,float *ra_int,float *ra2_int)
115{
116 int l,j,ir,i;
117 float rra,rra2,*ra,*ra2;
118
119 /* attention, Numerical reciepes prend des tableaux de 1 a n on remet
120 de 0 a n-1 en decramentant le pointeur du tableau d'entree*/
121 ra = ra_int-1;
122 ra2 = ra2_int-1;
123
124 l=(n >> 1)+1;
125 ir=n;
126 for (;;) {
127 if (l > 1)
128 {rra=ra[--l]; rra2=ra2[l];}
129 else {
130 rra=ra[ir]; rra2=ra2[ir];
131 ra[ir]=ra[1]; ra2[ir]=ra2[1];
132 if (--ir == 1) {
133 ra[1]=rra; ra2[1]=rra2;
134 return;
135 }
136 }
137 i=l;
138 j=l << 1;
139 while (j <= ir) {
140 if (j < ir && ra[j] < ra[j+1]) ++j;
141 if (rra < ra[j]) {
142 ra[i]=ra[j]; ra2[i]=ra2[j];
143 j += (i=j);
144 }
145 else j=ir+1;
146 }
147 ra[i]=rra; ra2[i]=rra2;
148 }
149}
150
151/*=========================================================================*/
152/*
153++
154 int_4 tri_double ( double *tab, int_4 *indx,int_4 N)
155 Methode de tri sans finesse (double boucles).
156| entree : tab -> tableau de double de longueur N (0 a N-1 )
157| sortie : indx -> tableau d'entiers de longueur N : la case i
158| contient le classement du ieme element dans tab
159| Retourne: 0 si tri impossible
160| 1 si tri reussi
161--
162*/
163int_4 tri_double ( double *tab, int_4 *indx,int_4 N)
164{
165int_4 i,j,k;
166
167if (N<=0) return (-1);
168
169for (i=0; i<N ; i++) indx[i]=i;
170if (N==1) return (N);
171
172/* classement dans l'ordre croissant */
173for (i=1; i< N ; i++)
174 {
175
176 if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
177 {
178 k = indx[i-1];
179 indx[i-1] = indx[i];
180 indx[i] = k;
181
182 if ( i > 1 )
183 for (j=i ; j>=1 ; j--)
184 { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
185 {
186 k = indx[j-1];
187 indx[j-1] = indx[j];
188 indx[j] = k;
189 } } } }
190
191return (N) ;
192}
193
194/*=========================================================================*/
195/*
196++
197 int_4 tri_float ( float *tab, int_4 *indx,int_4 N)
198 Methode de tri sans finesse (double boucles).
199| entree : tab -> tableau de flottant de longueur N (0 a N-1 )
200| sortie : indx -> tableau d'entiers de longueur N : la case i
201| contient le classement du ieme element dans tab
202| Retourne: 0 si tri impossible
203| 1 si tri reussi
204--
205*/
206int_4 tri_float ( float *tab, int_4 *indx,int_4 N)
207{
208int_4 i,j,k;
209
210if (N<=0) return (-1);
211
212for (i=0; i<N ; i++) indx[i]=i;
213if (N==1) return (N);
214
215/* classement dans l'ordre croissant */
216for (i=1; i< N ; i++)
217 {
218
219 if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
220 {
221 k = indx[i-1];
222 indx[i-1] = indx[i];
223 indx[i] = k;
224
225 if ( i > 1 )
226 for (j=i ; j>=1 ; j--)
227 { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
228 {
229 k = indx[j-1];
230 indx[j-1] = indx[j];
231 indx[j] = k;
232 } } } }
233
234return (N) ;
235}
236
237/*=========================================================================*/
238/*
239++
240 int_4 tri_entier ( int_4 *tab,int_4 *indx,int_4 N)
241 Methode de tri sans finesse (double boucles).
242| entree : tab -> tableau d'entiers de longueur N (0 a N-1 )
243| sortie : indx -> tableau d'entiers de longueur N : la case i
244| contient le classement du ieme element dans tab
245| Retourne: 0 si tri impossible
246| 1 si tri reussi
247--
248*/
249int_4 tri_entier ( int_4 *tab,int_4 *indx,int_4 N)
250{
251int_4 i,j,k;
252
253if (N<=0) return (-1);
254
255for (i=0; i<N ; i++) indx[i]=i;
256if (N==1) return (N);
257
258/* classement dans l'ordre croissant */
259for (i=1; i< N ; i++)
260 {
261
262 if ( *(tab+indx[i]) < *(tab+indx[i-1]) )
263 {
264 k = indx[i-1];
265 indx[i-1] = indx[i];
266 indx[i] = k;
267
268 if ( i > 1 )
269 for (j=i ; j>=1 ; j--)
270 { if ( *(tab+indx[j]) < *(tab+indx[j-1]) )
271 {
272 k = indx[j-1];
273 indx[j-1] = indx[j];
274 indx[j] = k;
275 } } } }
276
277return (N) ;
278}
279
280/*=========================================================================*/
281/*
282++
283 int_4 tri_rapide_I (int_4 *datum,int_4 *index,int_4 N)
284 algorythme de tri rapide
285 ( p114-119 THE ART OF COMPUTER PROGRAMMING, vol. 3 SORTING AND
286 SEARCHING de D.E. Knuth)
287| datum ( entree/sortie ) -> vecteur de dimension N contenant des ENTIERS
288| desordonnes en entree. Ces elements ressortent ordonnes.
289| index ( sortie ) -> vecteur de dimension N contenant des entiers.
290| Le contenu de la case i de index nous indique la place
291| du ieme element du datum originel.
292| tri_rapide = -1 si echec, 1 si reussite de l'operation
293--
294*/
295int_4 tri_rapide_I (int_4 *datum,int_4 *index,int_4 N)
296{
297/* 14 = nb max d'entrees que la pile (stack) peut contenir. Ce 14 limite la longueur
298 maximale possible pour les vecteurs a 32 768 ( = 2**15 ) */
299
300int_4 stklo[14], stkhi[14], hi, nstak, i, limlo, limhi, lo, ikey;
301int_4 dkey;
302
303/* initialisations */
304
305for (i=0; i<=N-1 ;i++) index[i]=i;
306
307nstak=0;
308limlo=0; limhi=N-1;
309
310grande_boucle:
311 dkey = *(datum+limlo);
312 ikey = *(index+limlo);
313
314/* compare tous les elements d'un ss-vecteur entre limlo et limhi avec la donnee-cle courante */
315
316 lo=limlo; hi=limhi;
317
318sous_boucle_1:
319 if ( lo == hi ) goto lo_egal_hi;
320
321 if ( *(datum+hi) <= dkey ) goto remplacement_lo;
322 hi = hi - 1;
323
324/* le pointeur hi doit pointer une donnee plus petite que la cle qui va etre remplacer */
325
326 goto sous_boucle_1;
327
328remplacement_lo:
329 *(datum+lo) = *(datum+hi);
330 *(index+lo) = *(index+hi);
331 lo = lo + 1;
332
333sous_boucle_2:
334 if ( lo == hi ) goto lo_egal_hi;
335
336 if ( *(datum+lo) >= dkey ) goto remplacement_hi;
337
338 lo = lo + 1;
339 goto sous_boucle_2;
340
341remplacement_hi:
342 *(datum+hi) = *(datum+lo);
343 *(index+hi) = *(index+lo);
344 hi = hi - 1;
345
346/* le pointeur lo doit pointer une donnee plus grande que la cle qui va etre remplacer */
347
348 goto sous_boucle_1;
349
350lo_egal_hi:
351
352/* 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 */
353
354 *(datum+lo) = dkey;
355 *(index+lo) = ikey;
356
357/* 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)
358 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. */
359
360 if ( (limhi-lo) > (lo-limlo) ) goto cas_1;
361
362/* 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 */
363
364 if ( (lo-limlo) <= 1) goto test_fin;
365
366/* 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. */
367
368 if ( (limhi-lo) >= 2) goto cas_1b;
369
370/* 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) */
371
372 limhi=lo-1;
373 goto grande_boucle;
374
375/* 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. */
376
377cas_1b:
378 nstak=nstak+1;
379 *(stklo+nstak)=limlo;
380 *(stkhi+nstak)=lo-1;
381 limlo=lo+1;
382 goto grande_boucle;
383
384/* 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. */
385
386cas_1:
387 if ( (limhi-lo) <= 1) goto test_fin;
388
389/* 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. */
390
391 if ( (lo-limlo) >= 2) goto cas_2b;
392
393/* 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) */
394
395 limlo=lo+1;
396 goto grande_boucle;
397
398/* 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. */
399
400cas_2b:
401 nstak=nstak+1;
402 *(stklo+nstak)=lo+1;
403 *(stkhi+nstak)=limhi;
404 limlo=lo-1;
405 goto grande_boucle;
406
407/* on prend l'intervalle le plus recent de la pile. Si le stack est vide, c'est fini !!! */
408
409test_fin:
410 if (nstak <= 0) return (1);
411 limlo = *(stklo+nstak);
412 limhi = *(stkhi+nstak);
413 nstak=nstak-1;
414 goto grande_boucle;
415
416}
417
418/*=========================================================================*/
419/*
420++
421 int_4 tri_rapide_F (float *datum,int_4 *index,int_4 N)
422 Idem tri_rapide_I mais pour un tableau de flottants
423 REMPLI avec des ENTIERS.
424--
425*/
426int_4 tri_rapide_F (float *datum,int_4 *index,int_4 N)
427{
428/* ATTENTION, Ce programme tri un tableau de flottants REMPLI avec des ENTIERS!!!!!
429 14 = nb max d'entrees que la pile (stack) peut contenir. Ce 14 limite la longueur
430 maximale possible pour les vecteurs a 32 768 ( = 2**15 ) */
431
432int_4 stklo[14], stkhi[14], hi, nstak, i, limlo, limhi, lo, ikey;
433float dkey;
434
435/* initialisations */
436
437for (i=0; i<=N-1 ; i++) index[i]=i;
438nstak=0;
439limlo=0; limhi=N-1;
440
441grande_boucle:
442 dkey = *(datum+limlo);
443 ikey = *(index+limlo);
444
445/* compare tous les elements d'un ss-vecteur entre limlo et limhi avec la donnee-cle courante */
446
447 lo=limlo; hi=limhi;
448
449sous_boucle_1:
450 if ( lo == hi ) goto lo_egal_hi;
451
452 if ( *(datum+hi) <= dkey ) goto remplacement_lo;
453 hi = hi - 1;
454
455/* le pointeur hi doit pointer une donnee plus petite que la cle qui va etre remplacer */
456
457 goto sous_boucle_1;
458
459remplacement_lo:
460 *(datum+lo) = *(datum+hi);
461 *(index+lo) = *(index+hi);
462 lo = lo + 1;
463
464sous_boucle_2:
465 if ( lo == hi ) goto lo_egal_hi;
466
467 if ( *(datum+lo) >= dkey ) goto remplacement_hi;
468
469 lo = lo + 1;
470 goto sous_boucle_2;
471
472remplacement_hi:
473 *(datum+hi) = *(datum+lo);
474 *(index+hi) = *(index+lo);
475 hi = hi - 1;
476
477/* le pointeur lo doit pointer une donnee plus grande que la cle qui va etre remplacer */
478
479 goto sous_boucle_1;
480
481lo_egal_hi:
482
483/* 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 */
484
485 *(datum+lo) = dkey;
486 *(index+lo) = ikey;
487
488/* 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)
489 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. */
490
491 if ( (limhi-lo) > (lo-limlo) ) goto cas_1;
492
493/* 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 */
494
495 if ( (lo-limlo) <= 1) goto test_fin;
496
497/* 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. */
498
499 if ( (limhi-lo) >= 2) goto cas_1b;
500
501/* 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) */
502
503 limhi=lo-1;
504 goto grande_boucle;
505
506/* 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. */
507
508cas_1b:
509 nstak=nstak+1;
510 *(stklo+nstak)=limlo;
511 *(stkhi+nstak)=lo-1;
512 limlo=lo+1;
513 goto grande_boucle;
514
515/* 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. */
516
517cas_1:
518 if ( (limhi-lo) <= 1) goto test_fin;
519
520/* 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. */
521
522 if ( (lo-limlo) >= 2) goto cas_2b;
523
524/* 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) */
525
526 limlo=lo+1;
527 goto grande_boucle;
528
529/* 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. */
530
531cas_2b:
532 nstak=nstak+1;
533 *(stklo+nstak)=lo+1;
534 *(stkhi+nstak)=limhi;
535 limlo=lo-1;
536 goto grande_boucle;
537
538/* on prend l'intervalle le plus recent de la pile. Si le stack est vide, c'est fini !!! */
539
540test_fin:
541 if (nstak <= 0) return (1);
542 limlo = *(stklo+nstak);
543 limhi = *(stkhi+nstak);
544 nstak=nstak-1;
545 goto grande_boucle;
546
547}
548
549/*===================================================================================*/
550/*
551++
552 int qSort_Float(const void *a1,const void *a2)
553 Fonction de tri de `float' a utiliser dans qsort.
554--
555*/
556int qSort_Float(const void *a1,const void *a2)
557{
558if( *((float *) a1) < *((float *) a2) ) return(-1);
559if( *((float *) a1) > *((float *) a2) ) return( 1);
560return(0);
561}
562
563/*===================================================================================*/
564/*
565++
566 int qSort_Dble(const void *a1,const void *a2)
567 Fonction de tri de `double' a utiliser dans qsort.
568--
569*/
570int qSort_Dble(const void *a1,const void *a2)
571{
572if( *((double *) a1) < *((double *) a2) ) return(-1);
573if( *((double *) a1) > *((double *) a2) ) return( 1);
574return(0);
575}
576
577/*===================================================================================*/
578/*
579++
580 int qSort_Int(const void *a1,const void *a2)
581 Fonction de tri de `int' a utiliser dans qsort.
582--
583*/
584int qSort_Int(const void *a1,const void *a2)
585{
586if( *((int *) a1) < *((int *) a2) ) return(-1);
587if( *((int *) a1) > *((int *) a2) ) return( 1);
588return(0);
589}
590
591/*===================================================================================*/
592/*
593++
594 int qSort_Ushort(const void *a1,const void *a2)
595 Fonction de tri de `unsigned short' a utiliser dans qsort.
596--
597*/
598int qSort_Ushort(const void *a1,const void *a2)
599{
600if( *((unsigned short *) a1) < *((unsigned short *) a2) ) return(-1);
601if( *((unsigned short *) a1) > *((unsigned short *) a2) ) return( 1);
602return(0);
603}
604
605/*===================================================================================*/
606/*
607++
608 int qSort_Short(const void *a1,const void *a2)
609 Fonction de tri de `short' a utiliser dans qsort.
610--
611*/
612int qSort_Short(const void *a1,const void *a2)
613{
614if( *((short *) a1) < *((short *) a2) ) return(-1);
615if( *((short *) a1) > *((short *) a2) ) return( 1);
616return(0);
617}
618
619/*===================================================================================*/
620/*
621++
622 void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
623 Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
624 such that arr[indx[j]] is in ascending order for j=1,2,,N.
625 The input quantities n and arr are not changed.
626--
627*/
628#define SWAP_INDEXR4(a,b) itemp=(a);(a)=(b);(b)=itemp;
629#define M 7
630#define NSTACK 50
631void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
632/* encore du Num.Rec. avec tableaux commencant a 1. */
633{
634 float *arr; int_4 *indx;
635 int_4 i,indxt,ir=n,itemp,j,k,l=1;
636 int_4 jstack=0;
637 float a;
638 int_4 istack[NSTACK+1];
639
640 arr = arr_c-1;
641 indx = indx_c-1;
642
643 for (j=1;j<=n;j++) indx[j]=j;
644 for (;;) {
645 if (ir-l < M) {
646 for (j=l+1;j<=ir;j++) {
647 indxt=indx[j];
648 a=arr[indxt];
649 for (i=j-1;i>=1;i--) {
650 if (arr[indx[i]] <= a) break;
651 indx[i+1]=indx[i];
652 }
653 indx[i+1]=indxt;
654 }
655 if (jstack == 0) break;
656 ir=istack[jstack--];
657 l=istack[jstack--];
658 } else {
659 k=(l+ir) >> 1;
660 SWAP_INDEXR4(indx[k],indx[l+1]);
661 if (arr[indx[l+1]] > arr[indx[ir]]) {
662 SWAP_INDEXR4(indx[l+1],indx[ir])
663 }
664 if (arr[indx[l]] > arr[indx[ir]]) {
665 SWAP_INDEXR4(indx[l],indx[ir])
666 }
667 if (arr[indx[l+1]] > arr[indx[l]]) {
668 SWAP_INDEXR4(indx[l+1],indx[l])
669 }
670 i=l+1;
671 j=ir;
672 indxt=indx[l];
673 a=arr[indxt];
674 for (;;) {
675 do i++; while (arr[indx[i]] < a);
676 do j--; while (arr[indx[j]] > a);
677 if (j < i) break;
678 SWAP_INDEXR4(indx[i],indx[j])
679 }
680 indx[l]=indx[j];
681 indx[j]=indxt;
682 jstack += 2;
683 if (jstack > NSTACK) {
684 printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
685 exit(-1);
686 }
687 if (ir-i+1 >= j-l) {
688 istack[jstack]=ir;
689 istack[jstack-1]=i;
690 ir=j-1;
691 } else {
692 istack[jstack]=j-1;
693 istack[jstack-1]=l;
694 l=i;
695 }
696 }
697 }
698}
699#undef SWAP_INDEXR4
700#undef M
701#undef NSTACK
702
703/*===================================================================================*/
704/*
705++
706 void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
707 Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
708 such that arr[indx[j]] is in ascending order for j=1,2,,N.
709 The input quantities n and arr are not changed.
710--
711*/
712#define SWAP_INDEXR8(a,b) itemp=(a);(a)=(b);(b)=itemp;
713#define M 7
714#define NSTACK 50
715void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
716/* encore du Num.Rec. avec tableaux commencant a 1. */
717{
718 double *arr; int_4 *indx;
719 int_4 i,indxt,ir=n,itemp,j,k,l=1;
720 int_4 jstack=0;
721 double a;
722 int_4 istack[NSTACK+1];
723
724 arr = arr_c-1;
725 indx = indx_c-1;
726
727 for (j=1;j<=n;j++) indx[j]=j;
728 for (;;) {
729 if (ir-l < M) {
730 for (j=l+1;j<=ir;j++) {
731 indxt=indx[j];
732 a=arr[indxt];
733 for (i=j-1;i>=1;i--) {
734 if (arr[indx[i]] <= a) break;
735 indx[i+1]=indx[i];
736 }
737 indx[i+1]=indxt;
738 }
739 if (jstack == 0) break;
740 ir=istack[jstack--];
741 l=istack[jstack--];
742 } else {
743 k=(l+ir) >> 1;
744 SWAP_INDEXR8(indx[k],indx[l+1]);
745 if (arr[indx[l+1]] > arr[indx[ir]]) {
746 SWAP_INDEXR8(indx[l+1],indx[ir])
747 }
748 if (arr[indx[l]] > arr[indx[ir]]) {
749 SWAP_INDEXR8(indx[l],indx[ir])
750 }
751 if (arr[indx[l+1]] > arr[indx[l]]) {
752 SWAP_INDEXR8(indx[l+1],indx[l])
753 }
754 i=l+1;
755 j=ir;
756 indxt=indx[l];
757 a=arr[indxt];
758 for (;;) {
759 do i++; while (arr[indx[i]] < a);
760 do j--; while (arr[indx[j]] > a);
761 if (j < i) break;
762 SWAP_INDEXR8(indx[i],indx[j])
763 }
764 indx[l]=indx[j];
765 indx[j]=indxt;
766 jstack += 2;
767 if (jstack > NSTACK) {
768 printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
769 exit(-1);
770 }
771 if (ir-i+1 >= j-l) {
772 istack[jstack]=ir;
773 istack[jstack-1]=i;
774 ir=j-1;
775 } else {
776 istack[jstack]=j-1;
777 istack[jstack-1]=l;
778 l=i;
779 }
780 }
781 }
782}
783#undef SWAP_INDEXR8
784#undef M
785#undef NSTACK
Note: See TracBrowser for help on using the repository browser.