# HG changeset patch
# User Gaël Bonithon <gael@xfce.org>
# Date 1645384872 -3600
#      Sun Feb 20 20:21:12 2022 +0100
# Node ID 5586228ac34a435b97b003cc7561fa1e5d7f287d
# Parent  496e5214c609c1927c21e1dd15d7b10445dc5e8d
idivide.m: Avoid floating-point representation issues (bug #61319)

* scripts/general/idivide.m:
We start from `z = x ./ y` which corresponds to the classical rounding, and
which we are sure is an integer value if `x` or `y` is an integer value. We
then add or subtract `1` if necessary according to the rounding rule and the
sign of `y`, keeping the starting integer type all along.

diff -r 496e5214c609 -r 5586228ac34a scripts/general/idivide.m
--- a/scripts/general/idivide.m	Sun Feb 13 12:47:36 2022 -0500
+++ b/scripts/general/idivide.m	Sun Feb 20 20:21:12 2022 +0100
@@ -84,33 +84,76 @@
     op = tolower (op);
   endif
 
-  if (isfloat (x))
-    if (isinteger (y))
-      typ = class (y);
-    else
-      error ("idivide: at least one input (X or Y) must be an integer type");
-    endif
-  elseif (isfloat (y))
-    if (isinteger (x))
-      typ = class (x);
-    else
-      error ("idivide: at least one input (X or Y) must be an integer type");
-    endif
-  else
-    typ = class (x);
-    if (! strcmp (typ, class (y)))
-      error ("idivide: integer type of X (%s) must match integer type of Y (%s)", typ, class (y));
-    endif
+  if (! isinteger (x) && ! isinteger (y))
+    error ("idivide: at least one input (X or Y) must be an integer type");
+  elseif (isinteger (x) && isinteger (y) && ! strcmp (class (x), class (y)))
+    error ("idivide: integer type of X (%s) must match integer type of Y (%s)",
+           class (x), class (y));
   endif
 
+  z = x ./ y;
   if (strcmp (op, "fix"))
-    z = cast (fix (double (x) ./ double (y)), typ);
+    ## The following is an optimized version of `z -= (z .* y > x) .* sign (y)`.
+    if (isscalar (y))
+      if (y > 0)
+        z -= (z * y > x);
+      else
+        z += (z * y > x);
+      endif
+    else
+      y_sel = (y > 0);
+      if (isscalar (x))
+        z(y_sel) -= (z(y_sel) .* y(y_sel) > x);
+        y_sel = ! y_sel;
+        z(y_sel) += (z(y_sel) .* y(y_sel) > x);
+      else
+        z(y_sel) -= (z(y_sel) .* y(y_sel) > x(y_sel));
+        y_sel = ! y_sel;
+        z(y_sel) += (z(y_sel) .* y(y_sel) > x(y_sel));
+      endif
+    endif
   elseif (strcmp (op, "round"))
-    z = x ./ y;
+    return;
   elseif (strcmp (op, "floor"))
-    z = cast (floor (double (x) ./ double (y)), typ);
+    ## The following is an optimized version of `z -= (z .* abs (y) > sign (y) .* x)`.
+    if (isscalar (y))
+      if (y > 0)
+        z -= (z * y > x);
+      else
+        z -= (z * y < x);
+      endif
+    else
+      y_sel = (y > 0);
+      if (isscalar (x))
+        z(y_sel) -= (z(y_sel) .* y(y_sel) > x);
+        y_sel = ! y_sel;
+        z(y_sel) -= (z(y_sel) .* y(y_sel) < x);
+      else
+        z(y_sel) -= (z(y_sel) .* y(y_sel) > x(y_sel));
+        y_sel = ! y_sel;
+        z(y_sel) -= (z(y_sel) .* y(y_sel) < x(y_sel));
+      endif
+    endif
   elseif (strcmp (op, "ceil"))
-    z = cast (ceil (double (x) ./ double (y)), typ);
+    ## The following is an optimized version of `z += (z .* abs (y) < sign (y) .* x)`.
+    if (isscalar (y))
+      if (y > 0)
+        z += (z * y < x);
+      else
+        z += (z * y > x);
+      endif
+    else
+      y_sel = (y > 0);
+      if (isscalar (x))
+        z(y_sel) += (z(y_sel) .* y(y_sel) < x);
+        y_sel = ! y_sel;
+        z(y_sel) += (z(y_sel) .* y(y_sel) > x);
+      else
+        z(y_sel) += (z(y_sel) .* y(y_sel) < x(y_sel));
+        y_sel = ! y_sel;
+        z(y_sel) += (z(y_sel) .* y(y_sel) > x(y_sel));
+      endif
+    endif
   else
     error ('idivide: unrecognized rounding type "%s"', op);
   endif
@@ -118,11 +161,13 @@
 endfunction
 
 
-%!shared a, af, b, bf
+%!shared a, af, b, bf, c, d
 %! a = int8 (3);
 %! af = 3;
 %! b = int8 ([-4, 4]);
 %! bf = [-4, 4];
+%! c = int64 (4e16);
+%! d = int64 ([-2e8, 2e8]);
 
 %!assert (idivide (a, b), int8 ([0, 0]))
 %!assert (idivide (a, b, "floor"), int8 ([-1, 0]))
@@ -139,6 +184,16 @@
 %!assert (idivide (a, bf, "ceil"), int8 ([0, 1]))
 %!assert (idivide (a, bf, "round"), int8 ([-1, 1]))
 
+%!assert (idivide (c, d + int64 (1)), d + int64 ([-1, -1]))
+%!assert (idivide (c, d + int64 (1), "floor"), d + int64 ([-2, -1]))
+%!assert (idivide (c, d + int64 (1), "ceil"), d + int64 ([-1, 0]))
+%!assert (idivide (c, d + int64 (1), "round"), d + int64 ([-1, -1]))
+
+%!assert (idivide (c + int64 (1), d), d)
+%!assert (idivide (c + int64 (1), d, "floor"), d + int64 ([-1, 0]))
+%!assert (idivide (c + int64 (1), d, "ceil"), d + int64 ([0, 1]))
+%!assert (idivide (c + int64 (1), d, "round"), d)
+
 ## Test input validation
 %!error idivide (uint8 (1))
 %!error idivide (uint8 (1), 2, 3)
