From 37a4c70bd4c5c74ac562072e450dc02e8cb4c150 Mon Sep 17 00:00:00 2001 From: "Paul E. Murphy" Date: Mon, 29 Feb 2016 13:27:36 -0600 Subject: [PATCH] Increase internal precision of ldbl-128ibm decimal printf [BZ #19853] When the signs differ, the precision of the conversion sometimes drops below 106 bits. This strategy is identical to the hexadecimal variant. I've refactored tst-sprintf3 to enable testing a value with more than 30 significant digits in order to demonstrate this failure and its solution. Additionally, this implicitly fixes a typo in the shift quantities when subtracting from the high mantissa to compute the difference. --- ChangeLog | 14 ++++++++++++++ stdio-common/tst-sprintf3.c | 14 ++++++++++---- sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c | 25 ++++++++++++++++++------- 3 files changed, 42 insertions(+), 11 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1a87d43e11..29dd164802 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,17 @@ +2016-03-31 Paul E. Murphy + + [BZ #19853] + * stdio-common/tst-sprintf3.c [TEST_N]: Refactor + TEST to take significant digits as second parameter. + [TEST]: Redefine in terms of TEST_N taking 30 + significant digits. + (do_test): Add test case to demonstrate precision + failure in the ldbl-128ibm printf. + * sysdeps/ieee754/ldbl-128ibm/ldbl2pm.c: + (__mpn_extract_long_double): Carry 7 extra intermediate + bits of precision to aide computing difference when + signs differ. + 2016-03-31 H.J. Lu [BZ #19881] diff --git a/stdio-common/tst-sprintf3.c b/stdio-common/tst-sprintf3.c index cffd1b69d0..dad216b941 100644 --- a/stdio-common/tst-sprintf3.c +++ b/stdio-common/tst-sprintf3.c @@ -38,11 +38,11 @@ do_test (void) # define COMPARE_LDBL(u, v) ((u).l == (v).l) #endif -#define TEST(val) \ +#define TEST_N(val, n) \ do \ { \ u.l = (val); \ - snprintf (buf, sizeof buf, "%.30LgL", u.l); \ + snprintf (buf, sizeof buf, "%." #n "LgL", u.l); \ if (strcmp (buf, #val) != 0) \ { \ printf ("Error on line %d: %s != %s\n", __LINE__, buf, #val); \ @@ -50,19 +50,25 @@ do_test (void) } \ if (sscanf (#val, "%Lg", &v.l) != 1 || !COMPARE_LDBL (u, v)) \ { \ - printf ("Error sscanf on line %d: %.30Lg != %.30Lg\n", __LINE__, \ - u.l, v.l); \ + printf ("Error sscanf on line %d: %." #n "Lg != %." #n "Lg\n", \ + __LINE__, u.l, v.l); \ result = 1; \ } \ /* printf ("%s %Lg %016Lx %016Lx\n", #val, u.l, u.x[0], u.x[1]); */ \ } \ while (0) +#define TEST(val) TEST_N (val,30) + #if LDBL_MANT_DIG >= 106 # if LDBL_MANT_DIG == 106 TEST (2.22507385850719347803989925739e-308L); TEST (2.22507385850719397210554509863e-308L); TEST (2.22507385850720088902458687609e-308L); + + /* Verify precision is not lost for long doubles + of the form +1.pN,-1.pM. */ + TEST_N (3.32306998946228968225951765070082e+35L, 34); # endif TEST (2.22507385850720138309023271733e-308L); TEST (2.22507385850720187715587855858e-308L); diff --git a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c index 4f550ef47c..30d9bc3443 100644 --- a/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c @@ -28,6 +28,12 @@ bits (106 for long double) and an integral power of two (MPN frexpl). */ + +/* When signs differ, the actual value is the difference between the + significant double and the less significant double. Sometimes a + bit can be lost when we borrow from the significant mantissa. */ +#define EXTRA_INTERNAL_PRECISION (7) + mp_size_t __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, int *expt, int *is_neg, @@ -45,10 +51,15 @@ __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 <<= EXTRA_INTERNAL_PRECISION; + /* 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 + EXTRA_INTERNAL_PRECISION); else lo = lo << 1; @@ -72,12 +83,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 + EXTRA_INTERNAL_PRECISION)) - 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 + EXTRA_INTERNAL_PRECISION)); + lo = 0x0fffffffffffffffLL & (lo << 1); (*expt)--; } else @@ -85,14 +96,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 >> EXTRA_INTERNAL_PRECISION; + res_ptr[1] = (hi << (53 - 32)) | (lo >> (32 + EXTRA_INTERNAL_PRECISION)); 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 >> EXTRA_INTERNAL_PRECISION); res_ptr[1] = hi >> 11; #define N 2 #else