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