From 4cf69995e26e16005d4e3843ad4d18c75cf21a04 Mon Sep 17 00:00:00 2001 From: Alan Modra Date: Sat, 17 Aug 2013 18:19:44 +0930 Subject: 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. --- sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c | 44 ++++++++++++++++------------------ 1 file changed, 21 insertions(+), 23 deletions(-) (limited to 'sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c') diff --git a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c index c3e42f2..c96852d 100644 --- a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c @@ -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; } -- cgit v1.1