PR64182: Fix rounding division and modulus

[Richard Sandiford]

As pointed out in PR 64182, wide-int rounded division gets the
ties-away-from-zero case wrong for odd-numbered dividends, while
double_int gets the unsigned case wrong by unconditionally treating
a dividend or remainder with the top bit set as negative.  As Jakub
says, the test used in double_int might also have overflow problems.

This patch uses:

   abs (remainder) >= abs (dividend) - abs (remainder)

for both wide-int and double_int and fixes the unsigned case in double_int.
I didn't know how to test the double_int change using input code so
resorted to doing some double_int arithmetic at the start of main.

Thanks to Joseph for the testcase.

Tested on x86_64-linux-gnu.  OK to install?

Thanks,
Richard


gcc/
        PR middle-end/64182
        * wide-int.h (wi::div_round, wi::mod_round): Fix rounding of tied
        cases.
        * double-int.c (div_and_round_double): Fix handling of unsigned
        cases.  Use same rounding approach as wide-int.h.

gc/testsuite/
2014-xx-xx  Joseph Myers  <jos...@codesourcery.com>

        PR middle-end/64182
        * gnat.dg/round_div.adb: New test.

Index: gcc/double-int.c
===================================================================
--- gcc/double-int.c    2014-12-11 10:45:44.430786435 +0000
+++ gcc/double-int.c    2014-12-11 10:46:10.570461030 +0000
@@ -569,24 +569,23 @@ div_and_round_double (unsigned code, int
       {
        unsigned HOST_WIDE_INT labs_rem = *lrem;
        HOST_WIDE_INT habs_rem = *hrem;
-       unsigned HOST_WIDE_INT labs_den = lden, ltwice;
-       HOST_WIDE_INT habs_den = hden, htwice;
+       unsigned HOST_WIDE_INT labs_den = lden, lnegabs_rem, ldiff;
+       HOST_WIDE_INT habs_den = hden, hnegabs_rem, hdiff;
 
        /* Get absolute values.  */
-       if (*hrem < 0)
+       if (!uns && *hrem < 0)
          neg_double (*lrem, *hrem, &labs_rem, &habs_rem);
-       if (hden < 0)
+       if (!uns && hden < 0)
          neg_double (lden, hden, &labs_den, &habs_den);
 
-       /* If (2 * abs (lrem) >= abs (lden)), adjust the quotient.  */
-       mul_double ((HOST_WIDE_INT) 2, (HOST_WIDE_INT) 0,
-                   labs_rem, habs_rem, &ltwice, &htwice);
+       /* If abs(rem) >= abs(den) - abs(rem), adjust the quotient.  */
+       neg_double (labs_rem, habs_rem, &lnegabs_rem, &hnegabs_rem);
+       add_double (labs_den, habs_den, lnegabs_rem, hnegabs_rem,
+                   &ldiff, &hdiff);
 
-       if (((unsigned HOST_WIDE_INT) habs_den
-            < (unsigned HOST_WIDE_INT) htwice)
-           || (((unsigned HOST_WIDE_INT) habs_den
-                == (unsigned HOST_WIDE_INT) htwice)
-               && (labs_den <= ltwice)))
+       if (((unsigned HOST_WIDE_INT) habs_rem
+            > (unsigned HOST_WIDE_INT) hdiff)
+           || (habs_rem == hdiff && labs_rem >= ldiff))
          {
            if (quo_neg)
              /* quo = quo - 1;  */
Index: gcc/testsuite/gnat.dg/round_div.adb
===================================================================
--- /dev/null   2014-11-19 08:41:51.310561007 +0000
+++ gcc/testsuite/gnat.dg/round_div.adb 2014-12-11 10:46:10.570461030 +0000
@@ -0,0 +1,17 @@
+-- { dg-do run }
+-- { dg-options "-O3" }
+procedure Round_Div is
+   type Fixed is delta 1.0 range -2147483648.0 .. 2147483647.0;
+   A : Fixed := 1.0;
+   B : Fixed := 3.0;
+   C : Integer;
+   function Divide (X, Y : Fixed) return Integer is
+   begin
+      return Integer (X / Y);
+   end;
+begin
+   C := Divide (A, B);
+   if C /= 0 then
+      raise Program_Error;
+   end if;
+end Round_Div;
Index: gcc/wide-int.h
===================================================================
--- gcc/wide-int.h      2014-12-11 10:45:44.434786385 +0000
+++ gcc/wide-int.h      2014-12-11 10:46:10.570461030 +0000
@@ -2616,8 +2616,8 @@ wi::div_round (const T1 &x, const T2 &y,
     {
       if (sgn == SIGNED)
        {
-         if (wi::ges_p (wi::abs (remainder),
-                        wi::lrshift (wi::abs (y), 1)))
+         WI_BINARY_RESULT (T1, T2) abs_remainder = wi::abs (remainder);
+         if (wi::geu_p (abs_remainder, wi::abs (y) - abs_remainder))
            {
              if (wi::neg_p (x, sgn) != wi::neg_p (y, sgn))
                return quotient - 1;
@@ -2627,7 +2627,7 @@ wi::div_round (const T1 &x, const T2 &y,
        }
       else
        {
-         if (wi::geu_p (remainder, wi::lrshift (y, 1)))
+         if (wi::geu_p (remainder, y - remainder))
            return quotient + 1;
        }
     }
@@ -2784,8 +2784,8 @@ wi::mod_round (const T1 &x, const T2 &y,
     {
       if (sgn == SIGNED)
        {
-         if (wi::ges_p (wi::abs (remainder),
-                        wi::lrshift (wi::abs (y), 1)))
+         WI_BINARY_RESULT (T1, T2) abs_remainder = wi::abs (remainder);
+         if (wi::geu_p (abs_remainder, wi::abs (y) - abs_remainder))
            {
              if (wi::neg_p (x, sgn) != wi::neg_p (y, sgn))
                return remainder + y;
@@ -2795,7 +2795,7 @@ wi::mod_round (const T1 &x, const T2 &y,
        }
       else
        {
-         if (wi::geu_p (remainder, wi::lrshift (y, 1)))
+         if (wi::geu_p (remainder, y - remainder))
            return remainder - y;
        }
     }