diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/e_logl.c')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_logl.c | 53 |
1 files changed, 34 insertions, 19 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_logl.c b/sysdeps/ieee754/ldbl-128ibm/e_logl.c index 14f47ebade..58d6bc6972 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_logl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_logl.c @@ -182,21 +182,26 @@ static const long double ln2a = 6.93145751953125e-1L, ln2b = 1.4286068203094172321214581765680755001344E-6L; +static const long double + ldbl_epsilon = 0x1p-106L; + long double __ieee754_logl(long double x) { - long double z, y, w; - ieee854_long_double_shape_type u, t; + long double z, y, w, t; unsigned int m; int k, e; + double xhi; + uint32_t hx, lx; - u.value = x; - m = u.parts32.w0; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + m = hx; /* Check for IEEE special cases. */ k = m & 0x7fffffff; /* log(0) = -infinity. */ - if ((k | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0) + if ((k | lx) == 0) { return -0.5L / ZERO; } @@ -216,7 +221,7 @@ __ieee754_logl(long double x) { z = x - 1.0L; k = 64; - t.value = 1.0L; + t = 1.0L; e = 0; } else @@ -224,6 +229,14 @@ __ieee754_logl(long double x) /* Extract exponent and reduce domain to 0.703125 <= u < 1.40625 */ unsigned int w0; e = (int) (m >> 20) - (int) 0x3fe; + if (e == -1022) + { + x *= 0x1p106L; + xhi = ldbl_high (x); + EXTRACT_WORDS (hx, lx, xhi); + m = hx; + e = (int) (m >> 20) - (int) 0x3fe - 106; + } m &= 0xfffff; w0 = m | 0x3fe00000; m |= 0x100000; @@ -233,10 +246,8 @@ __ieee754_logl(long double x) k = (m - 0xff000) >> 13; /* t is the argument 0.5 + (k+26)/128 of the nearest item to u in the lookup table. */ - t.parts32.w0 = 0x3ff00000 + (k << 13); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3ff00000 + (k << 13), 0); + t = xhi; w0 += 0x100000; e -= 1; k += 64; @@ -244,21 +255,24 @@ __ieee754_logl(long double x) else { k = (m - 0xfe000) >> 14; - t.parts32.w0 = 0x3fe00000 + (k << 14); - t.parts32.w1 = 0; - t.parts32.w2 = 0; - t.parts32.w3 = 0; + INSERT_WORDS (xhi, 0x3fe00000 + (k << 14), 0); + t = xhi; } - u.value = __scalbnl (u.value, ((int) ((w0 - u.parts32.w0) * 2)) >> 21); + x = __scalbnl (x, ((int) ((w0 - hx) * 2)) >> 21); /* log(u) = log( t u/t ) = log(t) + log(u/t) log(t) is tabulated in the lookup table. Express log(u/t) = log(1+z), where z = u/t - 1 = (u-t)/t. cf. Cody & Waite. */ - z = (u.value - t.value) / t.value; + z = (x - t) / t; } /* Series expansion of log(1+z). */ w = z * z; - y = ((((((((((((l15 * z + /* Avoid spurious underflows. */ + if (__glibc_unlikely(w <= ldbl_epsilon)) + y = 0.0L; + else + { + y = ((((((((((((l15 * z + l14) * z + l13) * z + l12) * z @@ -271,11 +285,12 @@ __ieee754_logl(long double x) + l5) * z + l4) * z + l3) * z * w; - y -= 0.5 * w; + y -= 0.5 * w; + } y += e * ln2b; /* Base 2 exponent offset times ln(2). */ y += z; y += logtbl[k-26]; /* log(t) - (t-1) */ - y += (t.value - 1.0L); + y += (t - 1.0L); y += e * ln2a; return y; } |