NaN fixes in pow, powf, modf, modff

NaN fixes in pow, powf, modf, modff

[Keith Packard]

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.

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

[Keith Packard]

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 <keithp@keithp.com>
---
 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

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

[Keith Packard]

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 <keithp@keithp.com>
---
 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

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

[Keith Packard]

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 <keithp@keithp.com>
---
 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

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

[Corinna Vinschen]

[-- Attachment #1: Type: text/plain, Size: 624 bytes --]

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

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 833 bytes --]