diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm')
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/s_frexpl.c | 146 |
1 files changed, 102 insertions, 44 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c index 7e40663..a644f92 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c @@ -31,57 +31,115 @@ static char rcsid[] = "$NetBSD: $"; #include <math_private.h> #include <math_ldbl_opt.h> -static const long double -two107 = 162259276829213363391578010288128.0; /* 0x4670000000000000, 0 */ - long double __frexpl(long double x, int *eptr) { - uint64_t hx, lx, ix, ixl; - int64_t explo; - double xhi, xlo; + uint64_t hx, lx, ix, ixl; + int64_t explo, expon; + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ixl = 0x7fffffffffffffffULL & lx; + ix = 0x7fffffffffffffffULL & hx; + expon = 0; + if (ix >= 0x7ff0000000000000ULL || ix == 0) + { + /* 0,inf,nan. */ + *eptr = expon; + return x; + } + expon = ix >> 52; + if (expon == 0) + { + /* Denormal high double, the low double must be 0.0. */ + int cnt; + + /* Normalize. */ + if (sizeof (ix) == sizeof (long)) + cnt = __builtin_clzl (ix); + else if ((ix >> 32) != 0) + cnt = __builtin_clzl ((long) (ix >> 32)); + else + cnt = __builtin_clzl ((long) ix) + 32; + cnt = cnt - 12; + expon -= cnt; + ix <<= cnt + 1; + } + expon -= 1022; + ix &= 0x000fffffffffffffULL; + hx &= 0x8000000000000000ULL; + hx |= (1022LL << 52) | ix; - ldbl_unpack (x, &xhi, &xlo); - EXTRACT_WORDS64 (hx, xhi); - EXTRACT_WORDS64 (lx, xlo); - ixl = 0x7fffffffffffffffULL&lx; - ix = 0x7fffffffffffffffULL&hx; - *eptr = 0; - if(ix>=0x7ff0000000000000ULL||ix==0) return x; /* 0,inf,nan */ - if (ix<0x0010000000000000ULL) { /* subnormal */ - x *= two107; - xhi = ldbl_high (x); - EXTRACT_WORDS64 (hx, xhi); - ix = hx&0x7fffffffffffffffULL; - *eptr = -107; + if (ixl != 0) + { + /* If the high double is an exact power of two and the low + double has the opposite sign, then the exponent calculated + from the high double is one too big. */ + if (ix == 0 + && (int64_t) (hx ^ lx) < 0) + { + hx += 1L << 52; + expon -= 1; } - *eptr += (ix>>52)-1022; - if (ixl != 0ULL) { - explo = (ixl>>52) - (ix>>52) + 0x3fe; - if ((ixl&0x7ff0000000000000ULL) == 0LL) { - /* the lower double is a denormal so we need to correct its - mantissa and perhaps its exponent. */ - int cnt; + explo = ixl >> 52; + if (explo == 0) + { + /* The low double started out as a denormal. Normalize its + mantissa and adjust the exponent. */ + int cnt; - if (sizeof (ixl) == sizeof (long)) - cnt = __builtin_clzl (ixl); - else if ((ixl >> 32) != 0) - cnt = __builtin_clzl ((long) (ixl >> 32)); - else - cnt = __builtin_clzl ((long) ixl) + 32; - cnt = cnt - 12; - lx = (lx&0x8000000000000000ULL) | ((explo-cnt)<<52) - | ((ixl<<(cnt+1))&0x000fffffffffffffULL); - } else - lx = (lx&0x800fffffffffffffULL) | (explo<<52); - } else - lx = 0ULL; + if (sizeof (ixl) == sizeof (long)) + cnt = __builtin_clzl (ixl); + else if ((ixl >> 32) != 0) + cnt = __builtin_clzl ((long) (ixl >> 32)); + else + cnt = __builtin_clzl ((long) ixl) + 32; + cnt = cnt - 12; + explo -= cnt; + ixl <<= cnt + 1; + } + + /* With variable precision we can't assume much about the + magnitude of the returned low double. It may even be a + denormal. */ + explo -= expon; + ixl &= 0x000fffffffffffffULL; + lx &= 0x8000000000000000ULL; + if (explo <= 0) + { + /* Handle denormal low double. */ + if (explo > -52) + { + ixl |= 1LL << 52; + ixl >>= 1 - explo; + } + else + { + ixl = 0; + lx = 0; + if ((hx & 0x7ff0000000000000ULL) == (1023LL << 52)) + { + /* Oops, the adjustment we made above for values a + little smaller than powers of two turned out to + be wrong since the returned low double will be + zero. This can happen if the input was + something weird like 0x1p1000 - 0x1p-1000. */ + hx -= 1L << 52; + expon += 1; + } + } + explo = 0; + } + lx |= (explo << 52) | ixl; + } - hx = (hx&0x800fffffffffffffULL) | 0x3fe0000000000000ULL; - INSERT_WORDS64 (xhi, hx); - INSERT_WORDS64 (xlo, lx); - x = ldbl_pack (xhi, xlo); - return x; + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); + *eptr = expon; + return x; } #ifdef IS_IN_libm long_double_symbol (libm, __frexpl, frexpl); |