Message ID | 56EB246E.4050504@linux.vnet.ibm.com |
---|---|
State | New |
Headers | show |
"Paul E. Murphy" <murphyp@linux.vnet.ibm.com> writes: > + /* Hold 7 extra bits of precision in the mantissa. This allows > + the normalizing shifts below to prevent losing precision when > + the signs differ and the exponents are sufficiently far apart. */ > + lo <<= 7; That magical number should probably get a name. > /* we have a borrow from the hidden bit, so shift left 1. */ > - hi = 0x0ffffffffffffeLL | (lo >> 51); > - lo = 0x1fffffffffffffLL & (lo << 1); > + hi = 0x000ffffffffffffeLL | (lo >> (52 + 7)); > + lo = 0x0fffffffffffffffLL & (lo << 1); Not 51 + 7? Andreas.
On Thu, 17 Mar 2016, Paul E. Murphy wrote: > When the signs differ, the precision of the conversion > could drop below 106 bits. This strategy is > identical to the hexadecimal variant. I assume this was user-visible in a release, in which case it should have a bug filed in Bugzilla (if it doesn't already have one - I take it this is distinct from bug 5268?). It should also have a testcase included in the patch (there are several tests in sysdeps/ieee754/ldbl-128ibm already for issues specific to this format - various of which use unions to test issues that can't be tested simply with 106-bit inputs GCC can represent, or where it's necessary to examine the two halves of a long double result separately).
On 03/17/2016 04:56 PM, Andreas Schwab wrote: > "Paul E. Murphy" <murphyp@linux.vnet.ibm.com> writes: > >> + /* Hold 7 extra bits of precision in the mantissa. This allows >> + the normalizing shifts below to prevent losing precision when >> + the signs differ and the exponents are sufficiently far apart. */ >> + lo <<= 7; > > That magical number should probably get a name. Suggestions? It was noted as an "IEEE 854 adjustment" constant in printf_fphex.c, and otherwise integrated into the shifts later on. INTERNAL_PRECISION_SHIFT? Setting it 1 should fix the I case noted in the comments. > >> /* we have a borrow from the hidden bit, so shift left 1. */ >> - hi = 0x0ffffffffffffeLL | (lo >> 51); >> - lo = 0x1fffffffffffffLL & (lo << 1); >> + hi = 0x000ffffffffffffeLL | (lo >> (52 + 7)); >> + lo = 0x0fffffffffffffffLL & (lo << 1); > > Not 51 + 7? That is intentional. It should be moving the high bit of lo lo into position 0. This is identical to how the ldbl-128ibm printf_fphex.c, which does round trip these numbers correctly. I.e ((1<<53)-1)>>51 == 3. The implicit bit is already factored into lo, so the existing shift appears wrong here. > > Andreas. >
On 03/17/2016 05:16 PM, Joseph Myers wrote: > On Thu, 17 Mar 2016, Paul E. Murphy wrote: > >> When the signs differ, the precision of the conversion >> could drop below 106 bits. This strategy is >> identical to the hexadecimal variant. > > I assume this was user-visible in a release, in which case it should have > a bug filed in Bugzilla (if it doesn't already have one - I take it this > is distinct from bug 5268?). > > It should also have a testcase included in the patch (there are several > tests in sysdeps/ieee754/ldbl-128ibm already for issues specific to this > format - various of which use unions to test issues that can't be tested > simply with 106-bit inputs GCC can represent, or where it's necessary to > examine the two halves of a long double result separately). > This came up internally. I think it is similar, the comments would indicate it is a 1 bit internal precision loss which this patch does address. I will figure out how to update the tests to catch this, and give the magic shift amount a name. Thanks, Paul
diff --git a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c index 4f550ef..afe6cf2 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c @@ -45,12 +45,17 @@ __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; + /* Hold 7 extra bits of precision in the mantissa. This allows + the normalizing shifts below to prevent losing precision when + the signs differ and the exponents are sufficiently far apart. */ + lo <<= 7; + /* 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; + lo |= 1ULL << (52 + 7); else - lo = lo << 1; + lo = lo << (1 + 7); /* The lower double is normalized separately from the upper. We may need to adjust the lower manitissa to reflect this. */ @@ -72,12 +77,12 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, if (u.d[0].ieee.negative != u.d[1].ieee.negative && lo != 0) { - lo = (1ULL << 53) - lo; + lo = (1ULL << (53 + 7)) - lo; if (hi == 0) { /* we have a borrow from the hidden bit, so shift left 1. */ - hi = 0x0ffffffffffffeLL | (lo >> 51); - lo = 0x1fffffffffffffLL & (lo << 1); + hi = 0x000ffffffffffffeLL | (lo >> (52 + 7)); + lo = 0x0fffffffffffffffLL & (lo << 1); (*expt)--; } else @@ -85,14 +90,14 @@ __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, } #if BITS_PER_MP_LIMB == 32 /* Combine the mantissas to be contiguous. */ - res_ptr[0] = lo; - res_ptr[1] = (hi << (53 - 32)) | (lo >> 32); + res_ptr[0] = lo >> 7; + res_ptr[1] = (hi << (53 - 32)) | (lo >> (32 + 7)); res_ptr[2] = hi >> 11; res_ptr[3] = hi >> (32 + 11); #define N 4 #elif BITS_PER_MP_LIMB == 64 /* Combine the two mantissas to be contiguous. */ - res_ptr[0] = (hi << 53) | lo; + res_ptr[0] = (hi << 53) | (lo >> 7); res_ptr[1] = hi >> 11; #define N 2 #else