libgcc: Fix up _BitInt division [PR113604]

The following testcase ends up with SIGFPE in __divmodbitint4.
The problem is a thinko in my attempt to implement Knuth's algorithm.

The algorithm does (where b is 65536, i.e. one larger than what
fits in their unsigned short word):
        // Compute estimate qhat of q[j].
        qhat = (un[j+n]*b + un[j+n-1])/vn[n-1];
        rhat = (un[j+n]*b + un[j+n-1]) - qhat*vn[n-1];
      again:
        if (qhat >= b || qhat*vn[n-2] > b*rhat + un[j+n-2])
        { qhat = qhat - 1;
          rhat = rhat + vn[n-1];
          if (rhat < b) goto again;
        }
The problem is that it uses a double-word / word -> double-word
division (and modulo), while all we have is udiv_qrnnd unless
we'd want to do further library calls, and udiv_qrnnd is a
double-word / word -> word division and modulo.
Now, as the algorithm description says, it can produce at most
word bits + 1 bit quotient.  And I believe that actually the
highest qhat the original algorithm can produce is
(1 << word_bits) + 1.  The algorithm performs earlier canonicalization
where both the divisor and dividend are shifted left such that divisor
has msb set.  If it has msb set already before, no shifting occurs but
we start with added 0 limb, so in the first uv1:uv0 double-word uv1
is 0 and so we can't get too high qhat, if shifting occurs, the first
limb of dividend is shifted right by UWtype bits - shift count into
a new limb, so again in the first iteration in the uv1:uv0 double-word
uv1 doesn't have msb set while vv1 does and qhat has to fit into word.
In the following iterations, previous iteration should guarantee that
the previous quotient digit is correct.  Even if the divisor was the
maximal possible vv1:all_ones_in_all_lower_limbs, if the old uv0:lower_limbs
would be larger or equal to the divisor, the previous quotient digit
would increase and another divisor would be subtracted, which I think
implies that in the next iteration in uv1:uv0 double-word uv1 <= vv1,
but uv0 could be up to all ones, e.g. in case of all lower limbs
of divisor being all ones and at least one dividend limb below uv0
being not all ones.  So, we can e.g. for 64-bit UWtype see
uv1:uv0 / vv1 0x8000000000000000UL:0xffffffffffffffffUL / 0x8000000000000000UL
or 0xffffffffffffffffUL:0xffffffffffffffffUL / 0xffffffffffffffffUL
In all these cases (when uv1 == vv1 && uv0 >= uv1), qhat is
0x10000000000000001UL, i.e. 2 more than fits into UWtype result,
if uv1 == vv1 && uv0 < uv1 it would be 0x10000000000000000UL, i.e.
1 more than fits into UWtype result.
Because we only have udiv_qrnnd which can't deal with those too large
cases (SIGFPEs or otherwise invokes undefined behavior on those), I've
tried to handle the uv1 >= vv1 case separately, but for one thing
I thought it would be at most 1 larger than what fits, and for two
have actually subtracted vv1:vv1 from uv1:uv0 instead of subtracting
0:vv1 from uv1:uv0.
For the uv1 < vv1 case, the implementation already performs roughly
what the algorithm does.
Now, let's see what happens with the two possible extra cases in
the original algorithm.
If uv1 == vv1 && uv0 < uv1, qhat above would be b, so we take
if (qhat >= b, decrement qhat by 1 (it becomes b - 1), add
vn[n-1] aka vv1 to rhat and goto again if rhat < b (but because
qhat already fits we can goto to the again label in the uv1 < vv1
code).  rhat in this case is uv0 and rhat + vv1 can but doesn't
have to overflow, say for uv0 42UL and vv1 0x8000000000000000UL
it will not (and so we should goto again), while for uv0
0x8000000000000000UL and vv1 0x8000000000000001UL it will (and
we shouldn't goto again).
If uv1 == vv1 && uv0 >= uv1, qhat above would be b + 1, so we
take if (qhat >= b, decrement qhat by 1 (it becomes b), add
vn[n-1] aka vv1 to rhat. But because vv1 has msb set and
rhat in this case is uv0 - vv1, the rhat + vv1 addition
certainly doesn't overflow, because (uv0 - vv1) + vv1 is uv0,
so in the algorithm we goto again, again take if (qhat >= b and
decrement qhat so it finally becomes b - 1, and add vn[n-1]
aka vv1 to rhat again.  But this time I believe it must always
overflow, simply because we added (uv0 - vv1) + vv1 + vv1 and
vv1 has msb set, so already vv1 + vv1 must overflow.  And
because it overflowed, it will not goto again.
So, I believe the following patch implements this correctly, by
subtracting vv1 from uv1:uv0 double-word once, then comparing
again if uv1 >= vv1.  If that is true, subtract vv1 from uv1:uv0
again and add 2 * vv1 to rhat, no __builtin_add_overflow is needed
as we know it always overflowed and so won't goto again.
If after the first subtraction uv1 < vv1, use __builtin_add_overflow
when adding vv1 to rhat, because it can but doesn't have to overflow.

I've added an extra testcase which tests the behavior of all the changed
cases, so it has a case where uv1:uv0 / vv1 is 1:1, where it is
1:0 and rhat + vv1 overflows and where it is 1:0 and rhat + vv1 does not
overflow, and includes tests also from Zdenek's other failing tests.

2024-02-02  Jakub Jelinek  <jakub@redhat.com>

	PR libgcc/113604
	* libgcc2.c (__divmodbitint4): If uv1 >= vv1, subtract
	vv1 from uv1:uv0 once or twice as needed, rather than
	subtracting vv1:vv1.

	* gcc.dg/torture/bitint-53.c: New test.
	* gcc.dg/torture/bitint-55.c: New test.

gcc/testsuite/gcc.dg/torture/bitint-53.c

@@ -0,0 +1,26 @@
+/* PR libgcc/113604 */
+/* { dg-do run { target bitint } } */
+/* { dg-options "-std=c23 -pedantic-errors" } */
+/* { dg-skip-if "" { ! run_expensive_tests }  { "*" } { "-O0" "-O2" } } */
+/* { dg-skip-if "" { ! run_expensive_tests } { "-flto" } { "" } } */
+
+#if __BITINT_MAXWIDTH__ >= 256
+unsigned _BitInt (256) x;
+
+void
+foo (unsigned _BitInt (256) a, unsigned _BitInt (128) b)
+{
+  x = a / b;
+}
+#endif
+
+int
+main ()
+{
+#if __BITINT_MAXWIDTH__ >= 256
+  foo (0xfffffffffffffffffffffc0000000000000000000004uwb, 0x7ffffffffffffffffffffffffffuwb);
+  if (x != 0x1fffffffffffffffffuwb)
+    __builtin_abort ();
+#endif
+  return 0;
+}

gcc/testsuite/gcc.dg/torture/bitint-55.c

@@ -0,0 +1,66 @@
+/* PR libgcc/113604 */
+/* { dg-do run { target bitint } } */
+/* { dg-options "-std=c23 -pedantic-errors" } */
+/* { dg-skip-if "" { ! run_expensive_tests }  { "*" } { "-O0" "-O2" } } */
+/* { dg-skip-if "" { ! run_expensive_tests } { "-flto" } { "" } } */
+
+#if __BITINT_MAXWIDTH__ >= 513
+signed _BitInt(513)
+foo (signed _BitInt(513) x, signed _BitInt(513) y)
+{
+  return x % y;
+}
+#endif
+
+#if __BITINT_MAXWIDTH__ >= 512
+unsigned _BitInt(512)
+bar (unsigned _BitInt(512) x, unsigned _BitInt(512) y)
+{
+  return x % y;
+}
+#endif
+
+#if __BITINT_MAXWIDTH__ >= 256
+unsigned _BitInt(256)
+baz (unsigned _BitInt(256) x, unsigned _BitInt(256) y)
+{
+  return x % y;
+}
+#endif
+
+int
+main ()
+{
+#if __BITINT_MAXWIDTH__ >= 513
+  if (foo (11155754932722990178552651944728825929130437979239421228991532051555943675wb,
+	   32783817256434357484609367438786815wb) != 0wb)
+    __builtin_abort ();
+  if (foo (542904728531209767665756029992981529373473101602268731408384wb,
+	   235447394450476261134537147263765988105wb)
+      != 235447394450476261116090403190056436489wb)
+    __builtin_abort ();
+  if (foo (542904728531209767665690117878483036079552477922114364506112wb,
+	   235447394450476261134537147263765988105wb)
+      != 169535279951982967195466723035689534217wb)
+    __builtin_abort ();
+  if (foo (542904728531209767665690117878483036079534031178040654954496wb,
+	   235447394450476261134537147263765988105wb)
+      != 169535279951982967177019978961979982601wb)
+    __builtin_abort ();
+  if (foo (542904728531209767665454670484032559818581851459155356811264wb,
+	   235447394450476261134537147263765988105wb)
+      != 169535279951982967359377407340447827474wb)
+    __builtin_abort ();
+#endif
+#if __BITINT_MAXWIDTH__ >= 512
+  if (bar (6703903964971298549787012499102923063739682910296196688861780721860882015036773488400937149083451713845015929093243025426876941405973284973216824503042048uwb,
+	   170141183460469231731687303715884105735uwb) != 19208uwb)
+    __builtin_abort ();
+#endif
+#if __BITINT_MAXWIDTH__ >= 256
+  if (baz (115792089237316195423570985008687907853269984665640564039457584007913129639926uwb,
+	   68056473384187692692674921486353642292uwb) != 6uwb)
+    __builtin_abort ();
+#endif
+  return 0;
+}

libgcc/libgcc2.c

@@ -1863,12 +1863,46 @@ __divmodbitint4 (UBILtype *q, SItype qprec,
 	      if (uv1 >= vv1)
 		{
 		  /* udiv_qrnnd doesn't support quotients which don't
-		     fit into UWtype, so subtract from uv1:uv0 vv1
-		     first.  */
-		  uv1 -= vv1 + __builtin_sub_overflow (uv0, vv1, &uv0);
-		  udiv_qrnnd (qhat, rhat, uv1, uv0, vv1);
-		  if (!__builtin_add_overflow (rhat, vv1, &rhat))
-		    goto again;
+		     fit into UWtype, while Knuth's algorithm originally
+		     uses a double-word by word to double-word division.
+		     Fortunately, the algorithm guarantees that uv1 <= vv1,
+		     because if uv1 > vv1, then even if v would have all
+		     bits in all words below vv1 set, the previous iteration
+		     would be supposed to use qhat larger by 1 and subtract
+		     v.  With uv1 == vv1 and uv0 >= vv1 the double-word
+		     qhat in Knuth's algorithm would be 1 in the upper word
+		     and 1 in the lower word, say for
+		     uv1 0x8000000000000000ULL
+		     uv0 0xffffffffffffffffULL
+		     vv1 0x8000000000000000ULL
+		     0x8000000000000000ffffffffffffffffuwb
+		     / 0x8000000000000000uwb == 0x10000000000000001uwb, and
+		     exactly like that also for any other value
+		     > 0x8000000000000000ULL in uv1 and vv1 and uv0 >= uv1.
+		     So we need to subtract one or at most two vv1s from
+		     uv1:uv0 (qhat because of that decreases by 1 or 2 and
+		     is then representable in UWtype) and need to increase
+		     rhat by vv1 once or twice because of that.  Now, if
+		     we need to subtract 2 vv1s, i.e. if
+		     uv1 == vv1 && uv0 >= vv1, then rhat (which is uv0 - vv1)
+		     + vv1 computation can't overflow, because it is equal
+		     to uv0 and therefore the original algorithm in that case
+		     performs goto again, but the second vv1 addition must
+		     overflow already because vv1 has msb set from the
+		     canonicalization.  */
+		  uv1 -= __builtin_sub_overflow (uv0, vv1, &uv0);
+		  if (uv1 >= vv1)
+		    {
+		      uv1 -= __builtin_sub_overflow (uv0, vv1, &uv0);
+		      udiv_qrnnd (qhat, rhat, uv1, uv0, vv1);
+		      rhat += 2 * vv1;
+		    }
+		  else
+		    {
+		      udiv_qrnnd (qhat, rhat, uv1, uv0, vv1);
+		      if (!__builtin_add_overflow (rhat, vv1, &rhat))
+			goto again;
+		    }
 		}
 	      else
 		{