Fix for [BZ #15680] IBM long double inaccuracy

http://sourceware.org/ml/libc-alpha/2013-06/msg00919.html

I discovered a number of places where denormals and other corner cases
were being handled wrongly.

- printf_fphex.c: Testing for the low double exponent being zero is
unnecessary.  If the difference in exponents is less than 53 then the
high double exponent must be nearing the low end of its range, and the
low double exponent hit rock bottom.

- ldbl2mpn.c: A denormal (ie. exponent of zero) value is treated as
if the exponent was one, so shift mantissa left by one.  Code handling
normalisation of the low double mantissa lacked a test for shift count
greater than bits in type being shifted, and lacked anything to handle
the case where the difference in exponents is less than 53 as in
printf_fphex.c.

- math_ldbl.h (ldbl_extract_mantissa): Same as above, but worse, with
code testing for exponent > 1 for some reason, probably a typo for >= 1.

- math_ldbl.h (ldbl_insert_mantissa): Round the high double as per
mpn2ldbl.c (hi is odd or explicit mantissas non-zero) so that the
number we return won't change when applying ldbl_canonicalize().
Add missing overflow checks and normalisation of high mantissa.
Correct misleading comment: "The hidden bit of the lo mantissa is
zero" is not always true as can be seen from the code rounding the hi
mantissa.  Also by inspection, lzcount can never be less than zero so
remove that test.  Lastly, masking bitfields to their widths can be
left to the compiler.

- mpn2ldbl.c: The overflow checks here on rounding of high double were
just plain wrong.  Incrementing the exponent must be accompanied by a
shift right of the mantissa to keep the value unchanged.  Above notes
for ldbl_insert_mantissa are also relevant.

	[BZ #15680]
	* sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c: Comment fix.
	* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c
	(PRINT_FPHEX_LONG_DOUBLE): Tidy code by moving -53 into ediff
	calculation.  Remove unnecessary test for denormal exponent.
	* sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double):
	Correct handling of denormals.  Avoid undefined shift behaviour.
	Correct normalisation of low mantissa when low double is denormal.
	* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
	(ldbl_extract_mantissa): Likewise.  Comment.  Use uint64_t* for hi64.
	(ldbl_insert_mantissa): Make both hi64 and lo64 parms uint64_t.
	Correct normalisation of low mantissa.  Test for overflow of high
	mantissa and normalise.
	(ldbl_nearbyint): Use more readable constant for two52.
	* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
	(__mpn_construct_long_double): Fix test for overflow of high
	mantissa and correct normalisation.  Avoid undefined shift.
This commit is contained in:
Alan Modra 2013-08-17 18:19:44 +09:30
parent 9605ca6c08
commit 4cf69995e2
6 changed files with 177 additions and 107 deletions

View File

@ -1,3 +1,23 @@
2013-10-04 Alan Modra <amodra@gmail.com>
[BZ #15680]
* sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c: Comment fix.
* sysdeps/ieee754/ldbl-128ibm/printf_fphex.c
(PRINT_FPHEX_LONG_DOUBLE): Tidy code by moving -53 into ediff
calculation. Remove unnecessary test for denormal exponent.
* sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double):
Correct handling of denormals. Avoid undefined shift behaviour.
Correct normalisation of low mantissa when low double is denormal.
* sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
(ldbl_extract_mantissa): Likewise. Comment. Use uint64_t* for hi64.
(ldbl_insert_mantissa): Make both hi64 and lo64 parms uint64_t.
Correct normalisation of low mantissa. Test for overflow of high
mantissa and normalise.
(ldbl_nearbyint): Use more readable constant for two52.
* sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
(__mpn_construct_long_double): Fix test for overflow of high
mantissa and correct normalisation. Avoid undefined shift.
2013-10-04 Alan Modra <amodra@gmail.com>
* sysdeps/ieee754/ldbl-128ibm/ieee754.h

View File

@ -243,7 +243,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
We split the 113 bits of the mantissa into 5 24bit integers
stored in a double array. */
/* Make the IBM extended format 105 bit mantissa look like the ieee854 112
bit mantissa so the next operatation will give the correct result. */
bit mantissa so the next operation will give the correct result. */
ldbl_extract_mantissa (&ixd, &lxd, &exp, x);
exp = exp - 23;
/* This is faster than doing this in floating point, because we

View File

@ -36,6 +36,7 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
union ibm_extended_long_double u;
unsigned long long hi, lo;
int ediff;
u.ld = value;
*is_neg = u.d[0].ieee.negative;
@ -43,27 +44,36 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
/* If the lower double is not a denomal or zero then set the hidden
53rd bit. */
if (u.d[1].ieee.exponent > 0)
{
lo |= 1LL << 52;
/* The lower double is normalized separately from the upper. We may
need to adjust the lower manitissa to reflect this. */
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;
if (ediff > 53)
lo = lo >> (ediff-53);
/* If the lower double is not a denormal or zero then set the hidden
53rd bit. */
if (u.d[1].ieee.exponent != 0)
lo |= 1ULL << 52;
else
lo = lo << 1;
/* The lower double is normalized separately from the upper. We may
need to adjust the lower manitissa to reflect this. */
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
if (ediff > 0)
{
if (ediff < 64)
lo = lo >> ediff;
else
lo = 0;
}
else if (ediff < 0)
lo = lo << -ediff;
/* The high double may be rounded and the low double reflects the
difference between the long double and the rounded high double
value. This is indicated by a differnce between the signs of the
high and low doubles. */
if ((u.d[0].ieee.negative != u.d[1].ieee.negative)
&& ((u.d[1].ieee.exponent != 0) && (lo != 0L)))
if (u.d[0].ieee.negative != u.d[1].ieee.negative
&& lo != 0)
{
lo = (1ULL << 53) - lo;
if (hi == 0LL)
if (hi == 0)
{
/* we have a borrow from the hidden bit, so shift left 1. */
hi = 0x0ffffffffffffeLL | (lo >> 51);

View File

@ -6,6 +6,10 @@
#include <ieee754.h>
#include <stdint.h>
/* To suit our callers we return *hi64 and *lo64 as if they came from
an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one
implicit bit) and 64 bits in *lo64. */
static inline void
ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
{
@ -14,77 +18,119 @@ ldbl_extract_mantissa (int64_t *hi64, uint64_t *lo64, int *exp, long double x)
the number before the decimal point and the second implicit bit
as bit 53 of the mantissa. */
uint64_t hi, lo;
int ediff;
union ibm_extended_long_double u;
u.ld = x;
*exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;
lo = ((uint64_t) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;
hi = ((uint64_t) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;
/* If the lower double is not a denomal or zero then set the hidden
53rd bit. */
if (u.d[1].ieee.exponent > 0x001)
{
lo |= (1ULL << 52);
lo = lo << 7; /* pre-shift lo to match ieee854. */
/* The lower double is normalized separately from the upper. We
may need to adjust the lower manitissa to reflect this. */
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;
if (ediff > 53)
lo = lo >> (ediff-53);
hi |= (1ULL << 52);
}
if ((u.d[0].ieee.negative != u.d[1].ieee.negative)
&& ((u.d[1].ieee.exponent != 0) && (lo != 0LL)))
if (u.d[0].ieee.exponent != 0)
{
hi--;
lo = (1ULL << 60) - lo;
if (hi < (1ULL << 52))
int ediff;
/* If not a denormal or zero then we have an implicit 53rd bit. */
hi |= (uint64_t) 1 << 52;
if (u.d[1].ieee.exponent != 0)
lo |= (uint64_t) 1 << 52;
else
/* A denormal is to be interpreted as having a biased exponent
of 1. */
lo = lo << 1;
/* We are going to shift 4 bits out of hi later, because we only
want 48 bits in *hi64. That means we want 60 bits in lo, but
we currently only have 53. Shift the value up. */
lo = lo << 7;
/* The lower double is normalized separately from the upper.
We may need to adjust the lower mantissa to reflect this.
The difference between the exponents can be larger than 53
when the low double is much less than 1ULP of the upper
(in which case there are significant bits, all 0's or all
1's, between the two significands). The difference between
the exponents can be less than 53 when the upper double
exponent is nearing its minimum value (in which case the low
double is denormal ie. has an exponent of zero). */
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;
if (ediff > 0)
{
/* we have a borrow from the hidden bit, so shift left 1. */
hi = (hi << 1) | (lo >> 59);
lo = 0xfffffffffffffffLL & (lo << 1);
*exp = *exp - 1;
if (ediff < 64)
lo = lo >> ediff;
else
lo = 0;
}
else if (ediff < 0)
lo = lo << -ediff;
if (u.d[0].ieee.negative != u.d[1].ieee.negative
&& lo != 0)
{
hi--;
lo = ((uint64_t) 1 << 60) - lo;
if (hi < (uint64_t) 1 << 52)
{
/* We have a borrow from the hidden bit, so shift left 1. */
hi = (hi << 1) | (lo >> 59);
lo = (((uint64_t) 1 << 60) - 1) & (lo << 1);
*exp = *exp - 1;
}
}
}
else
/* If the larger magnitude double is denormal then the smaller
one must be zero. */
hi = hi << 1;
*lo64 = (hi << 60) | lo;
*hi64 = hi >> 4;
}
static inline long double
ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64)
{
union ibm_extended_long_double u;
unsigned long hidden2, lzcount;
unsigned long long hi, lo;
int expnt2;
uint64_t hi, lo;
u.d[0].ieee.negative = sign;
u.d[1].ieee.negative = sign;
u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;
u.d[1].ieee.exponent = exp-53 + IEEE754_DOUBLE_BIAS;
u.d[1].ieee.exponent = 0;
expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS;
/* Expect 113 bits (112 bits + hidden) right justified in two longs.
The low order 53 bits (52 + hidden) go into the lower double */
lo = (lo64 >> 7)& ((1ULL << 53) - 1);
hidden2 = (lo64 >> 59) & 1ULL;
lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1);
/* The high order 53 bits (52 + hidden) go into the upper double */
hi = (lo64 >> 60) & ((1ULL << 11) - 1);
hi |= (hi64 << 4);
hi = lo64 >> 60;
hi |= hi64 << 4;
if (lo != 0LL)
if (lo != 0)
{
/* hidden2 bit of low double controls rounding of the high double.
If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)
int lzcount;
/* hidden bit of low double controls rounding of the high double.
If hidden is '1' and either the explicit mantissa is non-zero
or hi is odd, then round up hi and adjust lo (2nd mantissa)
plus change the sign of the low double to compensate. */
if (hidden2)
if ((lo & ((uint64_t) 1 << 52)) != 0
&& ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0))
{
hi++;
if ((hi & ((uint64_t) 1 << 53)) != 0)
{
hi = hi >> 1;
u.d[0].ieee.exponent++;
}
u.d[1].ieee.negative = !sign;
lo = (1ULL << 53) - lo;
lo = ((uint64_t) 1 << 53) - lo;
}
/* The hidden bit of the lo mantissa is zero so we need to
normalize the it for the low double. Shift it left until the
hidden bit is '1' then adjust the 2nd exponent accordingly. */
/* Normalize the low double. Shift the mantissa left until
the hidden bit is '1' and adjust the exponent accordingly. */
if (sizeof (lo) == sizeof (long))
lzcount = __builtin_clzl (lo);
@ -92,34 +138,30 @@ ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)
lzcount = __builtin_clzl ((long) (lo >> 32));
else
lzcount = __builtin_clzl ((long) lo) + 32;
lzcount = lzcount - 11;
if (lzcount > 0)
lzcount = lzcount - (64 - 53);
lo <<= lzcount;
expnt2 -= lzcount;
if (expnt2 >= 1)
/* Not denormal. */
u.d[1].ieee.exponent = expnt2;
else
{
int expnt2 = u.d[1].ieee.exponent - lzcount;
if (expnt2 >= 1)
{
/* Not denormal. Normalize and set low exponent. */
lo = lo << lzcount;
u.d[1].ieee.exponent = expnt2;
}
/* Is denormal. Note that biased exponent of 0 is treated
as if it was 1, hence the extra shift. */
if (expnt2 > -53)
lo >>= 1 - expnt2;
else
{
/* Is denormal. */
lo = lo << (lzcount + expnt2);
u.d[1].ieee.exponent = 0;
}
lo = 0;
}
}
else
{
u.d[1].ieee.negative = 0;
u.d[1].ieee.exponent = 0;
}
u.d[1].ieee.negative = 0;
u.d[1].ieee.mantissa1 = lo & ((1ULL << 32) - 1);
u.d[1].ieee.mantissa0 = (lo >> 32) & ((1ULL << 20) - 1);
u.d[0].ieee.mantissa1 = hi & ((1ULL << 32) - 1);
u.d[0].ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 1);
u.d[1].ieee.mantissa1 = lo;
u.d[1].ieee.mantissa0 = lo >> 32;
u.d[0].ieee.mantissa1 = hi;
u.d[0].ieee.mantissa0 = hi >> 32;
return u.ld;
}
@ -163,13 +205,13 @@ ldbl_canonicalize (double *a, double *aa)
*aa = xl;
}
/* Simple inline nearbyint (double) function .
/* Simple inline nearbyint (double) function.
Only works in the default rounding mode
but is useful in long double rounding functions. */
static inline double
ldbl_nearbyint (double a)
{
double two52 = 0x10000000000000LL;
double two52 = 0x1p52;
if (__builtin_expect ((__builtin_fabs (a) < two52), 1))
{

View File

@ -69,9 +69,9 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
else
lzcount = __builtin_clzl ((long) val) + 32;
if (hi)
lzcount = lzcount - 11;
lzcount = lzcount - (64 - 53);
else
lzcount = lzcount + 42;
lzcount = lzcount + 53 - (64 - 53);
if (lzcount > u.d[0].ieee.exponent)
{
@ -97,29 +97,27 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
}
}
if (lo != 0L)
if (lo != 0)
{
/* hidden2 bit of low double controls rounding of the high double.
If hidden2 is '1' and either the explicit mantissa is non-zero
/* hidden bit of low double controls rounding of the high double.
If hidden is '1' and either the explicit mantissa is non-zero
or hi is odd, then round up hi and adjust lo (2nd mantissa)
plus change the sign of the low double to compensate. */
if ((lo & (1LL << 52)) != 0
&& ((hi & 1) != 0 || (lo & ((1LL << 52) - 1))))
&& ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
{
hi++;
if ((hi & ((1LL << 52) - 1)) == 0)
if ((hi & (1LL << 53)) != 0)
{
if ((hi & (1LL << 53)) != 0)
hi -= 1LL << 52;
hi >>= 1;
u.d[0].ieee.exponent++;
}
u.d[1].ieee.negative = !sign;
lo = (1LL << 53) - lo;
}
/* The hidden bit of the lo mantissa is zero so we need to normalize
it for the low double. Shift it left until the hidden bit is '1'
then adjust the 2nd exponent accordingly. */
/* Normalize the low double. Shift the mantissa left until
the hidden bit is '1' and adjust the exponent accordingly. */
if (sizeof (lo) == sizeof (long))
lzcount = __builtin_clzl (lo);
@ -127,24 +125,24 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
lzcount = __builtin_clzl ((long) (lo >> 32));
else
lzcount = __builtin_clzl ((long) lo) + 32;
lzcount = lzcount - 11;
if (lzcount > 0)
{
lo = lo << lzcount;
exponent2 = exponent2 - lzcount;
}
lzcount = lzcount - (64 - 53);
lo <<= lzcount;
exponent2 -= lzcount;
if (exponent2 > 0)
u.d[1].ieee.exponent = exponent2;
else
else if (exponent2 > -53)
lo >>= 1 - exponent2;
else
lo = 0;
}
else
u.d[1].ieee.negative = 0;
u.d[1].ieee.mantissa1 = lo & 0xffffffffLL;
u.d[1].ieee.mantissa0 = (lo >> 32) & 0xfffff;
u.d[0].ieee.mantissa1 = hi & 0xffffffffLL;
u.d[0].ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);
u.d[1].ieee.mantissa1 = lo;
u.d[1].ieee.mantissa0 = lo >> 32;
u.d[0].ieee.mantissa1 = hi;
u.d[0].ieee.mantissa0 = hi >> 32;
return u.ld;
}

View File

@ -42,15 +42,15 @@ do { \
lo <<= 1; \
/* The lower double is normalized separately from the upper. We \
may need to adjust the lower manitissa to reflect this. */ \
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent; \
if (ediff > 53 + 63) \
ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53; \
if (ediff > 63) \
lo = 0; \
else if (ediff > 53) \
lo = lo >> (ediff - 53); \
else if (u.d[1].ieee.exponent == 0 && ediff < 53) \
lo = lo << (53 - ediff); \
else if (ediff > 0) \
lo = lo >> ediff; \
else if (ediff < 0) \
lo = lo << -ediff; \
if (u.d[0].ieee.negative != u.d[1].ieee.negative \
&& (u.d[1].ieee.exponent != 0 || lo != 0L)) \
&& lo != 0) \
{ \
lo = (1ULL << 60) - lo; \
if (hi == 0L) \