AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-22 21:10:07 +00:00
Code Issues Packages Projects Releases Wiki Activity

Correct IBM long double frexpl.

Browse Source 
Besides fixing the bugzilla, this also fixes corner-cases where the high
and low double differ greatly in magnitude, and handles a denormal
input without resorting to a fp rescale.

	[BZ #16740]
	[BZ #16619]
	* sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite.
	* math/libm-test.inc (frexp_test_data): Add tests.
This commit is contained in:
Alan Modra 2014-04-16 19:33:32 +09:30
parent bb9c256fb0
commit aa5f0ff11a
3 changed files with 118 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
ChangeLog
View File

@ -1,3 +1,10 @@
2014-04-16 Alan Modra <amodra@gmail.com>
[BZ #16740]
[BZ #16619]
* sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite.
* math/libm-test.inc (frexp_test_data): Add tests.
2014-04-15 Siddhesh Poyarekar <siddhesh@redhat.com>
* benchtests/Makefile: Depend on libraries in build directory.

9
math/libm-test.inc
View File

@ -7186,6 +7186,15 @@ static const struct test_f_f1_data frexp_test_data[] =
TEST_fI_f1 (frexp, 12.8L, 0.8L, 4, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, -27.34L, -0.854375L, 5, NO_INEXACT_EXCEPTION),
#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 106
TEST_fI_f1 (frexp, 1.0L-0x1p-106L, 1.0L-0x1p-106L, 0, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, 1.0L, 0.5L, 1, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, 1.0L+0x1p-105L, 0.5L+0x1p-106L, 1, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, -1.0L+0x1p-106L, -1.0L+0x1p-106L, 0, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, -1.0L, -0.5L, 1, NO_INEXACT_EXCEPTION),
TEST_fI_f1 (frexp, -1.0L-0x1p-105L, -0.5L-0x1p-106L, 1, NO_INEXACT_EXCEPTION),
#endif
};
static void

112
sysdeps/ieee754/ldbl-128ibm/s_frexpl.c
View File

@ -31,36 +31,63 @@ 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;
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;
*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;
ixl = 0x7fffffffffffffffULL & lx;
ix = 0x7fffffffffffffffULL & hx;
expon = 0;
if (ix >= 0x7ff0000000000000ULL || ix == 0)
{
/* 0,inf,nan. */
*eptr = expon;
return x;
}
*eptr += (ix>>52)-1022;
expon = ix >> 52;
if (expon == 0)
{
/* Denormal high double, the low double must be 0.0. */
int cnt;
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. */
/* 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;
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;
}
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))
@ -70,17 +97,48 @@ long double __frexpl(long double x, int *eptr)
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;
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);
*eptr = expon;
return x;
}
#ifdef IS_IN_libm