* [PATCH] libm/math: Use __math_xflow in obsolete math code @ 2020-07-30 23:42 Keith Packard 2020-07-31 15:19 ` Szabolcs Nagy 0 siblings, 1 reply; 9+ messages in thread From: Keith Packard @ 2020-07-30 23:42 UTC (permalink / raw) To: newlib 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 ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-07-30 23:42 [PATCH] libm/math: Use __math_xflow in obsolete math code Keith Packard @ 2020-07-31 15:19 ` Szabolcs Nagy 2020-08-01 22:40 ` Keith Packard 2020-08-03 18:07 ` Keith Packard 0 siblings, 2 replies; 9+ messages in thread From: Szabolcs Nagy @ 2020-07-31 15:19 UTC (permalink / raw) To: Keith Packard; +Cc: newlib 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). ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-07-31 15:19 ` Szabolcs Nagy @ 2020-08-01 22:40 ` Keith Packard 2020-08-03 9:21 ` Szabolcs Nagy 2020-08-03 18:07 ` Keith Packard 1 sibling, 1 reply; 9+ messages in thread From: Keith Packard @ 2020-08-01 22:40 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: newlib [-- 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 --] ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-08-01 22:40 ` Keith Packard @ 2020-08-03 9:21 ` Szabolcs Nagy 2020-08-03 11:32 ` Corinna Vinschen 2020-08-03 17:56 ` Keith Packard 0 siblings, 2 replies; 9+ messages in thread From: Szabolcs Nagy @ 2020-08-03 9:21 UTC (permalink / raw) To: Keith Packard; +Cc: newlib 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 -- ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-08-03 9:21 ` Szabolcs Nagy @ 2020-08-03 11:32 ` Corinna Vinschen 2020-08-03 17:56 ` Keith Packard 1 sibling, 0 replies; 9+ messages in thread From: Corinna Vinschen @ 2020-08-03 11:32 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: Keith Packard, newlib 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 ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-08-03 9:21 ` Szabolcs Nagy 2020-08-03 11:32 ` Corinna Vinschen @ 2020-08-03 17:56 ` Keith Packard 1 sibling, 0 replies; 9+ messages in thread From: Keith Packard @ 2020-08-03 17:56 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: newlib [-- 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 --] ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-07-31 15:19 ` Szabolcs Nagy 2020-08-01 22:40 ` Keith Packard @ 2020-08-03 18:07 ` Keith Packard 2020-08-03 21:35 ` Szabolcs Nagy 1 sibling, 1 reply; 9+ messages in thread From: Keith Packard @ 2020-08-03 18:07 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: newlib [-- 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 --] ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-08-03 18:07 ` Keith Packard @ 2020-08-03 21:35 ` Szabolcs Nagy 2020-08-03 23:21 ` Keith Packard 0 siblings, 1 reply; 9+ messages in thread From: Szabolcs Nagy @ 2020-08-03 21:35 UTC (permalink / raw) To: Keith Packard; +Cc: newlib 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) ^ permalink raw reply [flat|nested] 9+ messages in thread
* Re: [PATCH] libm/math: Use __math_xflow in obsolete math code 2020-08-03 21:35 ` Szabolcs Nagy @ 2020-08-03 23:21 ` Keith Packard 0 siblings, 0 replies; 9+ messages in thread From: Keith Packard @ 2020-08-03 23:21 UTC (permalink / raw) To: Szabolcs Nagy; +Cc: newlib [-- 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 --] ^ permalink raw reply [flat|nested] 9+ messages in thread
end of thread, other threads:[~2020-08-03 23:21 UTC | newest] Thread overview: 9+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2020-07-30 23:42 [PATCH] libm/math: Use __math_xflow in obsolete math code Keith Packard 2020-07-31 15:19 ` Szabolcs Nagy 2020-08-01 22:40 ` Keith Packard 2020-08-03 9:21 ` Szabolcs Nagy 2020-08-03 11:32 ` Corinna Vinschen 2020-08-03 17:56 ` Keith Packard 2020-08-03 18:07 ` Keith Packard 2020-08-03 21:35 ` Szabolcs Nagy 2020-08-03 23:21 ` Keith Packard
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).