mathl and NaN

[Bruno Haible]

Hi Paolo,

The mathl module uses the idiom   x != x
to test for a NaN. The IRIX 6.5 cc compiler, however, simplifies this expression
to false, both for 'double' and 'long double' variables. The two workarounds
that I found are:
  - Use the isnanl function,
  - Copy one of the operands into a 'volatile' variable:
       volatile double y = x;
       if (x != y) ...

Also, do you know which are the operations that can be performed on a
signalling NaN without getting an exception? Only an isnanl test? Or also
comparisons?

So, your opinion about the following patch?

2007-03-25  Bruno Haible  <address@hidden>

        * lib/atanl.c: Include isnanl.h.
        (atanl): Use isnanl, instead of the x != x idiom.
        * lib/ldexpl.c: Include isnanl.h.
        (ldexpl): Use isnanl, instead of the x != x idiom.
        * lib/logl.c: Include isnanl.h.
        (logl): Use isnanl, instead of the x != x idiom.
        * lib/sinl.c: Include isnanl.h.
        (sinl): Use isnanl, instead of the x != x idiom.
        * lib/sqrtl.c: Include isnanl.h.
        (sqrtl): Use isnanl, instead of the x != x idiom.
        * lib/tanl.c: Include isnanl.h.
        (tanl): Use isnanl, instead of the x != x idiom.
        * lib/trigl.c: Include isnanl.h.
        (ieee754_rem_pio2l): Use isnanl, instead of the x != x idiom.
        * modules/mathl (Depends-on): Add isnanl.

--- lib/atanl.c 18 Feb 2007 15:10:28 -0000      1.4
+++ lib/atanl.c 25 Mar 2007 20:34:50 -0000
@@ -64,6 +64,7 @@
  *
  */
 
+#include "isnanl.h"
 
 /* arctan(k/8), k = 0, ..., 82 */
 static const long double atantbl[84] = {
@@ -178,12 +179,12 @@
   int k, sign;
   long double t, u, p, q;
 
-  sign = x < 0.0;
-
   /* Check for zero or NaN.  */
-  if (x != x || x == 0.0)
+  if (isnanl (x) || x == 0.0)
     return x + x;
 
+  sign = x < 0.0;
+
   if (x + x == x)
     {
       /* Infinity. */
--- lib/ldexpl.c        18 Feb 2007 15:10:28 -0000      1.5
+++ lib/ldexpl.c        25 Mar 2007 20:34:50 -0000
@@ -25,6 +25,7 @@
 #include <math.h>
 
 #include <float.h>
+#include "isnanl.h"
 
 long double
 ldexpl(long double x, int exp)
@@ -33,7 +34,7 @@
   int bit;
 
   /* Check for zero, nan and infinity. */
-  if (x != x || x + x == x )
+  if (isnanl (x) || x + x == x)
     return x;
 
   if (exp < 0)
--- lib/logl.c  18 Feb 2007 15:10:28 -0000      1.4
+++ lib/logl.c  25 Mar 2007 20:34:50 -0000
@@ -64,6 +64,8 @@
  *
  */
 
+#include "isnanl.h"
+
 /* log(1+x) = x - .5 x^2 + x^3 l(x)
    -.0078125 <= x <= +.0078125
    peak relative error 1.2e-37 */
@@ -207,7 +209,7 @@
       return (x - x) / ZERO;
     }
   /* log (infinity or NaN) */
-  if (x + x == x || x != x)
+  if (isnanl (x) || x + x == x)
     {
       return x + x;
     }
--- lib/sinl.c  18 Feb 2007 15:10:28 -0000      1.3
+++ lib/sinl.c  25 Mar 2007 20:34:50 -0000
@@ -52,6 +52,7 @@
 #include "trigl.h"
 #include "trigl.c"
 #include "sincosl.c"
+#include "isnanl.h"
 
 long double
 sinl (long double x)
@@ -65,7 +66,7 @@
     return kernel_sinl (x, z, 0);
 
     /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */
-  else if (x + x == x || x != x)
+  else if (isnanl (x) || x + x == x)
     return x - x;              /* NaN */
 
   /* argument reduction needed */
--- lib/sqrtl.c 18 Feb 2007 15:10:28 -0000      1.4
+++ lib/sqrtl.c 25 Mar 2007 20:34:50 -0000
@@ -25,6 +25,7 @@
 #include <math.h>
 
 #include <float.h>
+#include "isnanl.h"
 
 /* A simple Newton-Raphson method. */
 long double
@@ -38,7 +39,7 @@
     return (long double) sqrt(-1);
 
   /* Check for zero, NANs and infinites */
-  if (x + x == x || x != x)
+  if (isnanl (x) || x + x == x)
     return x;
 
   frexpl (x, &exponent);
--- lib/tanl.c  18 Feb 2007 15:10:28 -0000      1.3
+++ lib/tanl.c  25 Mar 2007 20:34:50 -0000
@@ -55,6 +55,7 @@
 #include "trigl.c"
 #endif
 #endif
+#include "isnanl.h"
 
 /*
  * ====================================================
@@ -197,7 +198,7 @@
     return kernel_tanl (x, z, 1);
 
   /* tanl(Inf or NaN) is NaN, tanl(0) is 0 */
-  else if (x + x == x || x != x)
+  else if (isnanl (x) || x + x == x)
     return x - x;              /* NaN */
 
   /* argument reduction needed */
--- lib/trigl.c 25 Mar 2007 18:03:15 -0000      1.5
+++ lib/trigl.c 25 Mar 2007 20:34:50 -0000
@@ -23,6 +23,7 @@
 #include <math.h>
 
 #include <float.h>
+#include "isnanl.h"
 
 /* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */
 static const int two_over_pi[] = {
@@ -233,7 +234,7 @@
        return -1;
       }
 
-  if (x + x == x || x != x)    /* x is +=oo or NaN */
+  if (isnanl (x) || x + x == x)        /* x is +=oo or NaN */
     {
       y[0] = x - x;
       y[1] = y[0];
--- modules/mathl       23 Mar 2007 01:26:24 -0000      1.7
+++ modules/mathl       25 Mar 2007 20:34:50 -0000
@@ -22,6 +22,7 @@
 Depends-on:
 math
 frexpl
+isnanl
 
 configure.ac:
 gl_FUNC_LONG_DOUBLE_MATH