[PATCH] libm/math: Use __math_xflow in obsolete math code
Szabolcs Nagy
szabolcs.nagy@arm.com
Mon Aug 3 09:21:08 GMT 2020
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
--
More information about the Newlib
mailing list