NaN fixes in pow, powf, modf, modff

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

NaN fixes in pow, powf, modf, modff

Sourceware - newlib list mailing list
The first patch makes modf set *iptr to qnan when provided an snan
parameter, which wasn't happening on ARM softfp because GCC elided the
* 1.0 operation.

The second patch cleans up modf/modff by not re-converting the float
to bits before extracting the sign bit.

The third patch makes pow return qnan instead of 1.0 when the
parameter is snan.

These functions now match glibc for all of my tests, which were
derived from the IEEE requirements.


Reply | Threaded
Open this post in threaded view
|

[PATCH 1/3] newlib/libm/common: Fix modf/modff returning snan

Sourceware - newlib list mailing list
Recent GCC appears to elide multiplication by 1, which causes snan
parameters to be returned unchanged through *iptr. Use the existing
conversion of snan to qnan to also set the correct result in *iptr
instead.

Signed-off-by: Keith Packard <[hidden email]>
---
 newlib/libm/common/s_modf.c  | 10 ++--------
 newlib/libm/common/sf_modf.c | 10 ++--------
 2 files changed, 4 insertions(+), 16 deletions(-)

diff --git a/newlib/libm/common/s_modf.c b/newlib/libm/common/s_modf.c
index c948b8525..c826580b4 100644
--- a/newlib/libm/common/s_modf.c
+++ b/newlib/libm/common/s_modf.c
@@ -63,12 +63,6 @@ QUICKREF
 
 #ifndef _DOUBLE_IS_32BITS
 
-#ifdef __STDC__
-static const double one = 1.0;
-#else
-static double one = 1.0;
-#endif
-
 #ifdef __STDC__
  double modf(double x, double *iptr)
 #else
@@ -99,8 +93,8 @@ static double one = 1.0;
     }
  } else if (j0>51) { /* no fraction part */
     __uint32_t high;
-    *iptr = x*one;
-    if (__fpclassifyd(x) == FP_NAN) return x+x; /* x is NaN, return NaN */
+    *iptr = x;
+    if (__fpclassifyd(x) == FP_NAN) return *iptr = x+x; /* x is NaN, return NaN */
     GET_HIGH_WORD(high,x);
     INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
     return x;
diff --git a/newlib/libm/common/sf_modf.c b/newlib/libm/common/sf_modf.c
index ae970762b..e241e4612 100644
--- a/newlib/libm/common/sf_modf.c
+++ b/newlib/libm/common/sf_modf.c
@@ -15,12 +15,6 @@
 
 #include "fdlibm.h"
 
-#ifdef __STDC__
-static const float one = 1.0;
-#else
-static float one = 1.0;
-#endif
-
 #ifdef __STDC__
  float modff(float x, float *iptr)
 #else
@@ -51,8 +45,8 @@ static float one = 1.0;
     }
  } else { /* no fraction part */
     __uint32_t ix;
-    *iptr = x*one;
-    if (__fpclassifyf(x) == FP_NAN) return x+x; /* x is NaN, return NaN */
+    *iptr = x;
+    if (__fpclassifyf(x) == FP_NAN) return *iptr = x+x; /* x is NaN, return NaN */
     GET_FLOAT_WORD(ix,x);
     SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */
     return x;
--
2.25.1

Reply | Threaded
Open this post in threaded view
|

[PATCH 2/3] newlib/libm/common: Don't re-convert float to bits in modf/modff

Sourceware - newlib list mailing list
In reply to this post by Sourceware - newlib list mailing list
These functions shared a pattern of re-converting the argument to bits
when returning +/-0. Skip that as the initial conversion still has the
sign bit.

Signed-off-by: Keith Packard <[hidden email]>
---
 newlib/libm/common/s_modf.c  | 12 +++---------
 newlib/libm/common/sf_modf.c |  8 ++------
 2 files changed, 5 insertions(+), 15 deletions(-)

diff --git a/newlib/libm/common/s_modf.c b/newlib/libm/common/s_modf.c
index c826580b4..e552a9460 100644
--- a/newlib/libm/common/s_modf.c
+++ b/newlib/libm/common/s_modf.c
@@ -81,10 +81,8 @@ QUICKREF
     } else {
  i = (0x000fffff)>>j0;
  if(((i0&i)|i1)==0) { /* x is integral */
-    __uint32_t high;
     *iptr = x;
-    GET_HIGH_WORD(high,x);
-    INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
+    INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
     return x;
  } else {
     INSERT_WORDS(*iptr,i0&(~i),0);
@@ -92,19 +90,15 @@ QUICKREF
  }
     }
  } else if (j0>51) { /* no fraction part */
-    __uint32_t high;
     *iptr = x;
     if (__fpclassifyd(x) == FP_NAN) return *iptr = x+x; /* x is NaN, return NaN */
-    GET_HIGH_WORD(high,x);
-    INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
+    INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
     return x;
  } else { /* fraction part in low x */
     i = ((__uint32_t)(0xffffffff))>>(j0-20);
     if((i1&i)==0) { /* x is integral */
-        __uint32_t high;
  *iptr = x;
- GET_HIGH_WORD(high,x);
- INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
+ INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
  return x;
     } else {
         INSERT_WORDS(*iptr,i0,i1&(~i));
diff --git a/newlib/libm/common/sf_modf.c b/newlib/libm/common/sf_modf.c
index e241e4612..2994378bb 100644
--- a/newlib/libm/common/sf_modf.c
+++ b/newlib/libm/common/sf_modf.c
@@ -33,10 +33,8 @@
     } else {
  i = (0x007fffff)>>j0;
  if((i0&i)==0) { /* x is integral */
-    __uint32_t ix;
     *iptr = x;
-    GET_FLOAT_WORD(ix,x);
-    SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */
+    SET_FLOAT_WORD(x,i0&0x80000000); /* return +-0 */
     return x;
  } else {
     SET_FLOAT_WORD(*iptr,i0&(~i));
@@ -44,11 +42,9 @@
  }
     }
  } else { /* no fraction part */
-    __uint32_t ix;
     *iptr = x;
     if (__fpclassifyf(x) == FP_NAN) return *iptr = x+x; /* x is NaN, return NaN */
-    GET_FLOAT_WORD(ix,x);
-    SET_FLOAT_WORD(x,ix&0x80000000); /* return +-0 */
+    SET_FLOAT_WORD(x,i0&0x80000000); /* return +-0 */
     return x;
  }
 }
--
2.25.1

Reply | Threaded
Open this post in threaded view
|

[PATCH 3/3] newlib/libm/math: Make pow/powf return qnan for snan arg

Sourceware - newlib list mailing list
In reply to this post by Sourceware - newlib list mailing list
The IEEE spec for pow only has special case for x**0 and 1**y when x/y
are quiet NaN. For signaling NaN, the general case applies and these functions
should signal the invalid exception and return a quiet NaN.

Signed-off-by: Keith Packard <[hidden email]>
---
 newlib/libm/math/e_pow.c  | 13 +++++++++----
 newlib/libm/math/ef_pow.c | 10 +++++++---
 2 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index 5fd28e65f..4c450ec05 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -58,6 +58,8 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
+
 #if __OBSOLETE_MATH
 
 #ifndef _DOUBLE_IS_32BITS
@@ -116,14 +118,17 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
  EXTRACT_WORDS(hy,ly,y);
  ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
-    /* y==zero: x**0 = 1 */
- if((iy|ly)==0) return one;
+    /* y==zero: x**0 = 1 unless x is snan */
+ if((iy|ly)==0) {
+    if (issignaling_inline(x)) return x + y;
+    return one;
+ }
 
     /* x|y==NaN return NaN unless x==1 then return 1 */
  if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
    iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) {
-    if(((hx-0x3ff00000)|lx)==0) return one;
-    else return nan("");
+    if(((hx-0x3ff00000)|lx)==0 && !issignaling_inline(y)) return one;
+    else return x + y;
  }
 
     /* determine if y is an odd int when x < 0
diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
index d9e85a95e..d4ea4c5e8 100644
--- a/newlib/libm/math/ef_pow.c
+++ b/newlib/libm/math/ef_pow.c
@@ -14,6 +14,7 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #if __OBSOLETE_MATH
 #ifdef __v810__
@@ -74,13 +75,16 @@ ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
  ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
- if(FLT_UWORD_IS_ZERO(iy)) return one;
+ if(FLT_UWORD_IS_ZERO(iy)) {
+    if (issignalingf_inline(x)) return x + y;
+    return one;
+ }
 
     /* x|y==NaN return NaN unless x==1 then return 1 */
  if(FLT_UWORD_IS_NAN(ix) ||
    FLT_UWORD_IS_NAN(iy)) {
-    if(hx==0x3f800000) return one;
-    else return nanf("");
+    if(hx==0x3f800000 && !issignalingf_inline(y)) return one;
+    else return x + y;
  }
 
     /* determine if y is an odd int when x < 0
--
2.25.1

Reply | Threaded
Open this post in threaded view
|

Re: NaN fixes in pow, powf, modf, modff

Sourceware - newlib list mailing list
In reply to this post by Sourceware - newlib list mailing list
On Mar 25 17:18, Keith Packard via Newlib wrote:

> The first patch makes modf set *iptr to qnan when provided an snan
> parameter, which wasn't happening on ARM softfp because GCC elided the
> * 1.0 operation.
>
> The second patch cleans up modf/modff by not re-converting the float
> to bits before extracting the sign bit.
>
> The third patch makes pow return qnan instead of 1.0 when the
> parameter is snan.
>
> These functions now match glibc for all of my tests, which were
> derived from the IEEE requirements.
>
Pushed.


Thanks,
Corinna

--
Corinna Vinschen
Cygwin Maintainer
Red Hat

signature.asc (849 bytes) Download Attachment