Fix strtod multiple-precision division bug (bug 26137).

Bug 26137 reports spurious "inexact" exceptions from strtod, on 32-bit
systems only, for a decimal argument that is exactly 1 + 2^-32.  In
fact the same issue also appears for 1 + 2^-64 and 1 + 2^-96 as
arguments to strtof128 on 32-bit systems, and 1 + 2^-64 as an argument
to strtof128 on 64-bit systems.  In FE_DOWNWARD or FE_TOWARDZERO mode,
the return value is also incorrect.

The problem is in the multiple-precision division logic used in the
case of dividing by a denominator that occupies at least three GMP
limbs.  There was a comment "The division does not work if the upper
limb of the two-limb mumerator is greater than the denominator.", but
in fact there were problems for the case of equality (that is, where
the high limbs are equal, offset by some multiple of the GMP limb
size) as well.  In such cases, the code used "quot = ~(mp_limb_t) 0;"
(with subsequent correction if that is an overestimate), because
udiv_qrnnd does not support the case of equality, but it's possible
for the shifted numerator to be greater than or equal to the
denominator, in which case that is an underestimate.  To avoid that,
this patch changes the ">" condition to ">=", meaning the first
division is done with a zero high word.

The tests added are all 1 + 2^-n for n from 1 to 113 except for those
that were already present in tst-strtod-round-data.

Tested for x86_64 and x86.

1 files changed, 2 insertions, 2 deletions

diff --git a/stdlib/strtod_l.c b/stdlib/strtod_l.c
index a3836fc12e..158da787a2 100644
--- a/stdlib/strtod_l.c
+++ b/stdlib/strtod_l.c
@@ -1648,8 +1648,8 @@ ____STRTOF_INTERNAL (const STRING_TYPE *nptr, STRING_TYPE **endptr, int group,
 	  d1 = den[densize - 2];
 
 	  /* The division does not work if the upper limb of the two-limb
-	     numerator is greater than the denominator.  */
-	  if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
+	     numerator is greater or equal to than the denominator.  */
+	  if (__mpn_cmp (num, &den[densize - numsize], numsize) >= 0)
 	    num[numsize++] = 0;
 
 	  if (numsize < densize)