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

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

[Keith Packard]

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 <keithp@keithp.com>
---
 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..e68a76e06 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 s*__math_oflow(0);			/* overflow */
 	    else {
-		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
+		if(p_l+ovt>z-p_h) return s*__math_oflow(0);	/* overflow */
 	    }
 	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
 	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
-		return s*tiny*tiny;		/* underflow */
+		return s*__math_uflow(0);		/* underflow */
 	    else {
-		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
+		if(p_l<=z-p_h) return s*__math_uflow(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..3f5c5c7e1 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 s*__math_oflowf(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 s*__math_oflowf(0);	/* overflow */
         } else {
 	    if (i>FLT_UWORD_EXP_MIN)
-	        return s*tiny*tiny;			/* underflow */
+	        return s*__math_uflowf(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 s*__math_uflowf(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.rc2

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

[Szabolcs Nagy]

The 07/30/2020 16:42, Keith Packard via Newlib wrote:
> 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 <keithp@keithp.com>

this looks good.

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.

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

in the old style math functions errno is
handled in wrappers but presumably this fix
is for the case when error handling wrappers
are not in use.

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.

> --- a/newlib/libm/math/ef_pow.c
> +++ b/newlib/libm/math/ef_pow.c
...
> @@ -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 s*__math_oflowf(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 s*__math_oflowf(0);	/* overflow */
>          } else {
>  	    if (i>FLT_UWORD_EXP_MIN)
> -	        return s*tiny*tiny;			/* underflow */
> +	        return s*__math_uflowf(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 s*__math_uflowf(0);	/* underflow */

this is suboptimal: the reason there is
a 'sign' argument to these functions is
that you don't need mul with the sign
(which is wrong with upward or downward
rounding:

 s*huge*huge != (s*huge)*huge

but newlib may not care about rounding
modes, however it also prevents the
call to be a tail call and turns the
math function into a non-leaf function)

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).

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

[Keith Packard]

[-- Attachment #1.1: Type: text/plain, Size: 2453 bytes --]

Szabolcs Nagy <szabolcs.nagy@arm.com> 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.


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1.2: 0001-libm-math-Use-__math_xflow-in-obsolete-math-code-v2.patch --]
[-- Type: text/x-diff, Size: 10538 bytes --]

From 259c26b2a6697a813d5cd41923079eb873321b3b Mon Sep 17 00:00:00 2001
From: Keith Packard <keithp@keithp.com>
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 <keithp@keithp.com>

----

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


[-- Attachment #1.3: Type: text/plain, Size: 15 bytes --]


-- 
-keith

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

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

[Szabolcs Nagy]

The 08/01/2020 15:40, Keith Packard wrote:
> Szabolcs Nagy <szabolcs.nagy@arm.com> writes:
> > 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.

c99 code is fine with me, but if newlib
decides it requires a c99 compiler for
building then it might consider cleaning
up the pre-c89 compatibility code it has.

> > 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.

sounds good.

i think math errno is not particularly useful in newlib.

> 
> 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.

different translation units may be compiled with
different flags, so separate symbols may be needed.

but i think for now a newlib build time setting is ok.


> > 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.
> 

v2 looks good.


> From 259c26b2a6697a813d5cd41923079eb873321b3b Mon Sep 17 00:00:00 2001
> From: Keith Packard <keithp@keithp.com>
> 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 <keithp@keithp.com>
> 
> ----
> 
> 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
> 

> 
> -- 
> -keith




-- 

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

[Corinna Vinschen]

On Aug  3 10:21, Szabolcs Nagy wrote:
> The 08/01/2020 15:40, Keith Packard wrote:
> > You're right -- here's an updated patch which uses (s < 0) as that
> > doesn't require an extra constant.
> > 
> 
> v2 looks good.
> 
> 
> > From 259c26b2a6697a813d5cd41923079eb873321b3b Mon Sep 17 00:00:00 2001
> > From: Keith Packard <keithp@keithp.com>
> > 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 <keithp@keithp.com>
> > 
> > ----
> > 
> > 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(-)

Pushed.


Thanks,
Corinna

-- 
Corinna Vinschen
Cygwin Maintainer
Red Hat

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

[Keith Packard]

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

Szabolcs Nagy <szabolcs.nagy@arm.com> writes:

> i think math errno is not particularly useful in newlib.

For embedded systems, I happen to agree, which is why picolibc builds
without errno support by default.

Thanks much for review and helpful suggestions!

-- 
-keith

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

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

[Keith Packard]

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

Szabolcs Nagy <szabolcs.nagy@arm.com> writes:

> 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.

Hrm. I just learned about "#pragma STDC FENV_ACCESS on" this morning and
am wondering if we shouldn't just patch the old math library to use that
instead of calling the __math_xflow functions. Would be a simpler patch
and provide precise compatibility with the previous version.

Here's an example of how that works:

        #include <fenv.h>
        #include <stdio.h>

        #pragma STDC FENV_ACCESS on

        int main(void)
        {
                double x;
                int e;
                feclearexcept(FE_ALL_EXCEPT);
                x = 0.0 / 0.0;
                e = fetestexcept(FE_ALL_EXCEPT);
                printf("x %g e %x\n", x, e);
                return 0;
        }

$ cc except.c -lm
$ ./a.out
x -nan e 1
$

As you can see, even using constants still raise the appropriate
exception.

-- 
-keith

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

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

[Szabolcs Nagy]

The 08/03/2020 11:07, Keith Packard wrote:
> Szabolcs Nagy <szabolcs.nagy@arm.com> writes:
> 
> > 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.
> 
> Hrm. I just learned about "#pragma STDC FENV_ACCESS on" this morning and
> am wondering if we shouldn't just patch the old math library to use that
> instead of calling the __math_xflow functions. Would be a simpler patch
> and provide precise compatibility with the previous version.

i don't know of any c compiler that correctly
implements that pragma and some compilers
(clang) warn about it so i don't use it myself.

but the correct iso c way is to use it on every
code that may run in non-default fenv (non-nearest
rounding mode) or may run between exception
setting and getting operations. (i.e. not just for
math code that has fenv functions in it).

the drawback is that it means that any call may
access the fenv (i.e. if you mark a function
implementation with the pragma then any call
within that function may change the rounding mode
unless the compiler can prove otherwise), so all
extern calls become floating-point barriers which
prevents some useful optimizations (most calls
don't change the rounding mode). and this still
does not give complete fenv safety, because it's
not clear if inline asm or other language
extensions may access fenv or not and how
compilers will deal with that: currently float
operations are reordered across inline asm.

> 
> Here's an example of how that works:
> 
>         #include <fenv.h>
>         #include <stdio.h>
> 
>         #pragma STDC FENV_ACCESS on
> 
>         int main(void)
>         {
>                 double x;
>                 int e;
>                 feclearexcept(FE_ALL_EXCEPT);
>                 x = 0.0 / 0.0;
>                 e = fetestexcept(FE_ALL_EXCEPT);
>                 printf("x %g e %x\n", x, e);
>                 return 0;
>         }
> 
> $ cc except.c -lm
> $ ./a.out
> x -nan e 1
> $
> 
> As you can see, even using constants still raise the appropriate
> exception.

either you got lucky or you use another compiler
than i do.

(gcc gets close to exception safe operation if
you use -frounding-math: then it assumes the
compiled code itself does not change the rounding
mode, but its caller might, which e.g. means that
constant folding is not allowed: inexact arithmetic
operations are rounding mode dependent. so the
generated code mostly preserves fenv exceptions,
but this can break if the result does not depend
on rounding: nextafter(DBL_MIN,0) should raise
underflow but it may be const folded, or
relational operations may not preserve exceptions)

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

[Keith Packard]

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

Szabolcs Nagy <szabolcs.nagy@arm.com> writes:

> either you got lucky or you use another compiler
> than i do.

Yeah. Just lucky. Thanks for the warning. My fenv tests in picolibc are
using 'volatile' to force the desired order-of-operations.

        https://github.com/keith-packard/picolibc/blob/main/test/fenv.c

This gets tested in a huge range of configurations on RISC-V, x86 and
ARM and generates the expected results in all of them.

-- 
-keith

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