Fix fast_integer_divide for signed numerator with denominator == 1

src/FastIntegerDivide.cpp

@@ -168,6 +168,11 @@ Expr fast_integer_divide_impl(Expr numerator, Expr denominator, bool round_to_ze
 
     Type wide = t.widen();
 
+    // Preserve the original numerator: the signed non-round-to-zero branch
+    // below rewrites `numerator` in place, but the denominator == 1 fixup
+    // needs the unmodified value.
+    Expr original_numerator = numerator;
+
     Expr result;
     if (t.is_uint()) {
         Expr mul, shift = shift_for_denominator(denominator);
@@ -294,8 +299,12 @@ Expr fast_integer_divide_impl(Expr numerator, Expr denominator, bool round_to_ze
         result -= xsign;
     }
 
-    // The tables don't work for denominator == 1
-    result = select(std::move(denominator) == 1, std::move(numerator), result);
+    // The tables don't work for denominator == 1. Denominator == 0 is
+    // defined to give a result of 0, matching Halide's ordinary division
+    // operator.
+    Expr is_zero = denominator == 0;
+    Expr is_one = std::move(denominator) == 1;
+    result = select(is_zero, make_zero(t), select(is_one, std::move(original_numerator), std::move(result)));
 
     internal_assert(result.type() == t);
 
@@ -314,7 +323,12 @@ Expr fast_integer_divide(const Expr &numerator, const Expr &denominator) {
 
 Expr fast_integer_modulo(const Expr &numerator, const Expr &denominator) {
     Expr ratio = fast_integer_divide(numerator, denominator);
-    return numerator - ratio * denominator;
+    Expr result = numerator - ratio * denominator;
+    // Denominator == 0 gives a result of 0, matching Halide's ordinary
+    // modulo operator. fast_integer_divide already returns 0 for a zero
+    // denominator, but multiplying that back by denominator (also 0)
+    // would otherwise leave the numerator unchanged instead of zeroing it.
+    return select(denominator == 0, make_zero(numerator.type()), std::move(result));
 }
 
 }  // namespace Halide

src/FastIntegerDivide.h

@@ -22,9 +22,8 @@ namespace Halide {
  * bit vectors. For 32-bit vectors on x86 you're better off using
  * native integer division.
  *
- * Also, this routine treats division by zero as division by
- * 256. I.e. it interprets the uint8 divisor as a number from 1 to 256
- * inclusive.
+ * Like Halide's ordinary division operator, a denominator of zero gives
+ * a result of zero.
  */
 Expr fast_integer_divide(const Expr &numerator, const Expr &denominator);
 
@@ -33,7 +32,9 @@ Expr fast_integer_divide(const Expr &numerator, const Expr &denominator);
 Expr fast_integer_divide_round_to_zero(const Expr &numerator, const Expr &denominator);
 
 /** Use the fast integer division tables to implement a modulo
- * operation via the Euclidean identity: a%b = a - (a/b)*b
+ * operation via the Euclidean identity: a%b = a - (a/b)*b. Like
+ * Halide's ordinary modulo operator, a denominator of zero gives a
+ * result of zero.
  */
 Expr fast_integer_modulo(const Expr &numerator, const Expr &denominator);
 

test/performance/const_division.cpp

@@ -115,6 +115,72 @@ bool test(int w, bool div, bool round_to_zero) {
         }
     }
 
+    // Exhaustively check fast_integer_divide{,_round_to_zero}/fast_integer_modulo
+    // against Halide's native division/modulo across the complete space of
+    // uint8 denominators, 0 to 255 inclusive. A denominator of 0 gives a
+    // result of 0, matching Halide's ordinary / and %.
+    //
+    // The reference is computed in a signed type twice as wide as T: that's
+    // wide enough to hold both T's full range and the denominator's true
+    // (always non-negative) value, so for 8-bit signed T, a denominator
+    // above 127 isn't reinterpreted as negative the way casting it straight
+    // into T would. It also keeps native div_round_to_zero -- which, unlike
+    // / and %, isn't guaranteed safe for a zero divisor or for INT_MIN / -1
+    // -- away from a literal zero divisor.
+    {
+        const int denom_extent = 256;
+        Buffer<T> full_input(w, denom_extent);
+        for (int y = 0; y < denom_extent; y++) {
+            for (int x = 0; x < w; x++) {
+                full_input(x, y) = (T)(uint32_t)rng();
+                if (round_to_zero && full_input(x, y) == 0) {
+                    full_input(x, y) = 1;
+                }
+            }
+        }
+
+        Type wide = Int((int)bits * 2);
+        Expr wide_num = cast(wide, full_input(x, y));
+        Expr wide_den = cast(wide, y);
+
+        Func correct_full, fast_full;
+        if (div) {
+            if (round_to_zero) {
+                // Keep the actual hardware division away from a divisor of
+                // zero (unlike / and %, div_round_to_zero doesn't protect
+                // against this itself), and patch the result afterwards.
+                Expr safe_den = select(wide_den == 0, cast(wide, 1), wide_den);
+                correct_full(x, y) = cast<T>(select(wide_den == 0, cast(wide, 0),
+                                                    div_round_to_zero(wide_num, safe_den)));
+                fast_full(x, y) = Halide::fast_integer_divide_round_to_zero(full_input(x, y), cast<uint8_t>(y));
+            } else {
+                correct_full(x, y) = cast<T>(wide_num / wide_den);
+                fast_full(x, y) = Halide::fast_integer_divide(full_input(x, y), cast<uint8_t>(y));
+            }
+        } else {
+            correct_full(x, y) = cast<T>(wide_num % wide_den);
+            fast_full(x, y) = Halide::fast_integer_modulo(full_input(x, y), cast<uint8_t>(y));
+        }
+
+        correct_full.compile_jit(t);
+        fast_full.compile_jit(t);
+
+        Buffer<T> correct_full_result = correct_full.realize({w, denom_extent});
+        Buffer<T> fast_full_result = fast_full.realize({w, denom_extent});
+
+        for (int y = 0; y < denom_extent; y++) {
+            for (int x = 0; x < w; x++) {
+                if (fast_full_result(x, y) != correct_full_result(x, y)) {
+                    printf("fast_integer_divide/modulo(%lld, %d) = %lld instead of %lld\n",
+                           (long long int)full_input(x, y), y,
+                           (long long int)fast_full_result(x, y),
+                           (long long int)correct_full_result(x, y));
+                    return false;
+                }
+            }
+        }
+    }
+
     return true;
 }
 

test/performance/fast_inverse.cpp

@@ -12,6 +12,17 @@ int main(int argc, char **argv) {
         return 0;
     }
 
+    if (target.arch == Target::X86 && target.bits == 32) {
+        // The 32-bit x86 ABI only exposes 8 XMM registers, but this
+        // benchmark's vectorized recurrence needs 9 live values at once
+        // (8 independent lane groups plus the shared "+1" constant). That
+        // forces a stack spill/reload into the tightly-coupled loop-carried
+        // dependency chain every iteration, adding enough latency to erase
+        // rcpps's tiny edge over divps.
+        printf("[SKIP] x86-32 doesn't have enough XMM registers to avoid spilling.\n");
+        return 0;
+    }
+
     if (target.arch == Target::ARM &&
         target.os == Target::OSX) {
         // vrecpe, vrecps, fmul have inverse throughputs of 1, 0.25, 0.25