[PATCH] libm/math: Use __math_xflow in obsolete math code
Keith Packard
keithp@keithp.com
Thu Jul 30 23:42:36 GMT 2020
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
More information about the Newlib
mailing list