2928
|
1 /* |
|
2 |
|
3 Copyright (C) 1996, 1997 John W. Eaton |
4850
|
4 Copyright (C) 2004 David Bateman |
2928
|
5 |
|
6 This file is part of Octave. |
|
7 |
|
8 Octave is free software; you can redistribute it and/or modify it |
|
9 under the terms of the GNU General Public License as published by the |
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
11 later version. |
|
12 |
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
16 for more details. |
|
17 |
|
18 You should have received a copy of the GNU General Public License |
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
21 |
|
22 */ |
|
23 |
|
24 #ifdef HAVE_CONFIG_H |
|
25 #include <config.h> |
|
26 #endif |
|
27 |
3826
|
28 #include "lo-mappers.h" |
4153
|
29 #include "quit.h" |
3826
|
30 |
2928
|
31 #include "defun-dld.h" |
|
32 #include "error.h" |
|
33 #include "gripes.h" |
|
34 #include "oct-obj.h" |
4850
|
35 #include "lo-ieee.h" |
|
36 #include "data-conv.h" |
|
37 #include "ov-cx-mat.h" |
|
38 #include "oct-sort.cc" |
2928
|
39 |
4850
|
40 /* If we are IEEE 754 or IEEE 854 compliant, then we can use the trick of |
|
41 * casting doubles as unsigned eight byte integers, and with a little |
|
42 * bit of magic we can automatically sort the NaN's correctly. |
|
43 */ |
2928
|
44 |
4850
|
45 #if defined(HAVE_IEEE754_COMPLIANCE) && defined(EIGHT_BYTE_INT) |
|
46 |
|
47 static inline unsigned EIGHT_BYTE_INT FloatFlip(unsigned EIGHT_BYTE_INT f) |
2928
|
48 { |
4850
|
49 unsigned EIGHT_BYTE_INT mask = -(EIGHT_BYTE_INT)(f >> 63) | |
|
50 0x8000000000000000ULL; |
|
51 return f ^ mask; |
|
52 } |
2928
|
53 |
4850
|
54 inline unsigned EIGHT_BYTE_INT IFloatFlip(unsigned EIGHT_BYTE_INT f) |
|
55 { |
|
56 unsigned EIGHT_BYTE_INT mask = ((f >> 63) - 1) | 0x8000000000000000ULL; |
|
57 return f ^ mask; |
|
58 } |
2928
|
59 |
4850
|
60 struct vec_index |
|
61 { |
|
62 unsigned EIGHT_BYTE_INT vec; |
|
63 int indx; |
|
64 }; |
2928
|
65 |
4850
|
66 bool |
|
67 ieee754_compare (vec_index *a, vec_index *b) |
|
68 { |
|
69 return (a->vec < b->vec); |
2928
|
70 } |
|
71 |
4850
|
72 template octave_sort<unsigned EIGHT_BYTE_INT>; |
|
73 template octave_sort<vec_index *>; |
|
74 #else |
|
75 struct vec_index |
|
76 { |
|
77 double vec; |
|
78 int indx; |
|
79 }; |
2928
|
80 |
4850
|
81 bool |
|
82 double_compare (double a, double b) |
|
83 { |
|
84 return (xisnan(b) || (a < b)); |
|
85 } |
2928
|
86 |
4850
|
87 bool |
|
88 double_compare (vec_index *a, vec_index *b) |
|
89 { |
|
90 return (xisnan(b->vec) || (a->vec < b->vec)); |
|
91 } |
|
92 |
|
93 template octave_sort<double>; |
|
94 template octave_sort<vec_index *>; |
|
95 #endif |
2928
|
96 |
4850
|
97 struct complex_vec_index |
|
98 { |
|
99 Complex vec; |
|
100 int indx; |
|
101 }; |
2928
|
102 |
4850
|
103 bool |
|
104 complex_compare (complex_vec_index *a, complex_vec_index *b) |
|
105 { |
|
106 return (xisnan(b->vec) || (abs(a->vec) < abs(b->vec))); |
|
107 } |
2928
|
108 |
4850
|
109 template octave_sort<complex_vec_index *>; |
2928
|
110 |
|
111 static octave_value_list |
4850
|
112 mx_sort (NDArray &m, bool return_idx, int dim) |
2928
|
113 { |
|
114 octave_value_list retval; |
|
115 |
4850
|
116 if (m.length () < 1) |
|
117 return retval; |
|
118 |
|
119 dim_vector dv = m.dims (); |
4851
|
120 unsigned int ns = dv (dim); |
|
121 unsigned int iter = dv.numel () / ns; |
|
122 unsigned int stride = 1; |
|
123 for (unsigned int i = 0; i < (unsigned int)dim; i++) |
4850
|
124 stride *= dv(i); |
|
125 |
|
126 #if defined(HAVE_IEEE754_COMPLIANCE) && defined(EIGHT_BYTE_INT) |
|
127 double *v = m.fortran_vec (); |
|
128 |
|
129 unsigned EIGHT_BYTE_INT *p = (unsigned EIGHT_BYTE_INT *)v; |
|
130 |
|
131 if (return_idx) |
|
132 { |
|
133 octave_sort<vec_index *> indexed_ieee754_sort (ieee754_compare); |
|
134 |
|
135 OCTAVE_LOCAL_BUFFER (vec_index *, vi, ns); |
|
136 OCTAVE_LOCAL_BUFFER (vec_index, vix, ns); |
|
137 |
4851
|
138 for (unsigned int i = 0; i < ns; i++) |
4850
|
139 vi[i] = &vix[i]; |
|
140 |
|
141 NDArray idx (dv); |
|
142 |
4851
|
143 for (unsigned int j = 0; j < iter; j++) |
4850
|
144 { |
4851
|
145 unsigned int offset = j; |
|
146 unsigned int offset2 = 0; |
4850
|
147 while (offset >= stride) |
|
148 { |
|
149 offset -= stride; |
|
150 offset2++; |
|
151 } |
|
152 offset += offset2 * stride * ns; |
2928
|
153 |
4850
|
154 /* Flip the data in the vector so that int compares on |
|
155 * IEEE754 give the correct ordering |
|
156 */ |
4851
|
157 for (unsigned int i = 0; i < ns; i++) |
4850
|
158 { |
|
159 vi[i]->vec = FloatFlip (p[i*stride + offset]); |
|
160 vi[i]->indx = i + 1; |
|
161 } |
|
162 |
|
163 indexed_ieee754_sort.sort (vi, ns); |
|
164 |
|
165 /* Flip the data out of the vector so that int compares on |
|
166 * IEEE754 give the correct ordering |
|
167 */ |
4851
|
168 for (unsigned int i = 0; i < ns; i++) |
4850
|
169 { |
|
170 p[i*stride + offset] = IFloatFlip (vi[i]->vec); |
|
171 idx(i*stride + offset) = vi[i]->indx; |
|
172 } |
2928
|
173 |
4850
|
174 /* There are two representations of NaN. One will be sorted to |
|
175 * the beginning of the vector and the other to the end. If it |
|
176 * will be sorted to the beginning, fix things up. |
|
177 */ |
|
178 if (lo_ieee_signbit (octave_NaN)) |
|
179 { |
|
180 unsigned int i = 0; |
|
181 while (xisnan(v[i++*stride+offset]) && i < ns); |
|
182 OCTAVE_LOCAL_BUFFER (double, itmp, i - 1); |
|
183 for (unsigned int l = 0; l < i -1; l++) |
|
184 itmp[l] = idx(l*stride + offset); |
|
185 for (unsigned int l = 0; l < ns - i + 1; l++) |
|
186 { |
|
187 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
188 idx(l*stride + offset) = idx((l+i-1)*stride + offset); |
|
189 } |
|
190 for (unsigned int k = 0, l = ns - i + 1; l < ns; l++, k++) |
|
191 { |
|
192 v[l*stride + offset] = octave_NaN; |
|
193 idx(l*stride + offset) = itmp[k]; |
|
194 } |
|
195 } |
|
196 } |
2928
|
197 |
4850
|
198 retval (1) = idx; |
2928
|
199 } |
4850
|
200 else |
2928
|
201 { |
4850
|
202 octave_sort<unsigned EIGHT_BYTE_INT> ieee754_sort; |
|
203 |
|
204 if (stride == 1) |
2928
|
205 { |
4851
|
206 for (unsigned int j = 0; j < iter; j++) |
4850
|
207 { |
|
208 /* Flip the data in the vector so that int compares on |
|
209 * IEEE754 give the correct ordering |
|
210 */ |
4851
|
211 for (unsigned int i = 0; i < ns; i++) |
4850
|
212 p[i] = FloatFlip (p[i]); |
|
213 |
|
214 ieee754_sort.sort (p, ns); |
|
215 |
|
216 /* Flip the data out of the vector so that int compares on |
|
217 * IEEE754 give the correct ordering |
|
218 */ |
4851
|
219 for (unsigned int i = 0; i < ns; i++) |
4850
|
220 p[i] = IFloatFlip (p[i]); |
|
221 |
|
222 /* There are two representations of NaN. One will be sorted to |
|
223 * the beginning of the vector and the other to the end. If it |
|
224 * will be sorted to the beginning, fix things up. |
|
225 */ |
|
226 if (lo_ieee_signbit (octave_NaN)) |
|
227 { |
|
228 unsigned int i = 0; |
|
229 double *vtmp = (double *)p; |
|
230 while (xisnan(vtmp[i++]) && i < ns); |
|
231 for (unsigned int l = 0; l < ns - i + 1; l++) |
|
232 vtmp[l] = vtmp[l+i-1]; |
|
233 for (unsigned int l = ns - i + 1; l < ns; l++) |
|
234 vtmp[l] = octave_NaN; |
|
235 } |
|
236 |
|
237 p += ns; |
|
238 } |
2928
|
239 |
4850
|
240 } |
|
241 else |
|
242 { |
4851
|
243 OCTAVE_LOCAL_BUFFER (unsigned EIGHT_BYTE_INT, vi, ns); |
|
244 |
|
245 for (unsigned int j = 0; j < iter; j++) |
4850
|
246 { |
4851
|
247 unsigned int offset = j; |
|
248 unsigned int offset2 = 0; |
4850
|
249 while (offset >= stride) |
|
250 { |
|
251 offset -= stride; |
|
252 offset2++; |
|
253 } |
|
254 offset += offset2 * stride * ns; |
|
255 |
|
256 /* Flip the data in the vector so that int compares on |
|
257 * IEEE754 give the correct ordering |
|
258 */ |
4851
|
259 for (unsigned int i = 0; i < ns; i++) |
4850
|
260 vi[i] = FloatFlip (p[i*stride + offset]); |
2928
|
261 |
4850
|
262 ieee754_sort.sort (vi, ns); |
|
263 |
|
264 /* Flip the data out of the vector so that int compares on |
|
265 * IEEE754 give the correct ordering |
|
266 */ |
4851
|
267 for (unsigned int i = 0; i < ns; i++) |
4850
|
268 p[i*stride + offset] = IFloatFlip (vi[i]); |
|
269 |
|
270 /* There are two representations of NaN. One will be sorted to |
|
271 * the beginning of the vector and the other to the end. If it |
|
272 * will be sorted to the beginning, fix things up. |
|
273 */ |
|
274 if (lo_ieee_signbit (octave_NaN)) |
|
275 { |
|
276 unsigned int i = 0; |
|
277 while (xisnan(v[i++*stride + offset]) && i < ns); |
|
278 for (unsigned int l = 0; l < ns - i + 1; l++) |
|
279 v[l*stride + offset] = v[(l+i-1)*stride + offset]; |
|
280 for (unsigned int l = ns - i + 1; l < ns; l++) |
|
281 v[l*stride + offset] = octave_NaN; |
|
282 } |
|
283 } |
2928
|
284 } |
|
285 } |
4850
|
286 #else |
|
287 if (return_idx) |
|
288 { |
|
289 double *v = m.fortran_vec (); |
|
290 octave_sort<vec_index *> indexed_double_sort (double_compare); |
2928
|
291 |
4850
|
292 OCTAVE_LOCAL_BUFFER (vec_index *, vi, ns); |
|
293 OCTAVE_LOCAL_BUFFER (vec_index, vix, ns); |
|
294 |
4851
|
295 for (unsigned int i = 0; i < ns; i++) |
4850
|
296 vi[i] = &vix[i]; |
|
297 |
|
298 NDArray idx (dv); |
|
299 |
|
300 if (stride == 1) |
|
301 { |
4851
|
302 for (unsigned int j = 0; j < iter; j++) |
4850
|
303 { |
4851
|
304 unsigned int offset = j * ns; |
4850
|
305 |
4851
|
306 for (unsigned int i = 0; i < ns; i++) |
4850
|
307 { |
|
308 vi[i]->vec = v[i]; |
|
309 vi[i]->indx = i + 1; |
|
310 } |
|
311 |
|
312 indexed_double_sort.sort (vi, ns); |
|
313 |
4851
|
314 for (unsigned int i = 0; i < ns; i++) |
4850
|
315 { |
|
316 v[i] = vi[i]->vec; |
|
317 idx(i + offset) = vi[i]->indx; |
|
318 } |
|
319 v += ns; |
|
320 } |
|
321 } |
|
322 else |
|
323 { |
4851
|
324 for (unsigned int j = 0; j < iter; j++) |
4850
|
325 { |
4851
|
326 unsigned int offset = j; |
|
327 unsigned int offset2 = 0; |
4850
|
328 while (offset >= stride) |
|
329 { |
|
330 offset -= stride; |
|
331 offset2++; |
|
332 } |
|
333 offset += offset2 * stride * ns; |
|
334 |
4851
|
335 for (unsigned int i = 0; i < ns; i++) |
4850
|
336 { |
|
337 vi[i]->vec = v[i*stride + offset]; |
|
338 vi[i]->indx = i + 1; |
|
339 } |
2928
|
340 |
4850
|
341 indexed_double_sort.sort (vi, ns); |
|
342 |
4851
|
343 for (unsigned int i = 0; i < ns; i++) |
4850
|
344 { |
|
345 v[i*stride+offset] = vi[i]->vec; |
|
346 idx(i*stride+offset) = vi[i]->indx; |
|
347 } |
|
348 } |
|
349 } |
|
350 retval (1) = idx; |
|
351 } |
|
352 else |
|
353 { |
|
354 double *v = m.fortran_vec (); |
|
355 octave_sort<double> double_sort (double_compare); |
|
356 |
|
357 if (stride == 1) |
4851
|
358 for (unsigned int j = 0; j < iter; j++) |
4850
|
359 { |
|
360 double_sort.sort (v, ns); |
|
361 v += ns; |
|
362 } |
|
363 else |
|
364 { |
|
365 OCTAVE_LOCAL_BUFFER (double, vi, ns); |
4851
|
366 for (unsigned int j = 0; j < iter; j++) |
4850
|
367 { |
4851
|
368 unsigned int offset = j; |
|
369 unsigned int offset2 = 0; |
4850
|
370 while (offset >= stride) |
|
371 { |
|
372 offset -= stride; |
|
373 offset2++; |
|
374 } |
|
375 offset += offset2 * stride * ns; |
|
376 |
4851
|
377 for (unsigned int i = 0; i < ns; i++) |
4850
|
378 vi[i] = v[i*stride + offset]; |
|
379 |
|
380 double_sort.sort (vi, ns); |
|
381 |
4851
|
382 for (unsigned int i = 0; i < ns; i++) |
4850
|
383 v[i*stride + offset] = vi[i]; |
|
384 } |
|
385 } |
|
386 } |
|
387 #endif |
|
388 retval(0) = m; |
2928
|
389 return retval; |
|
390 } |
|
391 |
|
392 static octave_value_list |
4850
|
393 mx_sort (ComplexNDArray &m, bool return_idx, int dim) |
2928
|
394 { |
|
395 octave_value_list retval; |
|
396 |
4850
|
397 if (m.length () < 1) |
|
398 return retval; |
|
399 |
|
400 dim_vector dv = m.dims (); |
4851
|
401 unsigned int ns = dv (dim); |
|
402 unsigned int iter = dv.numel () / ns; |
|
403 unsigned int stride = 1; |
|
404 for (unsigned int i = 0; i < (unsigned int)dim; i++) |
4850
|
405 stride *= dv(i); |
2928
|
406 |
4850
|
407 octave_sort<complex_vec_index *> indexed_double_sort (complex_compare); |
|
408 |
|
409 Complex *v = m.fortran_vec (); |
2928
|
410 |
4850
|
411 OCTAVE_LOCAL_BUFFER (complex_vec_index *, vi, ns); |
|
412 OCTAVE_LOCAL_BUFFER (complex_vec_index, vix, ns); |
|
413 |
4851
|
414 for (unsigned int i = 0; i < ns; i++) |
4850
|
415 vi[i] = &vix[i]; |
2928
|
416 |
4850
|
417 NDArray idx (dv); |
|
418 |
|
419 if (stride == 1) |
2928
|
420 { |
4851
|
421 for (unsigned int j = 0; j < iter; j++) |
2928
|
422 { |
4851
|
423 unsigned int offset = j * ns; |
2928
|
424 |
4851
|
425 for (unsigned int i = 0; i < ns; i++) |
4153
|
426 { |
4850
|
427 vi[i]->vec = v[i]; |
|
428 vi[i]->indx = i + 1; |
|
429 } |
|
430 |
|
431 indexed_double_sort.sort (vi, ns); |
|
432 |
|
433 if (return_idx) |
|
434 { |
4851
|
435 for (unsigned int i = 0; i < ns; i++) |
4153
|
436 { |
4850
|
437 v[i] = vi[i]->vec; |
|
438 idx(i + offset) = vi[i]->indx; |
4153
|
439 } |
|
440 } |
4850
|
441 else |
|
442 { |
4851
|
443 for (unsigned int i = 0; i < ns; i++) |
4850
|
444 v[i] = vi[i]->vec; |
|
445 } |
|
446 v += ns; |
|
447 } |
|
448 } |
|
449 else |
|
450 { |
4851
|
451 for (unsigned int j = 0; j < iter; j++) |
4850
|
452 { |
4851
|
453 unsigned int offset = j; |
|
454 unsigned int offset2 = 0; |
4850
|
455 while (offset >= stride) |
|
456 { |
|
457 offset -= stride; |
|
458 offset2++; |
|
459 } |
|
460 offset += offset2 * stride * ns; |
2928
|
461 |
4851
|
462 for (unsigned int i = 0; i < ns; i++) |
4850
|
463 { |
|
464 vi[i]->vec = v[i*stride + offset]; |
|
465 vi[i]->indx = i + 1; |
|
466 } |
|
467 |
|
468 indexed_double_sort.sort (vi, ns); |
|
469 |
|
470 if (return_idx) |
|
471 { |
4851
|
472 for (unsigned int i = 0; i < ns; i++) |
4850
|
473 { |
|
474 v[i*stride + offset] = vi[i]->vec; |
|
475 idx(i*stride + offset) = vi[i]->indx; |
|
476 } |
|
477 } |
|
478 else |
|
479 { |
4851
|
480 for (unsigned int i = 0; i < ns; i++) |
4850
|
481 v[i*stride + offset] = vi[i]->vec; |
|
482 } |
2928
|
483 } |
|
484 } |
|
485 |
4850
|
486 if (return_idx) |
|
487 retval (1) = idx; |
|
488 |
|
489 retval(0) = m; |
2928
|
490 |
|
491 return retval; |
|
492 } |
|
493 |
|
494 DEFUN_DLD (sort, args, nargout, |
3369
|
495 "-*- texinfo -*-\n\ |
|
496 @deftypefn {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x})\n\ |
4850
|
497 @deftypefnx {Loadable Function} {[@var{s}, @var{i}] =} sort (@var{x}, @var{dim})\n\ |
3369
|
498 Return a copy of @var{x} with the elements elements arranged in\n\ |
|
499 increasing order. For matrices, @code{sort} orders the elements in each\n\ |
|
500 column.\n\ |
|
501 \n\ |
|
502 For example,\n\ |
|
503 \n\ |
|
504 @example\n\ |
|
505 @group\n\ |
|
506 sort ([1, 2; 2, 3; 3, 1])\n\ |
|
507 @result{} 1 1\n\ |
|
508 2 2\n\ |
|
509 3 3\n\ |
|
510 @end group\n\ |
|
511 @end example\n\ |
2928
|
512 \n\ |
3369
|
513 The @code{sort} function may also be used to produce a matrix\n\ |
|
514 containing the original row indices of the elements in the sorted\n\ |
|
515 matrix. For example,\n\ |
|
516 \n\ |
|
517 @example\n\ |
|
518 @group\n\ |
|
519 [s, i] = sort ([1, 2; 2, 3; 3, 1])\n\ |
|
520 @result{} s = 1 1\n\ |
|
521 2 2\n\ |
|
522 3 3\n\ |
|
523 @result{} i = 1 3\n\ |
|
524 2 1\n\ |
|
525 3 2\n\ |
|
526 @end group\n\ |
|
527 @end example\n\ |
4850
|
528 \n\ |
|
529 If the optional argument @var{dim} is given, then the matrix is sorted\n\ |
|
530 along the dimension defined by @var{dim}.\n\ |
|
531 \n\ |
|
532 For equal elements, the indices are such that the equal elements are listed\n\ |
|
533 in the order that appeared in the original list.\n\ |
|
534 \n\ |
|
535 The algorithm used in @code{sort} is optimized for the sorting of partially\n\ |
|
536 ordered lists.\n\ |
3369
|
537 @end deftypefn") |
2928
|
538 { |
|
539 octave_value_list retval; |
|
540 |
|
541 int nargin = args.length (); |
|
542 |
4850
|
543 if (nargin != 1 && nargin != 2) |
2928
|
544 { |
|
545 print_usage ("sort"); |
|
546 return retval; |
|
547 } |
|
548 |
4850
|
549 bool return_idx = nargout > 1; |
2928
|
550 |
|
551 octave_value arg = args(0); |
|
552 |
4850
|
553 int dim = 0; |
|
554 if (nargin == 2) |
|
555 dim = args(1).nint_value () - 1; |
|
556 |
|
557 dim_vector dv = ((const octave_complex_matrix&) arg) .dims (); |
|
558 if (error_state) |
|
559 { |
|
560 gripe_wrong_type_arg ("sort", arg); |
|
561 return retval; |
|
562 } |
|
563 if (nargin != 2) |
|
564 { |
|
565 // Find first non singleton dimension |
|
566 for (int i = 0; i < dv.length (); i++) |
|
567 if (dv(i) > 1) |
|
568 { |
|
569 dim = i; |
|
570 break; |
|
571 } |
|
572 } |
|
573 else |
|
574 { |
|
575 if (dim < 0 || dim > dv.length () - 1) |
|
576 { |
|
577 error ("sort: dim must be a valid dimension"); |
|
578 return retval; |
|
579 } |
|
580 } |
|
581 |
2928
|
582 if (arg.is_real_type ()) |
|
583 { |
4850
|
584 NDArray m = arg.array_value (); |
2928
|
585 |
|
586 if (! error_state) |
4850
|
587 retval = mx_sort (m, return_idx, dim); |
2928
|
588 } |
|
589 else if (arg.is_complex_type ()) |
|
590 { |
4850
|
591 ComplexNDArray cm = arg.complex_array_value (); |
2928
|
592 |
|
593 if (! error_state) |
4850
|
594 retval = mx_sort (cm, return_idx, dim); |
2928
|
595 } |
|
596 else |
|
597 gripe_wrong_type_arg ("sort", arg); |
|
598 |
|
599 return retval; |
|
600 } |
|
601 |
4850
|
602 |
2928
|
603 /* |
|
604 ;;; Local Variables: *** |
|
605 ;;; mode: C++ *** |
|
606 ;;; End: *** |
|
607 */ |