1 /*********************************************************************
2 qsort -- Functions used by qsort to sort an array.
3 This is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Mohammad Akhlaghi <mohammad@akhlaghi.org>
7 Contributing author(s):
8 Copyright (C) 2015-2021, Free Software Foundation, Inc.
9
10 Gnuastro is free software: you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by the
12 Free Software Foundation, either version 3 of the License, or (at your
13 option) any later version.
14
15 Gnuastro is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
22 **********************************************************************/
23 #include <config.h>
24
25 #include <math.h>
26 #include <stdlib.h>
27 #include <stdint.h>
28
29 #include <fitsio.h>
30
31 #include <gnuastro/qsort.h>
32
33
34 /*****************************************************************/
35 /********** Macros ****************/
36 /*****************************************************************/
37 /* When one or both elements are NaN, the simple comparison, like '(tb >
38 ta) - (tb < ta)', will give 0 (as if the elements are equal). However,
39 some preference has to be given to the NaN element in a comparison,
40 otherwise the output is not going to be reasonable. We also don't want
41 to check NaNs on every comparison (it will slow down the processing).
42
43 So we'll exploit the fact that when there comparison result doesn't
44 equal zero, we don't have any NaNs and this 'COMPARE_FLOAT_POSTPROCESS'
45 macro is called only when the comparison gives zero. Being larger or
46 smaller isn't defined for NaNs, so we'll just put them in the end of the
47 sorted list whether it is sorted by decreasing or increasing mode.*/
48 #define COMPARE_FLOAT_POSTPROCESS ( \
49 isnan(ta) && isnan(tb) \
50 ? 0 /* Both NaN, define as equal. */ \
51 /* One is NaN, one isn't. */ \
52 : ( isnan(ta) \
53 ? 1 /* First is NaN, set as smaller. */ \
54 : ( isnan(tb) \
55 ? -1 /* Second is NaN, set as larger. */ \
56 : 0 ) /* None are NaN, set as equal.*/ \
57 ) \
58 )
59
60
61
62 /*****************************************************************/
63 /********** Sorting of actual array ****************/
64 /*****************************************************************/
65 int
gal_qsort_uint8_d(const void * a,const void * b)66 gal_qsort_uint8_d(const void *a, const void *b)
67 {
68 return ( *(uint8_t *)b - *(uint8_t *)a );
69 }
70
71 int
gal_qsort_uint8_i(const void * a,const void * b)72 gal_qsort_uint8_i(const void *a, const void *b)
73 {
74 return ( *(uint8_t *)a - *(uint8_t *)b );
75 }
76
77 int
gal_qsort_int8_d(const void * a,const void * b)78 gal_qsort_int8_d(const void *a, const void *b)
79 {
80 return ( *(int8_t *)b - *(int8_t *)a );
81 }
82
83 int
gal_qsort_int8_i(const void * a,const void * b)84 gal_qsort_int8_i(const void *a, const void *b)
85 {
86 return ( *(int8_t *)a - *(int8_t *)b );
87 }
88
89 int
gal_qsort_uint16_d(const void * a,const void * b)90 gal_qsort_uint16_d(const void *a, const void *b)
91 {
92 return ( *(uint16_t *)b - *(uint16_t *)a );
93 }
94
95 int
gal_qsort_uint16_i(const void * a,const void * b)96 gal_qsort_uint16_i(const void *a, const void *b)
97 {
98 return ( *(uint16_t *)a - *(uint16_t *)b );
99 }
100
101 int
gal_qsort_int16_d(const void * a,const void * b)102 gal_qsort_int16_d(const void *a, const void *b)
103 {
104 return ( *(int16_t *)b - *(int16_t *)a );
105 }
106
107 int
gal_qsort_int16_i(const void * a,const void * b)108 gal_qsort_int16_i(const void *a, const void *b)
109 {
110 return ( *(int16_t *)a - *(int16_t *)b );
111 }
112
113 int
gal_qsort_uint32_d(const void * a,const void * b)114 gal_qsort_uint32_d(const void *a, const void *b)
115 {
116 return ( *(uint32_t *)b - *(uint32_t *)a );
117 }
118
119 int
gal_qsort_uint32_i(const void * a,const void * b)120 gal_qsort_uint32_i(const void *a, const void *b)
121 {
122 return ( *(uint32_t *)a - *(uint32_t *)b );
123 }
124
125 int
gal_qsort_int32_d(const void * a,const void * b)126 gal_qsort_int32_d(const void *a, const void *b)
127 {
128 return ( *(int32_t *)b - *(int32_t *)a );
129 }
130
131 int
gal_qsort_int32_i(const void * a,const void * b)132 gal_qsort_int32_i(const void *a, const void *b)
133 {
134 return ( *(int32_t *)a - *(int32_t *)b );
135 }
136
137 int
gal_qsort_uint64_d(const void * a,const void * b)138 gal_qsort_uint64_d(const void *a, const void *b)
139 {
140 return ( *(uint64_t *)b - *(uint64_t *)a );
141 }
142
143 int
gal_qsort_uint64_i(const void * a,const void * b)144 gal_qsort_uint64_i(const void *a, const void *b)
145 {
146 return ( *(uint64_t *)a - *(uint64_t *)b );
147 }
148
149 int
gal_qsort_int64_d(const void * a,const void * b)150 gal_qsort_int64_d(const void *a, const void *b)
151 {
152 return ( *(int64_t *)b - *(int64_t *)a );
153 }
154
155 int
gal_qsort_int64_i(const void * a,const void * b)156 gal_qsort_int64_i(const void *a, const void *b)
157 {
158 return ( *(int64_t *)a - *(int64_t *)b );
159 }
160
161 int
gal_qsort_float32_d(const void * a,const void * b)162 gal_qsort_float32_d(const void *a, const void *b)
163 {
164 float ta=*(float*)a;
165 float tb=*(float*)b;
166 int out=(tb > ta) - (tb < ta);
167 return out ? out : COMPARE_FLOAT_POSTPROCESS;
168 }
169
170 int
gal_qsort_float32_i(const void * a,const void * b)171 gal_qsort_float32_i(const void *a, const void *b)
172 {
173 float ta=*(float*)a;
174 float tb=*(float*)b;
175 int out=(ta > tb) - (ta < tb);
176 return out ? out : COMPARE_FLOAT_POSTPROCESS;
177 }
178
179 int
gal_qsort_float64_d(const void * a,const void * b)180 gal_qsort_float64_d(const void *a, const void *b)
181 {
182 double ta=*(double*)a;
183 double tb=*(double*)b;
184 int out=(tb > ta) - (tb < ta);
185 return out ? out : COMPARE_FLOAT_POSTPROCESS;
186 }
187
188 int
gal_qsort_float64_i(const void * a,const void * b)189 gal_qsort_float64_i(const void *a, const void *b)
190 {
191 double ta=*(double*)a;
192 double tb=*(double*)b;
193 int out=(ta > tb) - (ta < tb);
194 return out ? out : COMPARE_FLOAT_POSTPROCESS;
195 }
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216 /*****************************************************************/
217 /*************** Sorting indexs ******************/
218 /*****************************************************************/
219 /* Initialize the array for sorting indexs to NULL. */
220 void *gal_qsort_index_single=NULL;
221
222 int
gal_qsort_index_single_uint8_d(const void * a,const void * b)223 gal_qsort_index_single_uint8_d(const void *a, const void *b)
224 {
225 uint8_t ta=((uint8_t *)(gal_qsort_index_single))[ *(size_t *)a ];
226 uint8_t tb=((uint8_t *)(gal_qsort_index_single))[ *(size_t *)b ];
227 return (tb > ta) - (tb < ta);
228 }
229
230 int
gal_qsort_index_single_uint8_i(const void * a,const void * b)231 gal_qsort_index_single_uint8_i(const void *a, const void *b)
232 {
233 uint8_t ta=((uint8_t *)(gal_qsort_index_single))[ *(size_t *)a ];
234 uint8_t tb=((uint8_t *)(gal_qsort_index_single))[ *(size_t *)b ];
235 return (ta > tb) - (ta < tb);
236 }
237
238 int
gal_qsort_index_single_int8_d(const void * a,const void * b)239 gal_qsort_index_single_int8_d(const void *a, const void *b)
240 {
241 int8_t ta=((int8_t *)(gal_qsort_index_single))[ *(size_t *)a ];
242 int8_t tb=((int8_t *)(gal_qsort_index_single))[ *(size_t *)b ];
243 return (tb > ta) - (tb < ta);
244 }
245
246 int
gal_qsort_index_single_int8_i(const void * a,const void * b)247 gal_qsort_index_single_int8_i(const void *a, const void *b)
248 {
249 int8_t ta=((int8_t *)(gal_qsort_index_single))[ *(size_t *)a ];
250 int8_t tb=((int8_t *)(gal_qsort_index_single))[ *(size_t *)b ];
251 return (ta > tb) - (ta < tb);
252 }
253
254 int
gal_qsort_index_single_uint16_d(const void * a,const void * b)255 gal_qsort_index_single_uint16_d(const void *a, const void *b)
256 {
257 uint16_t ta=((uint16_t *)(gal_qsort_index_single))[ *(size_t *)a ];
258 uint16_t tb=((uint16_t *)(gal_qsort_index_single))[ *(size_t *)b ];
259 return (tb > ta) - (tb < ta);
260 }
261
262 int
gal_qsort_index_single_uint16_i(const void * a,const void * b)263 gal_qsort_index_single_uint16_i(const void *a, const void *b)
264 {
265 uint16_t ta=((uint16_t *)(gal_qsort_index_single))[ *(size_t *)a ];
266 uint16_t tb=((uint16_t *)(gal_qsort_index_single))[ *(size_t *)b ];
267 return (ta > tb) - (ta < tb);
268 }
269
270 int
gal_qsort_index_single_int16_d(const void * a,const void * b)271 gal_qsort_index_single_int16_d(const void *a, const void *b)
272 {
273 int16_t ta=((int16_t *)(gal_qsort_index_single))[ *(size_t *)a ];
274 int16_t tb=((int16_t *)(gal_qsort_index_single))[ *(size_t *)b ];
275 return (tb > ta) - (tb < ta);
276 }
277
278 int
gal_qsort_index_single_int16_i(const void * a,const void * b)279 gal_qsort_index_single_int16_i(const void *a, const void *b)
280 {
281 int16_t ta=((int16_t *)(gal_qsort_index_single))[ *(size_t *)a ];
282 int16_t tb=((int16_t *)(gal_qsort_index_single))[ *(size_t *)b ];
283 return (ta > tb) - (ta < tb);
284 }
285
286 int
gal_qsort_index_single_uint32_d(const void * a,const void * b)287 gal_qsort_index_single_uint32_d(const void *a, const void *b)
288 {
289 uint32_t ta=((uint32_t *)(gal_qsort_index_single))[ *(size_t *)a ];
290 uint32_t tb=((uint32_t *)(gal_qsort_index_single))[ *(size_t *)b ];
291 return (tb > ta) - (tb < ta);
292 }
293
294 int
gal_qsort_index_single_uint32_i(const void * a,const void * b)295 gal_qsort_index_single_uint32_i(const void *a, const void *b)
296 {
297 uint32_t ta=((uint32_t *)(gal_qsort_index_single))[ *(size_t *)a ];
298 uint32_t tb=((uint32_t *)(gal_qsort_index_single))[ *(size_t *)b ];
299 return (ta > tb) - (ta < tb);
300 }
301
302 int
gal_qsort_index_single_int32_d(const void * a,const void * b)303 gal_qsort_index_single_int32_d(const void *a, const void *b)
304 {
305 int32_t ta=((int32_t *)(gal_qsort_index_single))[ *(size_t *)a ];
306 int32_t tb=((int32_t *)(gal_qsort_index_single))[ *(size_t *)b ];
307 return (tb > ta) - (tb < ta);
308 }
309
310 int
gal_qsort_index_single_int32_i(const void * a,const void * b)311 gal_qsort_index_single_int32_i(const void *a, const void *b)
312 {
313 int32_t ta=((int32_t *)(gal_qsort_index_single))[ *(size_t *)a ];
314 int32_t tb=((int32_t *)(gal_qsort_index_single))[ *(size_t *)b ];
315 return (ta > tb) - (ta < tb);
316 }
317
318 int
gal_qsort_index_single_uint64_d(const void * a,const void * b)319 gal_qsort_index_single_uint64_d(const void *a, const void *b)
320 {
321 uint64_t ta=((uint64_t *)(gal_qsort_index_single))[ *(size_t *)a ];
322 uint64_t tb=((uint64_t *)(gal_qsort_index_single))[ *(size_t *)b ];
323 return (tb > ta) - (tb < ta);
324 }
325
326 int
gal_qsort_index_single_uint64_i(const void * a,const void * b)327 gal_qsort_index_single_uint64_i(const void *a, const void *b)
328 {
329 uint64_t ta=((uint64_t *)(gal_qsort_index_single))[ *(size_t *)a ];
330 uint64_t tb=((uint64_t *)(gal_qsort_index_single))[ *(size_t *)b ];
331 return (ta > tb) - (ta < tb);
332 }
333
334 int
gal_qsort_index_single_int64_d(const void * a,const void * b)335 gal_qsort_index_single_int64_d(const void *a, const void *b)
336 {
337 int64_t ta=((int64_t *)(gal_qsort_index_single))[ *(size_t *)a ];
338 int64_t tb=((int64_t *)(gal_qsort_index_single))[ *(size_t *)b ];
339 return (tb > ta) - (tb < ta);
340 }
341
342 int
gal_qsort_index_single_int64_i(const void * a,const void * b)343 gal_qsort_index_single_int64_i(const void *a, const void *b)
344 {
345 int64_t ta=((int64_t *)(gal_qsort_index_single))[ *(size_t *)a ];
346 int64_t tb=((int64_t *)(gal_qsort_index_single))[ *(size_t *)b ];
347 return (ta > tb) - (ta < tb);
348 }
349
350 int
gal_qsort_index_single_float32_d(const void * a,const void * b)351 gal_qsort_index_single_float32_d(const void *a, const void *b)
352 {
353 float ta=((float *)(gal_qsort_index_single))[ *(size_t *)a ];
354 float tb=((float *)(gal_qsort_index_single))[ *(size_t *)b ];
355 int out=(tb > ta) - (tb < ta);
356 return out ? out : COMPARE_FLOAT_POSTPROCESS;
357 }
358
359 int
gal_qsort_index_single_float32_i(const void * a,const void * b)360 gal_qsort_index_single_float32_i(const void *a, const void *b)
361 {
362 float ta=((float *)(gal_qsort_index_single))[ *(size_t *)a ];
363 float tb=((float *)(gal_qsort_index_single))[ *(size_t *)b ];
364 int out=(ta > tb) - (ta < tb);
365 return out ? out : COMPARE_FLOAT_POSTPROCESS;
366 }
367
368 int
gal_qsort_index_single_float64_d(const void * a,const void * b)369 gal_qsort_index_single_float64_d(const void *a, const void *b)
370 {
371 double ta=((double *)(gal_qsort_index_single))[ *(size_t *)a ];
372 double tb=((double *)(gal_qsort_index_single))[ *(size_t *)b ];
373 int out=(tb > ta) - (tb < ta);
374 return out ? out : COMPARE_FLOAT_POSTPROCESS;
375 }
376
377 int
gal_qsort_index_single_float64_i(const void * a,const void * b)378 gal_qsort_index_single_float64_i(const void *a, const void *b)
379 {
380 double ta=((double *)(gal_qsort_index_single))[ *(size_t *)a ];
381 double tb=((double *)(gal_qsort_index_single))[ *(size_t *)b ];
382 int out=(ta > tb) - (ta < tb);
383 return out ? out : COMPARE_FLOAT_POSTPROCESS;
384 }
385
386 int
gal_qsort_index_multi_d(const void * a,const void * b)387 gal_qsort_index_multi_d(const void *a, const void *b)
388 {
389 /* Define the structures. */
390 struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
391 struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
392
393 /* For easy reading. */
394 float ta=A->values[ A->index ];
395 float tb=B->values[ B->index ];
396
397 /* Return the result. */
398 int out=(tb > ta) - (tb < ta);
399 return out ? out : COMPARE_FLOAT_POSTPROCESS;
400 }
401
402 int
gal_qsort_index_multi_i(const void * a,const void * b)403 gal_qsort_index_multi_i(const void *a, const void *b)
404 {
405 /* Define the structures. */
406 struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
407 struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
408
409 /* For easy reading. */
410 float ta=A->values[ A->index ];
411 float tb=B->values[ B->index ];
412
413 /* Return the result. */
414 int out=(ta > tb) - (ta < tb);
415 return out ? out : COMPARE_FLOAT_POSTPROCESS;
416 }
417