Re: [PATCH] libm/math: Use __math_xflow in obsolete math code

Keith Packard via Newlib <[email protected]> Sat, 01 Aug 2020 15:40:28 -0700

Newsgroups	gmane.comp.lib.newlib
Message-ID	<[email protected]>

Szabolcs Nagy <[email protected]> writes:

> this looks good.

Thanks.

> note1: i used c99 code when i wrote the
> new math code and currently the old math
> code can be compiled with older compilers,
> this change might prevent that.

I don't really have any way to test that as I'm using current compilers
for this work.

> note2: xflow may also try to set errno,
> depending on the WANT_ERRNO setting (which
> is on by default).

I've got a separate patch ready which sets WANT_ERRNO based on
_IEEE_LIBM so that both halves of the math library agree on whether
errno should be set. That's sufficient for my uses (where I always set
_IEEE_LIBM, and so never set errno) but not going to support the
old distinction between the __ieee754 functions and the posix functions.

I hadn't considered that issue. In my environment, I'm not enabling
errno support. I can think of a couple of possible fixes:

 1) Copy the __math_xflow functions to libm/math and remove errno
    handling so that the existing wrapper functions control how
    errno is managed.

 2) Remove the __ieee754 entry points from the library so that the
    old library only offers the regular API, and the _IEEE_LIBM
    configuration changes those functions to not modify errno. Linking
    the new math code so that it also used the same configuration value
    would keep everyone in sync and offer matching errno behavior.

Because the new math code uses double-precision arithmetic, I'm using
the old math code on devices without double-precision hardware, like ARM
processors that don't have __ARM_FP & 0x8. This includes pretty much all
of the processors I build products out of.

> in that case note that some compilers may model math functions as pure
> (no sideeffect) and apply optimizations accrodingly, this may create
> unexpected errno clobbers if math functions actually set errno.
> (-fno-math-errno behaviour in gcc) but this affects other libm
> implementations too and so far we haven't run into issues.

It would be kinda cool if we could drive all of the errno support based
on this compiler flag and eliminate the library configuration
setting. We'd want this to take effect while building the application
instead of at runtime.

> i think e.g. return __math_uflowf(s == -1.0f)
> would work better (but it's better to set a
> sign flag where s is set).

You're right -- here's an updated patch which uses (s < 0) as that
doesn't require an extra constant.

-- 
-keith

0001-libm-math-Use-__math_xflow-in-obsolete-math-code-v2.patch (text/x-diff, 10 KB)

From 259c26b2a6697a813d5cd41923079eb873321b3b Mon Sep 17 00:00:00 2001
From: Keith Packard <[email protected]>
Date: Thu, 30 Jul 2020 16:41:05 -0700
Subject: [PATCH] libm/math: Use __math_xflow in obsolete math code [v2]

C compilers may fold const values at compile time, so expressions
which try to elicit underflow/overflow by performing simple
arithemetic on suitable values will not generate the required
exceptions.

Work around this by replacing code which does these arithmetic
operations with calls to the existing __math_xflow functions that are
designed to do this correctly.

Signed-off-by: Keith Packard <[email protected]>

----

v2:
	libm/math: Pass sign to __math_xflow instead of muliplying result
---
 newlib/libm/common/math_errf.c |  2 +-
 newlib/libm/math/e_cosh.c      |  9 +++++----
 newlib/libm/math/e_exp.c       |  5 +++--
 newlib/libm/math/e_pow.c       | 18 ++++++++----------
 newlib/libm/math/ef_cosh.c     |  7 ++++---
 newlib/libm/math/ef_exp.c      |  5 +++--
 newlib/libm/math/ef_pow.c      | 14 ++++++--------
 newlib/libm/math/s_erf.c       |  3 ++-
 newlib/libm/math/sf_erf.c      |  3 ++-
 9 files changed, 34 insertions(+), 32 deletions(-)

diff --git a/newlib/libm/common/math_errf.c b/newlib/libm/common/math_errf.c
index 53c68b1cf..bb8273b8d 100644
--- a/newlib/libm/common/math_errf.c
+++ b/newlib/libm/common/math_errf.c
@@ -51,13 +51,13 @@ xflowf (uint32_t sign, float y)
   return with_errnof (y, ERANGE);
 }
 
-#if !__OBSOLETE_MATH
 HIDDEN float
 __math_uflowf (uint32_t sign)
 {
   return xflowf (sign, 0x1p-95f);
 }
 
+#if !__OBSOLETE_MATH
 #if WANT_ERRNO_UFLOW
 /* Underflows to zero in some non-nearest rounding mode, setting errno
    is valid even if the result is non-zero, but in the subnormal range.  */
diff --git a/newlib/libm/math/e_cosh.c b/newlib/libm/math/e_cosh.c
index a6310bd07..7b258ffef 100644
--- a/newlib/libm/math/e_cosh.c
+++ b/newlib/libm/math/e_cosh.c
@@ -25,7 +25,7 @@
  *			       			          2
  *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
  *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
- *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
+ *	    ln2ovft  <  x	    :  cosh(x) := overflow
  *
  * Special cases:
  *	cosh(x) is |x| if x is +INF, -INF, or NaN.
@@ -33,13 +33,14 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #ifndef _DOUBLE_IS_32BITS
 
 #ifdef __STDC__
-static const double one = 1.0, half=0.5, huge = 1.0e300;
+static const double one = 1.0, half=0.5;
 #else
-static double one = 1.0, half=0.5, huge = 1.0e300;
+static double one = 1.0, half=0.5;
 #endif
 
 #ifdef __STDC__
@@ -87,7 +88,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 	}
 
     /* |x| > overflowthresold, cosh(x) overflow */
-	return huge*huge;
+	return __math_oflow(0);
 }
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/e_exp.c b/newlib/libm/math/e_exp.c
index d23b1162b..ec26c2099 100644
--- a/newlib/libm/math/e_exp.c
+++ b/newlib/libm/math/e_exp.c
@@ -75,6 +75,7 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 #if __OBSOLETE_MATH
 
 #ifndef _DOUBLE_IS_32BITS
@@ -126,8 +127,8 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 		     return x+x; 		/* NaN */
 		else return (xsb==0)? x:0.0;	/* exp(+-inf)={inf,0} */
 	    }
-	    if(x > o_threshold) return huge*huge; /* overflow */
-	    if(x < u_threshold) return twom1000*twom1000; /* underflow */
+	    if(x > o_threshold) return __math_oflow(0); /* overflow */
+	    if(x < u_threshold) return __math_uflow(0); /* underflow */
 	}
 
     /* argument reduction */
diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index 4c450ec05..258cca8dd 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -76,8 +76,6 @@ zero    =  0.0,
 one	=  1.0,
 two	=  2.0,
 two53	=  9007199254740992.0,	/* 0x43400000, 0x00000000 */
-huge	=  1.0e300,
-tiny    =  1.0e-300,
 	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
 L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
 L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
@@ -197,12 +195,12 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
     /* |y| is huge */
 	if(iy>0x41e00000) { /* if |y| > 2**31 */
 	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow */
-		if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
-		if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
+		if(ix<=0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
+		if(ix>=0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
 	    }
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
-	    if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
+	    if(ix<0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
+	    if(ix>0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
 	/* now |1-x| is tiny <= 2**-20, suffice to compute 
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
 	    t = ax-1;		/* t has 20 trailing zeros */
@@ -275,15 +273,15 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	EXTRACT_WORDS(j,i,z);
 	if (j>=0x40900000) {				/* z >= 1024 */
 	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
-		return s*huge*huge;			/* overflow */
+		return __math_oflow(s<0);			/* overflow */
 	    else {
-		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
+		if(p_l+ovt>z-p_h) return __math_oflow(s<0);	/* overflow */
 	    }
 	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
 	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
-		return s*tiny*tiny;		/* underflow */
+		return __math_uflow(s<0);		/* underflow */
 	    else {
-		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
+		if(p_l<=z-p_h) return __math_uflow(s<0);	/* underflow */
 	    }
 	}
     /*
diff --git a/newlib/libm/math/ef_cosh.c b/newlib/libm/math/ef_cosh.c
index bdce61a00..5690bd7a4 100644
--- a/newlib/libm/math/ef_cosh.c
+++ b/newlib/libm/math/ef_cosh.c
@@ -14,15 +14,16 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #ifdef __v810__
 #define const
 #endif
 
 #ifdef __STDC__
-static const float one = 1.0, half=0.5, huge = 1.0e30;
+static const float one = 1.0, half=0.5;
 #else
-static float one = 1.0, half=0.5, huge = 1.0e30;
+static float one = 1.0, half=0.5;
 #endif
 
 #ifdef __STDC__
@@ -67,5 +68,5 @@ static float one = 1.0, half=0.5, huge = 1.0e30;
 	}
 
     /* |x| > overflowthresold, cosh(x) overflow */
-	return huge*huge;
+	return __math_oflowf(0);
 }
diff --git a/newlib/libm/math/ef_exp.c b/newlib/libm/math/ef_exp.c
index fb3e2ffe6..0cd0c00b3 100644
--- a/newlib/libm/math/ef_exp.c
+++ b/newlib/libm/math/ef_exp.c
@@ -14,6 +14,7 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #if __OBSOLETE_MATH
 #ifdef __v810__
@@ -61,9 +62,9 @@ P5   =  4.1381369442e-08; /* 0x3331bb4c */
         if(FLT_UWORD_IS_INFINITE(hx))
 	    return (xsb==0)? x:0.0;		/* exp(+-inf)={inf,0} */
 	if(sx > FLT_UWORD_LOG_MAX)
-	    return huge*huge; /* overflow */
+	    return __math_oflowf(0); /* overflow */
 	if(sx < 0 && hx > FLT_UWORD_LOG_MIN)
-	    return twom100*twom100; /* underflow */
+	    return __math_uflow(0); /* underflow */
 	
     /* argument reduction */
 	if(hx > 0x3eb17218) {		/* if  |x| > 0.5 ln2 */ 
diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
index d4ea4c5e8..e3579f071 100644
--- a/newlib/libm/math/ef_pow.c
+++ b/newlib/libm/math/ef_pow.c
@@ -33,8 +33,6 @@ zero    =  0.0,
 one	=  1.0,
 two	=  2.0,
 two24	=  16777216.0,	/* 0x4b800000 */
-huge	=  1.0e30,
-tiny    =  1.0e-30,
 	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
 L1  =  6.0000002384e-01, /* 0x3f19999a */
 L2  =  4.2857143283e-01, /* 0x3edb6db7 */
@@ -140,8 +138,8 @@ ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
     /* |y| is huge */
 	if(iy>0x4d000000) { /* if |y| > 2**27 */
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3f7ffff8) return (hy<0)? huge*huge:tiny*tiny;
-	    if(ix>0x3f800007) return (hy>0)? huge*huge:tiny*tiny;
+	    if(ix<0x3f7ffff8) return (hy<0)? __math_oflowf(0):__math_uflowf(0);
+	    if(ix>0x3f800007) return (hy>0)? __math_oflowf(0):__math_uflowf(0);
 	/* now |1-x| is tiny <= 2**-20, suffice to compute 
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
 	    t = ax-1;		/* t has 20 trailing zeros */
@@ -219,14 +217,14 @@ ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
 	i = j&0x7fffffff;
 	if (j>0) {
 	    if (i>FLT_UWORD_EXP_MAX)
-	        return s*huge*huge;			/* overflow */
+	        return __math_oflowf(s<0);			/* overflow */
 	    else if (i==FLT_UWORD_EXP_MAX)
-	        if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
+	        if(p_l+ovt>z-p_h) return __math_oflowf(s<0);	/* overflow */
         } else {
 	    if (i>FLT_UWORD_EXP_MIN)
-	        return s*tiny*tiny;			/* underflow */
+	        return __math_uflowf(s<0);			/* underflow */
 	    else if (i==FLT_UWORD_EXP_MIN)
-	        if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
+		if(p_l<=z-p_h) return __math_uflowf(s<0);	/* underflow */
 	}
     /*
      * compute 2**(p_h+p_l)
diff --git a/newlib/libm/math/s_erf.c b/newlib/libm/math/s_erf.c
index eb288fc73..9e2333c11 100644
--- a/newlib/libm/math/s_erf.c
+++ b/newlib/libm/math/s_erf.c
@@ -152,6 +152,7 @@ PORTABILITY
 
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #ifndef _DOUBLE_IS_32BITS
 
@@ -352,7 +353,7 @@ sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 			__ieee754_exp((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
-	    if(hx>0) return tiny*tiny; else return two-tiny;
+	    if(hx>0) return __math_uflow(0); else return two-tiny;
 	}
 }
 
diff --git a/newlib/libm/math/sf_erf.c b/newlib/libm/math/sf_erf.c
index 0329c60fa..f3d0de97a 100644
--- a/newlib/libm/math/sf_erf.c
+++ b/newlib/libm/math/sf_erf.c
@@ -14,6 +14,7 @@
  */
 
 #include "fdlibm.h"
+#include "math_config.h"
 
 #ifdef __v810__
 #define const
@@ -217,7 +218,7 @@ sb7  = -2.2440952301e+01; /* 0xc1b38712 */
 			__ieee754_expf((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
-	    if(hx>0) return tiny*tiny; else return two-tiny;
+	    if(hx>0) return __math_uflow(0); else return two-tiny;
 	}
 }
 
-- 
2.28.0

signature.asc (application/pgp-signature, 832 B)