Mercurial > hg > octave-nkf

comparison liboctave/Array-d.cc @ 7433:402168152bb9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2008-01-31 18:59:09 by dbateman]
author dbateman
date Thu, 31 Jan 2008 18:59:11 +0000
parents a1dbe9d80eee
children d98dea7d16b0
comparison
equal deleted inserted replaced
7432:3c999b2b5de8 7433:402168152bb9
27 27
28 // Instantiate Arrays of double values. 28 // Instantiate Arrays of double values.
29 29
30 #include "Array.h" 30 #include "Array.h"
31 #include "Array.cc" 31 #include "Array.cc"
32 #include "oct-sort.cc"
33
34 #if defined (HAVE_IEEE754_DATA_FORMAT)
35
36 static inline uint64_t
37 FloatFlip (uint64_t f)
38 {
39 uint64_t mask
40 = -static_cast<int64_t>(f >> 63) | 0x8000000000000000ULL;
41
42 return f ^ mask;
43 }
44
45 static inline uint64_t
46 IFloatFlip (uint64_t f)
47 {
48 uint64_t mask = ((f >> 63) - 1) | 0x8000000000000000ULL;
49
50 return f ^ mask;
51 }
52
53 template <>
54 bool
55 ascending_compare (double a, double b)
56 {
57 return (xisnan (b) || (a < b));
58 }
59
60 template <>
61 bool
62 ascending_compare (vec_index<double> *a, vec_index<double> *b)
63 {
64 return (xisnan (b->vec) || (a->vec < b->vec));
65 }
66
67 template <>
68 bool
69 descending_compare (double a, double b)
70 {
71 return (xisnan (a) || (a > b));
72 }
73
74 template <>
75 bool
76 descending_compare (vec_index<double> *a, vec_index<double> *b)
77 {
78 return (xisnan (b->vec) || (a->vec > b->vec));
79 }
80
81 INSTANTIATE_ARRAY_SORT (uint64_t);
82
83 template <>
84 Array<double>
85 Array<double>::sort (octave_idx_type dim, sortmode mode) const
86 {
87 Array<double> m = *this;
88
89 dim_vector dv = m.dims ();
90
91 if (m.length () < 1)
92 return m;
93
94 octave_idx_type ns = dv(dim);
95 octave_idx_type iter = dv.numel () / ns;
96 octave_idx_type stride = 1;
97 for (int i = 0; i < dim; i++)
98 stride *= dv(i);
99
100 double *v = m.fortran_vec ();
101
102 uint64_t *p = reinterpret_cast<uint64_t *> (v);
103
104 octave_sort<uint64_t> lsort;
105
106 if (mode == ASCENDING)
107 lsort.set_compare (ascending_compare);
108 else if (mode == DESCENDING)
109 lsort.set_compare (descending_compare);
110
111 if (stride == 1)
112 {
113 for (octave_idx_type j = 0; j < iter; j++)
114 {
115 // Flip the data in the vector so that int compares on
116 // IEEE754 give the correct ordering.
117
118 for (octave_idx_type i = 0; i < ns; i++)
119 p[i] = FloatFlip (p[i]);
120
121 lsort.sort (p, ns);
122
123 // Flip the data out of the vector so that int compares
124 // on IEEE754 give the correct ordering.
125
126 for (octave_idx_type i = 0; i < ns; i++)
127 p[i] = IFloatFlip (p[i]);
128
129 // There are two representations of NaN. One will be
130 // sorted to the beginning of the vector and the other
131 // to the end. If it will be sorted incorrectly, fix
132 // things up.
133
134 if (lo_ieee_signbit (octave_NaN))
135 {
136 if (mode == UNDEFINED || mode == ASCENDING)
137 {
138 octave_idx_type i = 0;
139 double *vtmp = reinterpret_cast<double *> (p);
140 while (xisnan (vtmp[i++]) && i < ns);
141 for (octave_idx_type l = 0; l < ns - i + 1; l++)
142 vtmp[l] = vtmp[l+i-1];
143 for (octave_idx_type l = ns - i + 1; l < ns; l++)
144 vtmp[l] = octave_NaN;
145 }
146 else
147 {
148 octave_idx_type i = ns;
149 double *vtmp = reinterpret_cast<double *> (p);
150 while (xisnan (vtmp[--i]) && i > 0);
151 for (octave_idx_type l = i; l >= 0; l--)
152 vtmp[l-i+ns-1] = vtmp[l];
153 for (octave_idx_type l = 0; l < ns - i - 1; l++)
154 vtmp[l] = octave_NaN;
155 }
156 }
157
158 p += ns;
159 }
160 }
161 else
162 {
163 OCTAVE_LOCAL_BUFFER (uint64_t, vi, ns);
164
165 for (octave_idx_type j = 0; j < iter; j++)
166 {
167 octave_idx_type offset = j;
168 octave_idx_type offset2 = 0;
169 while (offset >= stride)
170 {
171 offset -= stride;
172 offset2++;
173 }
174 offset += offset2 * stride * ns;
175
176 // Flip the data in the vector so that int compares on
177 // IEEE754 give the correct ordering.
178
179 for (octave_idx_type i = 0; i < ns; i++)
180 vi[i] = FloatFlip (p[i*stride + offset]);
181
182 lsort.sort (vi, ns);
183
184 // Flip the data out of the vector so that int compares
185 // on IEEE754 give the correct ordering.
186
187 for (octave_idx_type i = 0; i < ns; i++)
188 p[i*stride + offset] = IFloatFlip (vi[i]);
189
190 // There are two representations of NaN. One will be
191 // sorted to the beginning of the vector and the other
192 // to the end. If it will be sorted to the beginning,
193 // fix things up.
194
195 if (lo_ieee_signbit (octave_NaN))
196 {
197 if (mode == UNDEFINED || mode == ASCENDING)
198 {
199 octave_idx_type i = 0;
200 while (xisnan (v[i++*stride + offset]) && i < ns);
201 for (octave_idx_type l = 0; l < ns - i + 1; l++)
202 v[l*stride + offset] = v[(l+i-1)*stride + offset];
203 for (octave_idx_type l = ns - i + 1; l < ns; l++)
204 v[l*stride + offset] = octave_NaN;
205 }
206 else
207 {
208 octave_idx_type i = ns;
209 while (xisnan (v[--i*stride + offset]) && i > 0);
210 for (octave_idx_type l = i; l >= 0; l--)
211 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
212 for (octave_idx_type l = 0; l < ns - i - 1; l++)
213 v[l*stride + offset] = octave_NaN;
214 }
215 }
216 }
217 }
218
219 return m;
220 }
221
222 template <>
223 Array<double>
224 Array<double>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim,
225 sortmode mode) const
226 {
227 Array<double> m = *this;
228
229 dim_vector dv = m.dims ();
230
231 if (m.length () < 1)
232 {
233 sidx = Array<octave_idx_type> (dv);
234 return m;
235 }
236
237 octave_idx_type ns = dv(dim);
238 octave_idx_type iter = dv.numel () / ns;
239 octave_idx_type stride = 1;
240 for (int i = 0; i < dim; i++)
241 stride *= dv(i);
242
243 double *v = m.fortran_vec ();
244
245 uint64_t *p = reinterpret_cast<uint64_t *> (v);
246
247 octave_sort<vec_index<uint64_t> *> indexed_sort;
248
249 if (mode == ASCENDING)
250 indexed_sort.set_compare (ascending_compare);
251 else if (mode == DESCENDING)
252 indexed_sort.set_compare (descending_compare);
253
254 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t> *, vi, ns);
255 OCTAVE_LOCAL_BUFFER (vec_index<uint64_t>, vix, ns);
256
257 for (octave_idx_type i = 0; i < ns; i++)
258 vi[i] = &vix[i];
259
260 sidx = Array<octave_idx_type> (dv);
261
262 for (octave_idx_type j = 0; j < iter; j++)
263 {
264 octave_idx_type offset = j;
265 octave_idx_type offset2 = 0;
266 while (offset >= stride)
267 {
268 offset -= stride;
269 offset2++;
270 }
271 offset += offset2 * stride * ns;
272
273 // Flip the data in the vector so that int compares on
274 // IEEE754 give the correct ordering.
275
276 for (octave_idx_type i = 0; i < ns; i++)
277 {
278 vi[i]->vec = FloatFlip (p[i*stride + offset]);
279 vi[i]->indx = i;
280 }
281
282 indexed_sort.sort (vi, ns);
283
284 // Flip the data out of the vector so that int compares on
285 // IEEE754 give the correct ordering
286
287 for (octave_idx_type i = 0; i < ns; i++)
288 {
289 p[i*stride + offset] = IFloatFlip (vi[i]->vec);
290 sidx(i*stride + offset) = vi[i]->indx;
291 }
292
293 // There are two representations of NaN. One will be sorted
294 // to the beginning of the vector and the other to the end.
295 // If it will be sorted to the beginning, fix things up.
296
297 if (lo_ieee_signbit (octave_NaN))
298 {
299 if (mode == UNDEFINED || mode == ASCENDING)
300 {
301 octave_idx_type i = 0;
302 while (xisnan (v[i++*stride+offset]) && i < ns);
303 OCTAVE_LOCAL_BUFFER (double, itmp, i - 1);
304 for (octave_idx_type l = 0; l < i -1; l++)
305 itmp[l] = sidx(l*stride + offset);
306 for (octave_idx_type l = 0; l < ns - i + 1; l++)
307 {
308 v[l*stride + offset] = v[(l+i-1)*stride + offset];
309 sidx(l*stride + offset) = sidx((l+i-1)*stride + offset);
310 }
311 for (octave_idx_type k = 0, l = ns - i + 1; l < ns; l++, k++)
312 {
313 v[l*stride + offset] = octave_NaN;
314 sidx(l*stride + offset) =
315 static_cast<octave_idx_type>(itmp[k]);
316 }
317 }
318 else
319 {
320 octave_idx_type i = ns;
321 while (xisnan (v[--i*stride+offset]) && i > 0);
322 OCTAVE_LOCAL_BUFFER (double, itmp, ns - i - 1);
323 for (octave_idx_type l = 0; l < ns - i -1; l++)
324 itmp[l] = sidx((l+i+1)*stride + offset);
325 for (octave_idx_type l = i; l >= 0; l--)
326 {
327 v[(l-i+ns-1)*stride + offset] = v[l*stride + offset];
328 sidx((l-i+ns-1)*stride + offset) = sidx(l*stride + offset);
329 }
330 for (octave_idx_type k = 0, l = 0; l < ns - i - 1; l++, k++)
331 {
332 v[l*stride + offset] = octave_NaN;
333 sidx(l*stride + offset) =
334 static_cast<octave_idx_type>(itmp[k]);
335 }
336 }
337 }
338 }
339
340 return m;
341 }
342
343 #else
344
345 template <>
346 bool
347 Array<double>::ascending_compare (double a, double b) const
348 {
349 return (xisnan (b) || (a < b));
350 }
351
352 template <>
353 bool
354 Array<double>::ascending_compare (vec_index<double> *a,
355 vec_index<double> *b) const
356 {
357 return (xisnan (b->vec) || (a->vec < b->vec));
358 }
359
360 template <>
361 bool
362 Array<double>::descending_compare (double a, double b) const
363 {
364 return (xisnan (a) || (a > b));
365 }
366
367 template <>
368 bool
369 Array<double>::descending_compare (vec_index<double> *a,
370 vec_index<double> *b) const
371 {
372 return (xisnan (b->vec) || (a->vec > b->vec));
373 }
374
375 INSTANTIATE_ARRAY_SORT (double);
376
377 #endif
32 378
33 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API); 379 INSTANTIATE_ARRAY_AND_ASSIGN (double, OCTAVE_API);
34 380
35 INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API); 381 INSTANTIATE_ARRAY_ASSIGN (double, int, OCTAVE_API);
36 INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API); 382 INSTANTIATE_ARRAY_ASSIGN (double, short, OCTAVE_API);