[gnuastro-commits] master 85a7991 3/3: qsort library: floating point NaN values placed in the end

[Mohammad Akhlaghi]

branch: master
commit 85a79914a1f75a17a2ba2a1cbaac64dcae00888a
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    qsort library: floating point NaN values placed in the end
    
    Until now, all the `gal_qsort' functions that worked with floating point
    types simply returned the standard check (like `(ta > tb) - (ta < tb)' for
    an increasing sort). However, if any, or both, of the elements were NaN, it
    would return 0 (meaning equal). This would create an incorrect sorted
    result.
    
    To fix the problem, with this commit a `COMPARE_FLOAT_POSTPROCESS' macro
    has been added in `lib/qsort.c' and all floating point comparison functions
    will use it when the result of the basic comparison is 0.
    
    This bug was reported by Roberto Baena Gallé.
    
    This fixes bug #56671.
---
 NEWS              |  1 +
 doc/gnuastro.texi |  9 +++++++++
 lib/qsort.c       | 57 ++++++++++++++++++++++++++++++++++++++++++++-----------
 3 files changed, 56 insertions(+), 11 deletions(-)

diff --git a/NEWS b/NEWS
index 98b850b..022d54d 100644
--- a/NEWS
+++ b/NEWS
@@ -147,6 +147,7 @@ See the end of the file for license conditions.
   bug #56641: MakeProfile's center position changes based on precision.
   bug #56635: Update tutorial 3 with bug-fixed NoiseChisel.
   bug #56662: Converting -R to -Wl,-R causes a crash in configure on macOS.
+  bug #56671: Bad sorting with asttable if nan is present.
 
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index bd49d1c..32a1961 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -29210,6 +29210,7 @@ both polygons have to be sorted in an anti-clock-wise 
manner.
 @node Qsort functions, Permutations, Polygons, Gnuastro library
 @subsection Qsort functions (@file{qsort.h})
 
+@cindex @code{qsort}
 When sorting a dataset is necessary, the C programming language provides
 the @code{qsort} (Quick sort) function. @code{qsort} is a generic function
 which allows you to sort any kind of data structure (not just a single
@@ -29217,6 +29218,14 @@ array of numbers). To define ``greater'' and 
``smaller'' (for sorting),
 @code{qsort} needs another function, even for simple numerical types. The
 functions introduced in this section are to passed onto @code{qsort}.
 
+@cindex NaN
+Note that larger and smaller operators are not defined on NaN
+elements. Therefore, if the input array is a floating point type, and
+contains NaN values, the relevant functions of this section are going to
+put the NaN elements at the end of the list (after the sorted non-NaN
+elements), irrespective of the requested sorting order (increasing or
+decreasing).
+
 The first class of functions below (with @code{TYPE} in their names) can be
 used for sorting a simple numeric array. Just replace @code{TYPE} with the
 dataset's numeric datatype. The second set of functions can be used to sort
diff --git a/lib/qsort.c b/lib/qsort.c
index 3bba8c6..20aa253 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -22,6 +22,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 **********************************************************************/
 #include <config.h>
 
+#include <math.h>
 #include <stdlib.h>
 #include <stdint.h>
 
@@ -30,6 +31,31 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/qsort.h>
 
 
+/*****************************************************************/
+/**********                  Macros               ****************/
+/*****************************************************************/
+/* When one or both elements are NaN, the simple comparison, like `(tb >
+   ta) - (tb < ta)', will give 0 (as if the elements are equal). However,
+   some preference has to be given to the NaN element in a comparison,
+   otherwise the output is not going to be reasonable. We also don't want
+   to check NaNs on every comparison (it will slow down the processing).
+
+   So we'll exploit the fact that when there comparison result doesn't
+   equal zero, we don't have any NaNs and this `COMPARE_FLOAT_POSTPROCESS'
+   macro is called only when the comparison gives zero. Being larger or
+   smaller isn't defined for NaNs, so we'll just put them in the end of the
+   sorted list whether it is sorted by decreasing or increasing mode.*/
+#define COMPARE_FLOAT_POSTPROCESS (                                     \
+   isnan(ta) && isnan(tb)                                               \
+   ? 0                                /* Both NaN, define as equal. */  \
+   /* One is NaN, one isn't. */                                         \
+   : ( isnan(ta)                                                        \
+       ? 1                         /* First is NaN, set as smaller. */  \
+       : ( isnan(tb)                                                    \
+           ? -1                    /* Second is NaN, set as larger. */  \
+           : 0 )                      /* None are NaN, set as equal.*/  \
+       )                                                                \
+)
 
 
 
@@ -120,7 +146,6 @@ gal_qsort_uint64_i(const void *a, const void *b)
   return ( *(uint64_t *)a - *(uint64_t *)b );
 }
 
-
 int
 gal_qsort_int64_d(const void *a, const void *b)
 {
@@ -138,7 +163,8 @@ gal_qsort_float32_d(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
-  return (tb > ta) - (tb < ta);
+  int out=(tb > ta) - (tb < ta);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -146,7 +172,8 @@ gal_qsort_float32_i(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
-  return (ta > tb) - (ta < tb);
+  int out=(ta > tb) - (ta < tb);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -154,7 +181,8 @@ gal_qsort_float64_d(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
-  return (tb > ta) - (tb < ta);
+  int out=(tb > ta) - (tb < ta);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -162,7 +190,8 @@ gal_qsort_float64_i(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
-  return (ta > tb) - (ta < tb);
+  int out=(ta > tb) - (ta < tb);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 
@@ -323,7 +352,8 @@ gal_qsort_index_single_float32_d(const void *a, const void 
*b)
 {
   float ta=((float *)(gal_qsort_index_single))[ *(size_t *)a ];
   float tb=((float *)(gal_qsort_index_single))[ *(size_t *)b ];
-  return (tb > ta) - (tb < ta);
+  int out=(tb > ta) - (tb < ta);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -331,7 +361,8 @@ gal_qsort_index_single_float32_i(const void *a, const void 
*b)
 {
   float ta=((float *)(gal_qsort_index_single))[ *(size_t *)a ];
   float tb=((float *)(gal_qsort_index_single))[ *(size_t *)b ];
-  return (ta > tb) - (ta < tb);
+  int out=(ta > tb) - (ta < tb);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -339,7 +370,8 @@ gal_qsort_index_single_float64_d(const void *a, const void 
*b)
 {
   double ta=((double *)(gal_qsort_index_single))[ *(size_t *)a ];
   double tb=((double *)(gal_qsort_index_single))[ *(size_t *)b ];
-  return (tb > ta) - (tb < ta);
+  int out=(tb > ta) - (tb < ta);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -347,7 +379,8 @@ gal_qsort_index_single_float64_i(const void *a, const void 
*b)
 {
   double ta=((double *)(gal_qsort_index_single))[ *(size_t *)a ];
   double tb=((double *)(gal_qsort_index_single))[ *(size_t *)b ];
-  return (ta > tb) - (ta < tb);
+  int out=(ta > tb) - (ta < tb);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -362,7 +395,8 @@ gal_qsort_index_multi_d(const void *a, const void *b)
   float tb=B->values[ B->index ];
 
   /* Return the result. */
-  return (tb > ta) - (tb < ta);
+  int out=(tb > ta) - (tb < ta);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }
 
 int
@@ -377,5 +411,6 @@ gal_qsort_index_multi_i(const void *a, const void *b)
   float tb=B->values[ B->index ];
 
   /* Return the result. */
-  return (ta > tb) - (ta < tb);
+  int out=(ta > tb) - (ta < tb);
+  return out ? out : COMPARE_FLOAT_POSTPROCESS;
 }