source: Sophya/trunk/SophyaLib/NTools/nbtrixx.cc@ 3678

Last change on this file since 3678 was 3678, checked in by cmv, 16 years ago

ajout routines template C++ pout trier, cmv 16/11/2009
File size: 7.3 KB
Line 
1#include "machdefs.h"
2#include <stdio.h>
3#include <stdlib.h>
4
5#include "nbtrixx.h"
6
7namespace SOPHYA {
8
9// attention, Numerical recipes prend des tableaux de 1 a n on remet
10// de 0 a n-1 en decramentant le pointeur du tableau d'entree
11
12//-------------------------------------------------------------
13#define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
14#define M 7
15#define NSTACK 50
16template <class T>
17void TabSort(uint_4 n, T* arr_c)
18{
19 uint_4 i,ir=n,j,k,l=1,*istack;
20 int jstack=0;
21 T a,temp;
22 T *arr = arr_c - 1; // MODIF
23
24 //istack=lvector(1,NSTACK);
25 istack = new uint_4[NSTACK+1]; // MODIF
26 for (;;) {
27 if (ir-l < M) {
28 for (j=l+1;j<=ir;j++) {
29 a=arr[j];
30 for (i=j-1;i>=l;i--) {
31 if (arr[i] <= a) break;
32 arr[i+1]=arr[i];
33 }
34 arr[i+1]=a;
35 }
36 if (jstack == 0) break;
37 ir=istack[jstack--];
38 l=istack[jstack--];
39 } else {
40 k=(l+ir) >> 1;
41 SWAP_NBTRI(arr[k],arr[l+1])
42 if (arr[l] > arr[ir]) {
43 SWAP_NBTRI(arr[l],arr[ir])
44 }
45 if (arr[l+1] > arr[ir]) {
46 SWAP_NBTRI(arr[l+1],arr[ir])
47 }
48 if (arr[l] > arr[l+1]) {
49 SWAP_NBTRI(arr[l],arr[l+1])
50 }
51 i=l+1;
52 j=ir;
53 a=arr[l+1];
54 for (;;) {
55 do i++; while (arr[i] < a);
56 do j--; while (arr[j] > a);
57 if (j < i) break;
58 SWAP_NBTRI(arr[i],arr[j]);
59 }
60 arr[l+1]=arr[j];
61 arr[j]=a;
62 jstack += 2;
63 if (jstack > NSTACK) throw("NSTACK too small in sort.");
64 if (ir-i+1 >= j-l) {
65 istack[jstack]=ir;
66 istack[jstack-1]=i;
67 ir=j-1;
68 } else {
69 istack[jstack]=j-1;
70 istack[jstack-1]=l;
71 l=i;
72 }
73 }
74 }
75 //free_lvector(istack,1,NSTACK);
76 delete [] istack; // MODIF
77}
78#undef M
79#undef NSTACK
80#undef SWAP_NBTRI
81
82//-------------------------------------------------------------
83#define SWAP_NBTRI(a,b) temp=(a);(a)=(b);(b)=temp;
84#define M 7
85#define NSTACK 50
86template <class T>
87void TabSort2(uint_4 n, T *arr_c, T *brr_c)
88{
89 uint_4 i,ir=n,j,k,l=1,*istack;
90 int jstack=0;
91 T a,b,temp;
92 T *arr = arr_c - 1; // MODIF
93 T *brr = brr_c - 1; // MODIF
94
95 //istack=lvector(1,NSTACK);
96 istack = new uint_4[NSTACK+1]; // MODIF
97 for (;;) {
98 if (ir-l < M) {
99 for (j=l+1;j<=ir;j++) {
100 a=arr[j];
101 b=brr[j];
102 for (i=j-1;i>=l;i--) {
103 if (arr[i] <= a) break;
104 arr[i+1]=arr[i];
105 brr[i+1]=brr[i];
106 }
107 arr[i+1]=a;
108 brr[i+1]=b;
109 }
110 if (!jstack) {
111 //free_lvector(istack,1,NSTACK);
112 delete [] istack; // MODIF
113 return;
114 }
115 ir=istack[jstack];
116 l=istack[jstack-1];
117 jstack -= 2;
118 } else {
119 k=(l+ir) >> 1;
120 SWAP_NBTRI(arr[k],arr[l+1])
121 SWAP_NBTRI(brr[k],brr[l+1])
122 if (arr[l] > arr[ir]) {
123 SWAP_NBTRI(arr[l],arr[ir])
124 SWAP_NBTRI(brr[l],brr[ir])
125 }
126 if (arr[l+1] > arr[ir]) {
127 SWAP_NBTRI(arr[l+1],arr[ir])
128 SWAP_NBTRI(brr[l+1],brr[ir])
129 }
130 if (arr[l] > arr[l+1]) {
131 SWAP_NBTRI(arr[l],arr[l+1])
132 SWAP_NBTRI(brr[l],brr[l+1])
133 }
134 i=l+1;
135 j=ir;
136 a=arr[l+1];
137 b=brr[l+1];
138 for (;;) {
139 do i++; while (arr[i] < a);
140 do j--; while (arr[j] > a);
141 if (j < i) break;
142 SWAP_NBTRI(arr[i],arr[j])
143 SWAP_NBTRI(brr[i],brr[j])
144 }
145 arr[l+1]=arr[j];
146 arr[j]=a;
147 brr[l+1]=brr[j];
148 brr[j]=b;
149 jstack += 2;
150 if (jstack > NSTACK) throw("NSTACK too small in sort2.");
151 if (ir-i+1 >= j-l) {
152 istack[jstack]=ir;
153 istack[jstack-1]=i;
154 ir=j-1;
155 } else {
156 istack[jstack]=j-1;
157 istack[jstack-1]=l;
158 l=i;
159 }
160 }
161 }
162}
163#undef M
164#undef NSTACK
165#undef SWAP_NBTRI
166
167//-------------------------------------------------------------
168#define SWAP_NBTRI(a,b) itemp=(a);(a)=(b);(b)=itemp;
169#define M 7
170#define NSTACK 50
171template <class T>
172void TabSortInd(uint_4 n, T *arr_c, uint_4 *indx_c)
173{
174 uint_4 i,indxt,ir=n,itemp,j,k,l=1;
175 int jstack=0,*istack;
176 T a;
177 T *arr = arr_c - 1; // MODIF
178 uint_4 *indx = indx_c - 1; // MODIF
179
180 //istack=ivector(1,NSTACK);
181 istack = new int[NSTACK+1]; // MODIF
182 for (j=1;j<=n;j++) indx[j]=j;
183 for (;;) {
184 if (ir-l < M) {
185 for (j=l+1;j<=ir;j++) {
186 indxt=indx[j];
187 a=arr[indxt];
188 for (i=j-1;i>=l;i--) {
189 if (arr[indx[i]] <= a) break;
190 indx[i+1]=indx[i];
191 }
192 indx[i+1]=indxt;
193 }
194 if (jstack == 0) break;
195 ir=istack[jstack--];
196 l=istack[jstack--];
197 } else {
198 k=(l+ir) >> 1;
199 SWAP_NBTRI(indx[k],indx[l+1]);
200 if (arr[indx[l]] > arr[indx[ir]]) {
201 SWAP_NBTRI(indx[l],indx[ir])
202 }
203 if (arr[indx[l+1]] > arr[indx[ir]]) {
204 SWAP_NBTRI(indx[l+1],indx[ir])
205 }
206 if (arr[indx[l]] > arr[indx[l+1]]) {
207 SWAP_NBTRI(indx[l],indx[l+1])
208 }
209 i=l+1;
210 j=ir;
211 indxt=indx[l+1];
212 a=arr[indxt];
213 for (;;) {
214 do i++; while (arr[indx[i]] < a);
215 do j--; while (arr[indx[j]] > a);
216 if (j < i) break;
217 SWAP_NBTRI(indx[i],indx[j])
218 }
219 indx[l+1]=indx[j];
220 indx[j]=indxt;
221 jstack += 2;
222 if (jstack > NSTACK) throw("NSTACK too small in indexx.");
223 if (ir-i+1 >= j-l) {
224 istack[jstack]=ir;
225 istack[jstack-1]=i;
226 ir=j-1;
227 } else {
228 istack[jstack]=j-1;
229 istack[jstack-1]=l;
230 l=i;
231 }
232 }
233 }
234 //free_ivector(istack,1,NSTACK);
235 delete [] istack; // MODIF
236 /* pour avoir un retour d'indice entre 0 et n-1 */
237 for(j=0;j<n;j++) indx_c[j]--; // MODIF
238}
239#undef M
240#undef NSTACK
241#undef SWAP_NBTRI
242
243//-------------------------------------------------------------
244#ifdef __CXX_PRAGMA_TEMPLATES__
245#pragma define_template TabSort<int_2>
246#pragma define_template TabSort<uint_2>
247#pragma define_template TabSort<int_4>
248#pragma define_template TabSort<uint_4>
249#pragma define_template TabSort<int_8>
250#pragma define_template TabSort<uint_8>
251#pragma define_template TabSort<r_4>
252#pragma define_template TabSort<r_8>
253
254#pragma define_template TabSort2<int_2>
255#pragma define_template TabSort2<uint_2>
256#pragma define_template TabSort2<int_4>
257#pragma define_template TabSort2<uint_4>
258#pragma define_template TabSort2<int_8>
259#pragma define_template TabSort2<uint_8>
260#pragma define_template TabSort2<r_4>
261#pragma define_template TabSort2<r_8>
262
263#pragma define_template TabSortInd<int_2>
264#pragma define_template TabSortInd<uint_2>
265#pragma define_template TabSortInd<int_4>
266#pragma define_template TabSortInd<uint_4>
267#pragma define_template TabSortInd<int_8>
268#pragma define_template TabSortInd<uint_8>
269#pragma define_template TabSortInd<r_4>
270#pragma define_template TabSortInd<r_8>
271#endif
272
273#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
274template void TabSort(uint_4, int_2*);
275template void TabSort(uint_4, uint_2*);
276template void TabSort(uint_4, int_4*);
277template void TabSort(uint_4, uint_4*);
278template void TabSort(uint_4, int_8*);
279template void TabSort(uint_4, uint_8*);
280template void TabSort(uint_4, r_4*);
281template void TabSort(uint_4, r_8*);
282
283template void TabSort2(uint_4, int_2*, int_2*);
284template void TabSort2(uint_4, uint_2*, uint_2*);
285template void TabSort2(uint_4, int_4*, int_4*);
286template void TabSort2(uint_4, uint_4*, uint_4*);
287template void TabSort2(uint_4, int_8*, int_8*);
288template void TabSort2(uint_4, uint_8*, uint_8*);
289template void TabSort2(uint_4, r_4*, r_4*);
290template void TabSort2(uint_4, r_8*, r_8*);
291
292template void TabSortInd(uint_4, int_2*, uint_4*);
293template void TabSortInd(uint_4, uint_2*, uint_4*);
294template void TabSortInd(uint_4, int_4*, uint_4*);
295template void TabSortInd(uint_4, uint_4*, uint_4*);
296template void TabSortInd(uint_4, int_8*, uint_4*);
297template void TabSortInd(uint_4, uint_8*, uint_4*);
298template void TabSortInd(uint_4, r_4*, uint_4*);
299template void TabSortInd(uint_4, r_8*, uint_4*);
300#endif
301
302} // FIN namespace SOPHYA
Note: See TracBrowser for help on using the repository browser.