totalorder, totalorderf: Fix handling of SNaN on i386 and x86_64 CPUs

[Bruno Haible]

On NetBSD 10.0/i386 and OpenBSD 7.5/i386, I see these test failures:

FAIL: test-totalorder
=====================

Failed: i=0 j=1
Failed: i=1 j=0
Failed: i=12 j=13
Failed: i=13 j=12
FAIL test-totalorder (exit status: 1)

FAIL: test-totalorderf
======================

Failed: i=0 j=1
Failed: i=1 j=0
Failed: i=12 j=13
Failed: i=13 j=12
FAIL test-totalorderf (exit status: 1)

I can reproduce them also on Linux, with glibc, by use of a testdir for the
modules
  totalorder
  totalorderf
  totalorderl
configured with
  gl_cv_func_totalorder_no_libm=no gl_cv_func_totalorder_in_libm=no \
  gl_cv_func_totalorderf_no_libm=no gl_cv_func_totalorderf_in_libm=no \
  gl_cv_func_totalorderl_no_libm=no gl_cv_func_totalorderl_in_libm=no \
  CFLAGS=-ggdb ./configure
both
  - on i386,
  - on x86_64 with CC="gcc -mfpmath=387".

The problem is that in the lines
54        xu.f = *x;
55        yu.f = *y;

an SNaN value gets converted to a quiet NaN, through an 'flds' or
'fldl' instruction, respectively.

The long explanation is here:
  https://lists.gnu.org/archive/html/bug-gnulib/2023-10/msg00060.html
We need to follow the rule

  Programs that need to distinguish QNaN and SNaN, but don't do
  floating-point operations:
  In these cases it is sufficient to
    - either pass values by reference, not by value,
    - or apply the 'union' workaround.

This patch does it and thus fixes the test failures on Linux and NetBSD.


2024-04-09  Bruno Haible  <bruno@clisp.org>

        totalorder, totalorderf: Fix handling of SNaN on i386 and x86_64 CPUs.
        * lib/totalorder.c: Include <string.h>.
        (totalorder): Use memcpy to copy the 'double' values into the union.
        Drop 'volatile'.
        * lib/totalorderf.c: Include <string.h>.
        (totalorderf): Use memcpy to copy the 'float' values into the union.
        Drop 'volatile'.

diff --git a/lib/totalorder.c b/lib/totalorder.c
index 1fb8f0d24d..635e3cb276 100644
--- a/lib/totalorder.c
+++ b/lib/totalorder.c
@@ -21,6 +21,8 @@
 /* Specification.  */
 #include <math.h>
 
+#include <string.h>
+
 int
 totalorder (double const *x, double const *y)
 {
@@ -50,8 +52,18 @@ totalorder (double const *x, double const *y)
   /* Invert the most significant bit of the mantissa field.  Cf. snan.h.  */
   extended_sign ^= (1ULL << 51);
 #endif
-  union { unsigned long long i; double f; } volatile xu = {0}, yu = {0};
+  union { unsigned long long i; double f; } xu = {0}, yu = {0};
+#if 0
   xu.f = *x;
   yu.f = *y;
+#else
+  /* On 32-bit x86 processors, as well as on x86_64 processors with
+     CC="gcc -mfpmath=387", the evaluation of *x and *y above is done through
+     an 'fldl' instruction, which converts a signalling NaN to a quiet NaN. See
+     <https://lists.gnu.org/archive/html/bug-gnulib/2023-10/msg00060.html>
+     for details.  Use memcpy to avoid this.  */
+  memcpy (&xu.f, x, sizeof (double));
+  memcpy (&yu.f, y, sizeof (double));
+#endif
   return (xu.i ^ extended_sign) <= (yu.i ^ extended_sign);
 }
diff --git a/lib/totalorderf.c b/lib/totalorderf.c
index d235bc98bc..75024b6839 100644
--- a/lib/totalorderf.c
+++ b/lib/totalorderf.c
@@ -21,6 +21,8 @@
 /* Specification.  */
 #include <math.h>
 
+#include <string.h>
+
 int
 totalorderf (float const *x, float const *y)
 {
@@ -50,8 +52,18 @@ totalorderf (float const *x, float const *y)
   /* Invert the most significant bit of the mantissa field.  Cf. snan.h.  */
   extended_sign ^= (1U << 22);
 #endif
-  union { unsigned int i; float f; } volatile xu = {0}, yu = {0};
+  union { unsigned int i; float f; } xu = {0}, yu = {0};
+#if 0
   xu.f = *x;
   yu.f = *y;
+#else
+  /* On 32-bit x86 processors, as well as on x86_64 processors with
+     CC="gcc -mfpmath=387", the evaluation of *x and *y above is done through
+     an 'flds' instruction, which converts a signalling NaN to a quiet NaN. See
+     <https://lists.gnu.org/archive/html/bug-gnulib/2023-10/msg00060.html>
+     for details.  Use memcpy to avoid this.  */
+  memcpy (&xu.f, x, sizeof (float));
+  memcpy (&yu.f, y, sizeof (float));
+#endif
   return (xu.i ^ extended_sign) <= (yu.i ^ extended_sign);
 }