aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--math/exp.c7
-rw-r--r--math/exp2.c5
-rw-r--r--math/log.c1
-rw-r--r--math/log2.c1
-rw-r--r--math/math_config.h26
-rw-r--r--math/pow.c11
6 files changed, 50 insertions, 1 deletions
diff --git a/math/exp.c b/math/exp.c
index e73f2a6..ce47b01 100644
--- a/math/exp.c
+++ b/math/exp.c
@@ -33,6 +33,10 @@
#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
+/* Handle inputs that may overflow or underflow when computing the result
+ that is scale*(1+tmp), the exponent bits of scale might have overflown
+ into the sign bit so that needs correction before sbits is used as a
+ double, ki is only used to determine the sign of the exponent. */
static inline double
specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
{
@@ -71,12 +75,15 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
return check_uflow (y);
}
+/* Top 12 bits of a double (sign and exponent bits). */
static inline uint32_t
top12 (double x)
{
return asuint64 (x) >> 52;
}
+/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+ If hastail is 0 then xtail is assumed to be 0 too. */
static inline double
exp_inline (double x, double xtail, int hastail)
{
diff --git a/math/exp2.c b/math/exp2.c
index 3d8aadf..a64a26c 100644
--- a/math/exp2.c
+++ b/math/exp2.c
@@ -31,6 +31,10 @@
#define C5 __exp_data.exp2_poly[4]
#define C6 __exp_data.exp2_poly[5]
+/* Handle inputs that may overflow or underflow when computing the result
+ that is scale*(1+tmp), the exponent bits of scale might have overflown
+ into the sign bit so that needs correction before sbits is used as a
+ double, ki is only used to determine the sign of the exponent. */
static inline double
specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
{
@@ -69,6 +73,7 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
return check_uflow (y);
}
+/* Top 12 bits of a double (sign and exponent bits). */
static inline uint32_t
top12 (double x)
{
diff --git a/math/log.c b/math/log.c
index 324936d..29e0cc1 100644
--- a/math/log.c
+++ b/math/log.c
@@ -30,6 +30,7 @@
#define N (1 << LOG_TABLE_BITS)
#define OFF 0x3fe6000000000000
+/* Top 16 bits of a double. */
static inline uint32_t
top16 (double x)
{
diff --git a/math/log2.c b/math/log2.c
index 2b3960d..60c2e44 100644
--- a/math/log2.c
+++ b/math/log2.c
@@ -30,6 +30,7 @@
#define N (1 << LOG2_TABLE_BITS)
#define OFF 0x3fe6000000000000
+/* Top 16 bits of a double. */
static inline uint32_t
top16 (double x)
{
diff --git a/math/math_config.h b/math/math_config.h
index 3986d48..ee44468 100644
--- a/math/math_config.h
+++ b/math/math_config.h
@@ -66,12 +66,16 @@
#if HAVE_FAST_ROUND
# define TOINT_INTRINSICS 1
+/* Round x to nearest int, ties have to be rounded consistently with
+ converttoint. The result is in [-2^31+1, 2^31-1]. */
static inline double_t
roundtoint (double_t x)
{
return round (x);
}
+/* Convert x to nearest int, ties have to be rounded consistently with
+ roundtoint. The result is in [-2^31+1, 2^31-1]. */
static inline uint64_t
converttoint (double_t x)
{
@@ -240,28 +244,48 @@ eval_as_double (double x)
# define unlikely(x) (x)
#endif
-/* Error handling tail calls for special cases, with sign argument. */
+/* Error handling tail calls for special cases, with a sign argument.
+ The sign of the return value is set if the argument is non-zero. */
+
+/* The result overflows. */
HIDDEN float __math_oflowf (uint32_t);
+/* The result underflows to 0 in nearest rounding mode. */
HIDDEN float __math_uflowf (uint32_t);
+/* The result underflows to 0 in some directed rounding mode only. */
HIDDEN float __math_may_uflowf (uint32_t);
+/* Division by zero. */
HIDDEN float __math_divzerof (uint32_t);
+/* The result overflows. */
HIDDEN double __math_oflow (uint32_t);
+/* The result underflows to 0 in nearest rounding mode. */
HIDDEN double __math_uflow (uint32_t);
+/* The result underflows to 0 in some directed rounding mode only. */
HIDDEN double __math_may_uflow (uint32_t);
+/* Division by zero. */
HIDDEN double __math_divzero (uint32_t);
+
/* Error handling using input checking. */
+
+/* Invalid input unless it is a quiet NaN. */
HIDDEN float __math_invalidf (float);
+/* Invalid input unless it is a quiet NaN. */
HIDDEN double __math_invalid (double);
+
/* Error handling using output checking, only for errno setting. */
+
+/* Check if the result overflowed to infinity. */
HIDDEN double __math_check_oflow (double);
+/* Check if the result underflowed to 0. */
HIDDEN double __math_check_uflow (double);
+/* Check if the result overflowed to infinity. */
static inline double
check_oflow (double x)
{
return WANT_ERRNO ? __math_check_oflow (x) : x;
}
+/* Check if the result underflowed to 0. */
static inline double
check_uflow (double x)
{
diff --git a/math/pow.c b/math/pow.c
index 90f2976..d391509 100644
--- a/math/pow.c
+++ b/math/pow.c
@@ -34,12 +34,16 @@ ulperr_exp: 0.509 ULP (ULP error of exp, 0.511 ULP without fma)
#define N (1 << POW_LOG_TABLE_BITS)
#define OFF 0x3fe6955500000000
+/* Top 12 bits of a double (sign and exponent bits). */
static inline uint32_t
top12 (double x)
{
return asuint64 (x) >> 52;
}
+/* Compute y+tail = log(x) where the rounded result is y and tail has about
+ additional 15 bits precision. The bit representation of x if in ix, but
+ normalized in the subnormal range using sign bit too for the exponent. */
static inline double_t
log_inline (uint64_t ix, double_t *tail)
{
@@ -122,6 +126,10 @@ log_inline (uint64_t ix, double_t *tail)
#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
#define C6 __exp_data.poly[9 - EXP_POLY_ORDER]
+/* Handle inputs that may overflow or underflow when computing the result
+ that is scale*(1+tmp), the exponent bits of scale might have overflown
+ into the sign bit so that needs correction before sbits is used as a
+ double, ki is only used to determine the sign of the exponent. */
static inline double
specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
{
@@ -165,6 +173,8 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
#define SIGN_BIAS (0x800 << EXP_TABLE_BITS)
+/* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+ The sign_bias argument is SIGN_BIAS or 0 and sets the sign to -1 or 1. */
static inline double
exp_inline (double x, double xtail, uint32_t sign_bias)
{
@@ -257,6 +267,7 @@ checkint (uint64_t iy)
return 2;
}
+/* Returns 1 if input is the bit representation of 0, infinity or nan. */
static inline int
zeroinfnan (uint64_t i)
{