Re: libquadmath fix for 94756 and 87204

[Jakub Jelinek <jakub@redhat.com>]

On Sat, Jan 21, 2023 at 04:31:52PM +0000, i.nixman--- via Gcc-patches wrote:
> > Why?
> 
> could you explain which of the nine lines are you talking about?

All the uselessly changed ones.

> > As for the rest, it would help if you could list the exact glibc commits
> > which you've ported to libquadmath and indicate if it is solely those
> > and nothing else.
> 
> I'm sorry but it was not my intention to remember exactly which commits I
> was looking at...
> I didn't think about it.

Because this isn't a complete update to latest glibc source, I think it is
quite important to document what has been merged for reproducibility etc.

> --- a/libquadmath/printf/gmp-impl.h
> +++ b/libquadmath/printf/gmp-impl.h
> @@ -33,15 +33,24 @@ MA 02111-1307, USA. */
>  #define MAX(h,i) ((h) > (i) ? (h) : (i))
>  #endif
>  
> -#define BITS_PER_MP_LIMB (__SIZEOF_LONG__ * __CHAR_BIT__)
> -#define BYTES_PER_MP_LIMB (BITS_PER_MP_LIMB / __CHAR_BIT__)
> -typedef unsigned long int	mp_limb_t;
> -typedef long int		mp_limb_signed_t;
> +#ifdef __MINGW32__ && defined(__x86_64__)

This isn't valid in C, so it is unclear how it was actually tested.
Furthermore, as I said earlier, I think we should change only the
absolute minimum, not reformat stuff, etc.
And, I really don't like to base this on mingw32, doesn't say cygwin
want to do the same?  So, I think it is much better to base this
on LLP64.

> --- a/libquadmath/strtod/strtod_l.c
> +++ b/libquadmath/strtod/strtod_l.c

As for the strtod_l.c changes, I went through all of them, found
the corresponding glibc commits and backported them one by one
(plus one follow-up).  Changes in this file from your version to
mine are in the second hunk which was weirdly indented compared to
surrounding code, and some comments.

So, here is what I'd like to commit instead.
I've tested on x86_64-linux -m32/-m64 the testcases from both of the
PRs without/with the patch and while without the patch the output
differed between -m32/-m64, with the patch it doesn't.

I'll bootstrap/regtest this tonight on x86_64-linux and i686-linux,
though I don't have a way to test this on mingw32 nor cygwin.

2023-03-01  niXman  <i.nixman@autistici.org>
	    Jakub Jelinek  <jakub@redhat.com>

	PR libquadmath/87204
	PR libquadmath/94756
	* printf/gmp-impl.h (mp_limb_t, mp_limb_signed_t, BITS_PER_MP_LIMB):
	Use 64-bit limbs on LLP64 targets.
	* strtod/strtod_l.c (round_and_return): Cherry-pick glibc
	9310c284ae9 BZ #16151, 4406c41c1d6 BZ #16965 and fcd6b5ac36a
	BZ #23279 fixes.
	(____STRTOF_INTERNAL): Cherry-pick glibc b0debe14fcf BZ #23007,
	5556d30caee BZ #18247, 09555b9721d and c6aac3bf366 BZ #26137 and
	d84f25c7d87 fixes.

--- libquadmath/printf/gmp-impl.h.jj	2020-01-12 11:54:39.787362505 +0100
+++ libquadmath/printf/gmp-impl.h	2023-03-01 18:40:07.696374919 +0100
@@ -33,10 +33,18 @@ MA 02111-1307, USA. */
 #define MAX(h,i) ((h) > (i) ? (h) : (i))
 #endif
 
+#if __SIZEOF_LONG__ == 4 && __SIZEOF_LONG_LONG__ == 8 \
+    && __SIZEOF_POINTER__ == 8
+/* Use 64-bit limbs on LLP64 targets.  */
+#define BITS_PER_MP_LIMB (__SIZEOF_LONG_LONG__ * __CHAR_BIT__)
+typedef unsigned long long int	mp_limb_t;
+typedef long long int		mp_limb_signed_t;
+#else
 #define BITS_PER_MP_LIMB (__SIZEOF_LONG__ * __CHAR_BIT__)
-#define BYTES_PER_MP_LIMB (BITS_PER_MP_LIMB / __CHAR_BIT__)
 typedef unsigned long int	mp_limb_t;
 typedef long int		mp_limb_signed_t;
+#endif
+#define BYTES_PER_MP_LIMB (BITS_PER_MP_LIMB / __CHAR_BIT__)
 
 typedef mp_limb_t *             mp_ptr;
 typedef const mp_limb_t *	mp_srcptr;
--- libquadmath/strtod/strtod_l.c.jj	2020-01-12 11:54:39.788362490 +0100
+++ libquadmath/strtod/strtod_l.c	2023-03-01 18:59:49.490151048 +0100
@@ -200,7 +200,7 @@ round_and_return (mp_limb_t *retval, int
 
 	  round_limb = retval[RETURN_LIMB_SIZE - 1];
 	  round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
-	  for (i = 0; i < RETURN_LIMB_SIZE; ++i)
+	  for (i = 0; i < RETURN_LIMB_SIZE - 1; ++i)
 	    more_bits |= retval[i] != 0;
 	  MPN_ZERO (retval, RETURN_LIMB_SIZE);
 	}
@@ -215,9 +215,14 @@ round_and_return (mp_limb_t *retval, int
 	  more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
 			!= 0);
 
-	  (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
-			     RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
-			     shift % BITS_PER_MP_LIMB);
+	  /* mpn_rshift requires 0 < shift < BITS_PER_MP_LIMB.  */
+	  if ((shift % BITS_PER_MP_LIMB) != 0)
+	    (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
+			       RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
+			       shift % BITS_PER_MP_LIMB);
+	  else
+	    for (i = 0; i < RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB); i++)
+	      retval[i] = retval[i + (shift / BITS_PER_MP_LIMB)];
 	  MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
 		    shift / BITS_PER_MP_LIMB);
 	}
@@ -276,7 +281,7 @@ round_and_return (mp_limb_t *retval, int
 	}
     }
 
-  if (exponent > MAX_EXP)
+  if (exponent >= MAX_EXP)
     goto overflow;
 
 #ifdef HAVE_FENV_H
@@ -308,7 +313,7 @@ round_and_return (mp_limb_t *retval, int
     }
 #endif
 
-  if (exponent > MAX_EXP)
+  if (exponent >= MAX_EXP)
   overflow:
     return overflow_value (negative);
 
@@ -688,7 +693,7 @@ ____STRTOF_INTERNAL (nptr, endptr, group
 	  if (endptr != NULL)
 	    *endptr = (STRING_TYPE *) cp;
 
-	  return retval;
+	  return negative ? -retval : retval;
 	}
 
       /* It is really a text we do not recognize.  */
@@ -1193,7 +1198,16 @@ ____STRTOF_INTERNAL (nptr, endptr, group
   if (__builtin_expect (exponent > MAX_10_EXP + 1 - (intmax_t) int_no, 0))
     return overflow_value (negative);
 
-  if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0))
+  /* 10^(MIN_10_EXP-1) is not normal.  Thus, 10^(MIN_10_EXP-1) /
+     2^MANT_DIG is below half the least subnormal, so anything with a
+     base-10 exponent less than the base-10 exponent (which is
+     MIN_10_EXP - 1 - ceil(MANT_DIG*log10(2))) of that value
+     underflows.  DIG is floor((MANT_DIG-1)log10(2)), so an exponent
+     below MIN_10_EXP - (DIG + 3) underflows.  But EXPONENT is
+     actually an exponent multiplied only by a fractional part, not an
+     integer part, so an exponent below MIN_10_EXP - (DIG + 2)
+     underflows.  */
+  if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 2), 0))
     return underflow_value (negative);
 
   if (int_no > 0)
@@ -1360,7 +1374,7 @@ ____STRTOF_INTERNAL (nptr, endptr, group
 
     assert (dig_no > int_no
 	    && exponent <= 0
-	    && exponent >= MIN_10_EXP - (DIG + 1));
+	    && exponent >= MIN_10_EXP - (DIG + 2));
 
     /* We need to compute MANT_DIG - BITS fractional bits that lie
        within the mantissa of the result, the following bit for
@@ -1651,8 +1665,8 @@ ____STRTOF_INTERNAL (nptr, endptr, 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 than or equal to the denominator.  */
+	  if (mpn_cmp (num, &den[densize - numsize], numsize) >= 0)
 	    num[numsize++] = 0;
 
 	  if (numsize < densize)
@@ -1761,7 +1775,7 @@ ____STRTOF_INTERNAL (nptr, endptr, group
 	      got_limb;
 	    }
 
-	  for (i = densize; num[i] == 0 && i >= 0; --i)
+	  for (i = densize; i >= 0 && num[i] == 0; --i)
 	    ;
 	  return round_and_return (retval, exponent - 1, negative,
 				   quot, BITS_PER_MP_LIMB - 1 - used,


	Jakub