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

The 08/01/2020 15:40, Keith Packard wrote:
> Szabolcs Nagy <[email protected]> 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 <[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
> 

> 
> -- 
> -keith




--