drew
/
acadia-newlib
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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
This commit is contained in:
Keith Packard 2020-07-30 16:41:05 -07:00 committed by Corinna Vinschen
parent 5717262b8e
commit 12ad9a46df
9 changed files with 34 additions and 32 deletions

2
newlib/libm/common/math_errf.c
View File

@ -51,13 +51,13 @@ xflowf (uint32_t sign, float y)
return with_errnof (y, ERANGE); return with_errnof (y, ERANGE);
} }
#if !__OBSOLETE_MATH
HIDDEN float HIDDEN float
__math_uflowf (uint32_t sign) __math_uflowf (uint32_t sign)
{ {
return xflowf (sign, 0x1p-95f); return xflowf (sign, 0x1p-95f);
} }
#if !__OBSOLETE_MATH
#if WANT_ERRNO_UFLOW #if WANT_ERRNO_UFLOW
/* Underflows to zero in some non-nearest rounding mode, setting errno /* 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. */ is valid even if the result is non-zero, but in the subnormal range. */

9
newlib/libm/math/e_cosh.c
View File

@ -25,7 +25,7 @@
* 2 * 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2 * 22 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * 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: * Special cases:
* cosh(x) is |x| if x is +INF, -INF, or NaN. * cosh(x) is |x| if x is +INF, -INF, or NaN.
@ -33,13 +33,14 @@
*/ */
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#ifndef _DOUBLE_IS_32BITS #ifndef _DOUBLE_IS_32BITS
#ifdef __STDC__ #ifdef __STDC__
static const double one = 1.0, half=0.5, huge = 1.0e300; static const double one = 1.0, half=0.5;
#else #else
static double one = 1.0, half=0.5, huge = 1.0e300; static double one = 1.0, half=0.5;
#endif #endif
#ifdef __STDC__ #ifdef __STDC__
@ -87,7 +88,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
} }
/* |x| > overflowthresold, cosh(x) overflow */ /* |x| > overflowthresold, cosh(x) overflow */
return huge*huge; return __math_oflow(0);
} }
#endif /* defined(_DOUBLE_IS_32BITS) */ #endif /* defined(_DOUBLE_IS_32BITS) */

5
newlib/libm/math/e_exp.c
View File

@ -75,6 +75,7 @@
*/ */
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#if __OBSOLETE_MATH #if __OBSOLETE_MATH
#ifndef _DOUBLE_IS_32BITS #ifndef _DOUBLE_IS_32BITS
@ -126,8 +127,8 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
return x+x; /* NaN */ return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
} }
if(x > o_threshold) return huge*huge; /* overflow */ if(x > o_threshold) return __math_oflow(0); /* overflow */
if(x < u_threshold) return twom1000*twom1000; /* underflow */ if(x < u_threshold) return __math_uflow(0); /* underflow */
} }
/* argument reduction */ /* argument reduction */

18
newlib/libm/math/e_pow.c
View File

@ -76,8 +76,6 @@ zero = 0.0,
one = 1.0, one = 1.0,
two = 2.0, two = 2.0,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
huge = 1.0e300,
tiny = 1.0e-300,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
@ -197,12 +195,12 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
/* |y| is huge */ /* |y| is huge */
if(iy>0x41e00000) { /* if |y| > 2**31 */ if(iy>0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix<=0x3fefffff) 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)? huge*huge:tiny*tiny; if(ix>=0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
} }
/* over/underflow if x is not close to one */ /* over/underflow if x is not close to one */
if(ix<0x3fefffff) 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)? huge*huge:tiny*tiny; if(ix>0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
/* now |1-x| is tiny <= 2**-20, suffice to compute /* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */ log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-1; /* t has 20 trailing zeros */ 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); EXTRACT_WORDS(j,i,z);
if (j>=0x40900000) { /* z >= 1024 */ if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */ if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */ return __math_oflow(s<0); /* overflow */
else { 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 */ } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
return s*tiny*tiny; /* underflow */ return __math_uflow(s<0); /* underflow */
else { else {
if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ if(p_l<=z-p_h) return __math_uflow(s<0); /* underflow */
} }
} }
/* /*

7
newlib/libm/math/ef_cosh.c
View File

@ -14,15 +14,16 @@
*/ */
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#ifdef __v810__ #ifdef __v810__
#define const #define const
#endif #endif
#ifdef __STDC__ #ifdef __STDC__
static const float one = 1.0, half=0.5, huge = 1.0e30; static const float one = 1.0, half=0.5;
#else #else
static float one = 1.0, half=0.5, huge = 1.0e30; static float one = 1.0, half=0.5;
#endif #endif
#ifdef __STDC__ #ifdef __STDC__
@ -67,5 +68,5 @@ static float one = 1.0, half=0.5, huge = 1.0e30;
} }
/* |x| > overflowthresold, cosh(x) overflow */ /* |x| > overflowthresold, cosh(x) overflow */
return huge*huge; return __math_oflowf(0);
} }

5
newlib/libm/math/ef_exp.c
View File

@ -14,6 +14,7 @@
*/ */
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#if __OBSOLETE_MATH #if __OBSOLETE_MATH
#ifdef __v810__ #ifdef __v810__
@ -61,9 +62,9 @@ P5 = 4.1381369442e-08; /* 0x3331bb4c */
if(FLT_UWORD_IS_INFINITE(hx)) if(FLT_UWORD_IS_INFINITE(hx))
return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
if(sx > FLT_UWORD_LOG_MAX) if(sx > FLT_UWORD_LOG_MAX)
return huge*huge; /* overflow */ return __math_oflowf(0); /* overflow */
if(sx < 0 && hx > FLT_UWORD_LOG_MIN) if(sx < 0 && hx > FLT_UWORD_LOG_MIN)
return twom100*twom100; /* underflow */ return __math_uflow(0); /* underflow */
/* argument reduction */ /* argument reduction */
if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */

14
newlib/libm/math/ef_pow.c
View File

@ -33,8 +33,6 @@ zero = 0.0,
one = 1.0, one = 1.0,
two = 2.0, two = 2.0,
two24 = 16777216.0, /* 0x4b800000 */ two24 = 16777216.0, /* 0x4b800000 */
huge = 1.0e30,
tiny = 1.0e-30,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 6.0000002384e-01, /* 0x3f19999a */ L1 = 6.0000002384e-01, /* 0x3f19999a */
L2 = 4.2857143283e-01, /* 0x3edb6db7 */ L2 = 4.2857143283e-01, /* 0x3edb6db7 */
@ -140,8 +138,8 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
/* |y| is huge */ /* |y| is huge */
if(iy>0x4d000000) { /* if |y| > 2**27 */ if(iy>0x4d000000) { /* if |y| > 2**27 */
/* over/underflow if x is not close to one */ /* over/underflow if x is not close to one */
if(ix<0x3f7ffff8) 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)? huge*huge:tiny*tiny; if(ix>0x3f800007) return (hy>0)? __math_oflowf(0):__math_uflowf(0);
/* now |1-x| is tiny <= 2**-20, suffice to compute /* now |1-x| is tiny <= 2**-20, suffice to compute
log(x) by x-x^2/2+x^3/3-x^4/4 */ log(x) by x-x^2/2+x^3/3-x^4/4 */
t = ax-1; /* t has 20 trailing zeros */ t = ax-1; /* t has 20 trailing zeros */
@ -219,14 +217,14 @@ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
i = j&0x7fffffff; i = j&0x7fffffff;
if (j>0) { if (j>0) {
if (i>FLT_UWORD_EXP_MAX) if (i>FLT_UWORD_EXP_MAX)
return s*huge*huge; /* overflow */ return __math_oflowf(s<0); /* overflow */
else if (i==FLT_UWORD_EXP_MAX) 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 { } else {
if (i>FLT_UWORD_EXP_MIN) if (i>FLT_UWORD_EXP_MIN)
return s*tiny*tiny; /* underflow */ return __math_uflowf(s<0); /* underflow */
else if (i==FLT_UWORD_EXP_MIN) 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) * compute 2**(p_h+p_l)

3
newlib/libm/math/s_erf.c
View File

@ -152,6 +152,7 @@ PORTABILITY
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#ifndef _DOUBLE_IS_32BITS #ifndef _DOUBLE_IS_32BITS
@ -352,7 +353,7 @@ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
__ieee754_exp((z-x)*(z+x)+R/S); __ieee754_exp((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x; if(hx>0) return r/x; else return two-r/x;
} else { } else {
if(hx>0) return tiny*tiny; else return two-tiny; if(hx>0) return __math_uflow(0); else return two-tiny;
} }
} }

3
newlib/libm/math/sf_erf.c
View File

@ -14,6 +14,7 @@
*/ */
#include "fdlibm.h" #include "fdlibm.h"
#include "math_config.h"
#ifdef __v810__ #ifdef __v810__
#define const #define const
@ -217,7 +218,7 @@ sb7 = -2.2440952301e+01; /* 0xc1b38712 */
__ieee754_expf((z-x)*(z+x)+R/S); __ieee754_expf((z-x)*(z+x)+R/S);
if(hx>0) return r/x; else return two-r/x; if(hx>0) return r/x; else return two-r/x;
} else { } else {
if(hx>0) return tiny*tiny; else return two-tiny; if(hx>0) return __math_uflow(0); else return two-tiny;
} }
} }