From fbb569315a291d2d5b32ad0fdaf0c42da9f5e93b Mon Sep 17 00:00:00 2001 From: Jakub Jelinek Date: Fri, 2 Feb 2024 22:14:33 +0100 Subject: [PATCH] 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 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 | 26 +++++++++++++ gcc/testsuite/gcc.dg/torture/bitint-55.c | 66 ++++++++++++++++++++++++++++++++ libgcc/libgcc2.c | 46 +++++++++++++++++++--- 3 files changed, 132 insertions(+), 6 deletions(-) create mode 100644 gcc/testsuite/gcc.dg/torture/bitint-53.c create mode 100644 gcc/testsuite/gcc.dg/torture/bitint-55.c diff --git a/gcc/testsuite/gcc.dg/torture/bitint-53.c b/gcc/testsuite/gcc.dg/torture/bitint-53.c new file mode 100644 index 00000000000..8ead40e450b --- /dev/null +++ b/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; +} diff --git a/gcc/testsuite/gcc.dg/torture/bitint-55.c b/gcc/testsuite/gcc.dg/torture/bitint-55.c new file mode 100644 index 00000000000..7d4bff3fd12 --- /dev/null +++ b/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; +} diff --git a/libgcc/libgcc2.c b/libgcc/libgcc2.c index e04d169f888..ef46153731f 100644 --- a/libgcc/libgcc2.c +++ b/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 { -- 2.11.4.GIT