diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_log1p.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_log1p.c | 23 |
1 files changed, 13 insertions, 10 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index ea1dc6ca76..cff555b0aa 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -78,6 +78,7 @@ * See HP-15C Advanced Functions Handbook, p.193. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -107,18 +108,25 @@ __log1p (double x) k = 1; if (hx < 0x3FDA827A) /* x < 0.41422 */ { - if (__builtin_expect (ax >= 0x3ff00000, 0)) /* x <= -1.0 */ + if (__glibc_unlikely (ax >= 0x3ff00000)) /* x <= -1.0 */ { if (x == -1.0) - return -two54 / (x - x); /* log1p(-1)=+inf */ + return -two54 / zero; /* log1p(-1)=-inf */ else return (x - x) / (x - x); /* log1p(x<-1)=NaN */ } - if (__builtin_expect (ax < 0x3e200000, 0)) /* |x| < 2**-29 */ + if (__glibc_unlikely (ax < 0x3e200000)) /* |x| < 2**-29 */ { math_force_eval (two54 + x); /* raise inexact */ if (ax < 0x3c900000) /* |x| < 2**-54 */ - return x; + { + if (fabs (x) < DBL_MIN) + { + double force_underflow = x * x; + math_force_eval (force_underflow); + } + return x; + } else return x - x * x * 0.5; } @@ -127,7 +135,7 @@ __log1p (double x) k = 0; f = x; hu = 1; } /* -0.2929<x<0.41422 */ } - else if (__builtin_expect (hx >= 0x7ff00000, 0)) + else if (__glibc_unlikely (hx >= 0x7ff00000)) return x + x; if (k != 0) { @@ -189,8 +197,3 @@ __log1p (double x) else return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); } -weak_alias (__log1p, log1p) -#ifdef NO_LONG_DOUBLE -strong_alias (__log1p, __log1pl) -weak_alias (__log1p, log1pl) -#endif |