Re: mathl and NaN

[Bruno Haible]

Paolo Bonzini wrote:
> I think ==, !=, isordered, isunordered comparisons are
> the sole operations you can safely perform on a SNaN.

Even these operations will trap on a SNaN, according to [1]:

   Signaling NaNs shall be reserved operands that signal the invalid
   operation exception (7.1) for every operation listed in
   Section 5. Whether copying a signaling NaN without a change of
   format signals the invalid operation exception is the implementor's
   option.

   (And section 5.7 contains the comparison operations.)

Kahan [2] writes this about SNaNs:

   An SNaN may be moved ( copied ) without incident, but any other arithmetic
   operation upon an SNaN is an INVALID operation ( and so is loading one onto
   the ix87's stack ) that must trap or else produce a new nonsignaling NaN.
   ( Another way to turn an SNaN into a NaN is to turn 0xxx...xxx into
   1xxx...xxx with a logical OR.)

So, in functions that should trap when they get an sNaN argument, we need
to perform at least one operation before bailing out due to isnanl(x).

But in POSIX with MX extension, the functions are specified:
    If x is NaN, a NaN shall be returned.
So we need to write the isnanl(x) test as first thing, before any other
operations.

I'm applying the appended patch.

Bruno

[1] http://754r.ucbtest.org/standards/754.pdf
[2] http://www.cs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF


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

        Better support of signalling NaNs.
        * lib/atanl.c: Include isnanl.h.
        (atanl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/cosl.c: Include isnanl.h.
        (cosl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/ldexpl.c: Include isnanl.h.
        (ldexpl): Perform test for NaN through a call to isnanl.
        * lib/logl.c: Include isnanl.h.
        (logl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/sinl.c: Include isnanl.h.
        (sinl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/sqrtl.c: Include isnanl.h.
        (sqrtl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/tanl.c: Include isnanl.h.
        (tanl): Perform test for NaN at the beginning of the function and
        through a call to isnanl.
        * lib/trigl.c (ieee754_rem_pio2l): Remove test for NaN.

--- lib/atanl.c 18 Feb 2007 15:10:28 -0000      1.4
+++ lib/atanl.c 26 Mar 2007 22:18:33 -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/cosl.c  18 Feb 2007 15:10:28 -0000      1.3
+++ lib/cosl.c  26 Mar 2007 22:18:33 -0000
@@ -54,19 +54,24 @@
 #include "trigl.c"
 #include "sincosl.c"
 #endif
+#include "isnanl.h"
 
 long double cosl(long double x)
 {
        long double y[2],z=0.0L;
        int n;
 
+    /* cosl(NaN) is NaN */
+        if (isnanl (x))
+          return x;
+
     /* |x| ~< pi/4 */
         if(x >= -0.7853981633974483096156608458198757210492 &&
            x <= 0.7853981633974483096156608458198757210492)
           return kernel_cosl(x, z);
 
-    /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */
-        else if ((x + x == x && x != 0.0) || x != x)
+    /* cosl(Inf) is NaN, cosl(0) is 1 */
+        else if (x + x == x && x != 0.0)
           return x-x;           /* NaN */
 
     /* argument reduction needed */
--- lib/ldexpl.c        18 Feb 2007 15:10:28 -0000      1.5
+++ lib/ldexpl.c        26 Mar 2007 22:18:33 -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  26 Mar 2007 22:18:33 -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 */
@@ -196,6 +198,11 @@
 
   /* Check for IEEE special cases.  */
 
+  /* log(NaN) = NaN. */
+  if (isnanl (x))
+    {
+      return x;
+    }
   /* log(0) = -infinity. */
   if (x == 0.0L)
     {
@@ -206,8 +213,8 @@
     {
       return (x - x) / ZERO;
     }
-  /* log (infinity or NaN) */
-  if (x + x == x || x != x)
+  /* log (infinity) */
+  if (x + x == x)
     {
       return x + x;
     }
--- lib/sinl.c  18 Feb 2007 15:10:28 -0000      1.3
+++ lib/sinl.c  26 Mar 2007 22:18:33 -0000
@@ -52,6 +52,7 @@
 #include "trigl.h"
 #include "trigl.c"
 #include "sincosl.c"
+#include "isnanl.h"
 
 long double
 sinl (long double x)
@@ -59,13 +60,17 @@
   long double y[2], z = 0.0L;
   int n;
 
+  /* sinl(NaN) is NaN */
+  if (isnanl (x))
+    return x;
+
   /* |x| ~< pi/4 */
   if (x >= -0.7853981633974483096156608458198757210492 &&
       x <= 0.7853981633974483096156608458198757210492)
     return kernel_sinl (x, z, 0);
 
-    /* sinl(Inf or NaN) is NaN, sinl(0) is 0 */
-  else if (x + x == x || x != x)
+    /* sinl(Inf) is NaN, sinl(0) is 0 */
+  else if (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 26 Mar 2007 22:18:33 -0000
@@ -25,6 +25,7 @@
 #include <math.h>
 
 #include <float.h>
+#include "isnanl.h"
 
 /* A simple Newton-Raphson method. */
 long double
@@ -33,12 +34,16 @@
   long double delta, y;
   int exponent;
 
+  /* Check for NaN */
+  if (isnanl (x))
+    return x;
+
   /* Check for negative numbers */
   if (x < 0.0L)
     return (long double) sqrt(-1);
 
-  /* Check for zero, NANs and infinites */
-  if (x + x == x || x != x)
+  /* Check for zero and infinites */
+  if (x + x == x)
     return x;
 
   frexpl (x, &exponent);
--- lib/tanl.c  25 Mar 2007 20:36:17 -0000      1.4
+++ lib/tanl.c  26 Mar 2007 22:18:33 -0000
@@ -55,6 +55,7 @@
 #include "trigl.c"
 #endif
 #endif
+#include "isnanl.h"
 
 /*
  * ====================================================
@@ -191,13 +192,17 @@
   long double y[2], z = 0.0L;
   int n;
 
+  /* tanl(NaN) is NaN */
+  if (isnanl (x))
+    return x;
+
   /* |x| ~< pi/4 */
   if (x >= -0.7853981633974483096156608458198757210492 &&
       x <= 0.7853981633974483096156608458198757210492)
     return kernel_tanl (x, z, 1);
 
-  /* tanl(Inf or NaN) is NaN, tanl(0) is 0 */
-  else if (x + x == x || x != x)
+  /* tanl(Inf) is NaN, tanl(0) is 0 */
+  else if (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 26 Mar 2007 22:18:34 -0000
@@ -233,7 +233,7 @@
        return -1;
       }
 
-  if (x + x == x || x != x)    /* x is +=oo or NaN */
+  if (x + x == x)      /* x is ±oo */
     {
       y[0] = x - x;
       y[1] = y[0];