source: Sophya/trunk/SophyaLib/NTools/nbtri.c

Last change on this file was 1474, checked in by cmv, 24 years ago

2 routines de tri bugguee enlevees cmv 19/4/2001
File size: 11.8 KB
RevLine 
[220]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 qSort_Float(const void *a1,const void *a2)
284 Fonction de tri de `float' a utiliser dans qsort.
285--
286*/
287int qSort_Float(const void *a1,const void *a2)
288{
289if( *((float *) a1) < *((float *) a2) ) return(-1);
290if( *((float *) a1) > *((float *) a2) ) return( 1);
291return(0);
292}
293
294/*===================================================================================*/
295/*
296++
297 int qSort_Dble(const void *a1,const void *a2)
298 Fonction de tri de `double' a utiliser dans qsort.
299--
300*/
301int qSort_Dble(const void *a1,const void *a2)
302{
303if( *((double *) a1) < *((double *) a2) ) return(-1);
304if( *((double *) a1) > *((double *) a2) ) return( 1);
305return(0);
306}
307
308/*===================================================================================*/
309/*
310++
311 int qSort_Int(const void *a1,const void *a2)
312 Fonction de tri de `int' a utiliser dans qsort.
313--
314*/
315int qSort_Int(const void *a1,const void *a2)
316{
317if( *((int *) a1) < *((int *) a2) ) return(-1);
318if( *((int *) a1) > *((int *) a2) ) return( 1);
319return(0);
320}
321
322/*===================================================================================*/
323/*
324++
325 int qSort_Ushort(const void *a1,const void *a2)
326 Fonction de tri de `unsigned short' a utiliser dans qsort.
327--
328*/
329int qSort_Ushort(const void *a1,const void *a2)
330{
331if( *((unsigned short *) a1) < *((unsigned short *) a2) ) return(-1);
332if( *((unsigned short *) a1) > *((unsigned short *) a2) ) return( 1);
333return(0);
334}
335
336/*===================================================================================*/
337/*
338++
339 int qSort_Short(const void *a1,const void *a2)
340 Fonction de tri de `short' a utiliser dans qsort.
341--
342*/
343int qSort_Short(const void *a1,const void *a2)
344{
345if( *((short *) a1) < *((short *) a2) ) return(-1);
346if( *((short *) a1) > *((short *) a2) ) return( 1);
347return(0);
348}
349
350/*===================================================================================*/
351/*
352++
353 void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
354 Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
355 such that arr[indx[j]] is in ascending order for j=1,2,,N.
356 The input quantities n and arr are not changed.
357--
358*/
359#define SWAP_INDEXR4(a,b) itemp=(a);(a)=(b);(b)=itemp;
360#define M 7
361#define NSTACK 50
362void IndexR4(int_4 n, float* arr_c, int_4* indx_c)
363/* encore du Num.Rec. avec tableaux commencant a 1. */
364{
365 float *arr; int_4 *indx;
366 int_4 i,indxt,ir=n,itemp,j,k,l=1;
367 int_4 jstack=0;
368 float a;
369 int_4 istack[NSTACK+1];
370
371 arr = arr_c-1;
372 indx = indx_c-1;
373
374 for (j=1;j<=n;j++) indx[j]=j;
375 for (;;) {
376 if (ir-l < M) {
377 for (j=l+1;j<=ir;j++) {
378 indxt=indx[j];
379 a=arr[indxt];
380 for (i=j-1;i>=1;i--) {
381 if (arr[indx[i]] <= a) break;
382 indx[i+1]=indx[i];
383 }
384 indx[i+1]=indxt;
385 }
386 if (jstack == 0) break;
387 ir=istack[jstack--];
388 l=istack[jstack--];
389 } else {
390 k=(l+ir) >> 1;
391 SWAP_INDEXR4(indx[k],indx[l+1]);
392 if (arr[indx[l+1]] > arr[indx[ir]]) {
393 SWAP_INDEXR4(indx[l+1],indx[ir])
394 }
395 if (arr[indx[l]] > arr[indx[ir]]) {
396 SWAP_INDEXR4(indx[l],indx[ir])
397 }
398 if (arr[indx[l+1]] > arr[indx[l]]) {
399 SWAP_INDEXR4(indx[l+1],indx[l])
400 }
401 i=l+1;
402 j=ir;
403 indxt=indx[l];
404 a=arr[indxt];
405 for (;;) {
406 do i++; while (arr[indx[i]] < a);
407 do j--; while (arr[indx[j]] > a);
408 if (j < i) break;
409 SWAP_INDEXR4(indx[i],indx[j])
410 }
411 indx[l]=indx[j];
412 indx[j]=indxt;
413 jstack += 2;
414 if (jstack > NSTACK) {
415 printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
416 exit(-1);
417 }
418 if (ir-i+1 >= j-l) {
419 istack[jstack]=ir;
420 istack[jstack-1]=i;
421 ir=j-1;
422 } else {
423 istack[jstack]=j-1;
424 istack[jstack-1]=l;
425 l=i;
426 }
427 }
428 }
[1474]429 /* pour avoir un retour d'indice entre 0 et n-1 */
430 for (j=0;j<n;j++) indx_c[j]=indx_c[j]-1;
[220]431}
432#undef SWAP_INDEXR4
433#undef M
434#undef NSTACK
435
436/*===================================================================================*/
437/*
438++
439 void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
440 Indexes an array arr[1..n], i.e., outputs the array indx[1..n]
441 such that arr[indx[j]] is in ascending order for j=1,2,,N.
442 The input quantities n and arr are not changed.
443--
444*/
445#define SWAP_INDEXR8(a,b) itemp=(a);(a)=(b);(b)=itemp;
446#define M 7
447#define NSTACK 50
448void IndexR8(int_4 n, double* arr_c, int_4* indx_c)
449/* encore du Num.Rec. avec tableaux commencant a 1. */
450{
451 double *arr; int_4 *indx;
452 int_4 i,indxt,ir=n,itemp,j,k,l=1;
453 int_4 jstack=0;
454 double a;
455 int_4 istack[NSTACK+1];
456
457 arr = arr_c-1;
458 indx = indx_c-1;
459
460 for (j=1;j<=n;j++) indx[j]=j;
461 for (;;) {
462 if (ir-l < M) {
463 for (j=l+1;j<=ir;j++) {
464 indxt=indx[j];
465 a=arr[indxt];
466 for (i=j-1;i>=1;i--) {
467 if (arr[indx[i]] <= a) break;
468 indx[i+1]=indx[i];
469 }
470 indx[i+1]=indxt;
471 }
472 if (jstack == 0) break;
473 ir=istack[jstack--];
474 l=istack[jstack--];
475 } else {
476 k=(l+ir) >> 1;
477 SWAP_INDEXR8(indx[k],indx[l+1]);
478 if (arr[indx[l+1]] > arr[indx[ir]]) {
479 SWAP_INDEXR8(indx[l+1],indx[ir])
480 }
481 if (arr[indx[l]] > arr[indx[ir]]) {
482 SWAP_INDEXR8(indx[l],indx[ir])
483 }
484 if (arr[indx[l+1]] > arr[indx[l]]) {
485 SWAP_INDEXR8(indx[l+1],indx[l])
486 }
487 i=l+1;
488 j=ir;
489 indxt=indx[l];
490 a=arr[indxt];
491 for (;;) {
492 do i++; while (arr[indx[i]] < a);
493 do j--; while (arr[indx[j]] > a);
494 if (j < i) break;
495 SWAP_INDEXR8(indx[i],indx[j])
496 }
497 indx[l]=indx[j];
498 indx[j]=indxt;
499 jstack += 2;
500 if (jstack > NSTACK) {
501 printf("NSTACK too small in indexx. %d>%d",jstack,NSTACK);
502 exit(-1);
503 }
504 if (ir-i+1 >= j-l) {
505 istack[jstack]=ir;
506 istack[jstack-1]=i;
507 ir=j-1;
508 } else {
509 istack[jstack]=j-1;
510 istack[jstack-1]=l;
511 l=i;
512 }
513 }
514 }
[1474]515 /* pour avoir un retour d'indice entre 0 et n-1 */
516 for (j=0;j<n;j++) indx_c[j]=indx_c[j]-1;
[220]517}
518#undef SWAP_INDEXR8
519#undef M
520#undef NSTACK
Note: See TracBrowser for help on using the repository browser.