[To Discuss] Python floor division and modulo lower to C semantics

dace/codegen/cppunparse.py

@@ -851,14 +851,24 @@ def _UnaryOp(self, t):
         "Sub": "-",
         "Mult": "*",
         "Div": "/",
-        "Mod": "%",
         "LShift": "<<",
         "RShift": ">>",
         "BitOr": "|",
         "BitXor": "^",
         "BitAnd": "&"
     }
-    funcops = {"FloorDiv": (" /", "dace::math::ifloor"), "MatMult": (",", "dace::gemm")}
+    # ``//`` and ``%`` are Python's, hence numpy's: the quotient rounds toward negative infinity and
+    # the remainder therefore takes the divisor's sign (``-32 // 7 == -5``, ``-32 % 7 == 3``). C
+    # rounds toward zero and gives the remainder the dividend's sign, so neither can be written
+    # infix. ``ifloor(a / b)`` was wrong on integers, where ``a / b`` has already truncated and
+    # flooring an integer changes nothing; a bare ``%`` was wrong the same way on integers and does
+    # not compile at all on floats, which C has no ``%`` for. ``py_floor`` / ``py_mod`` dispatch on
+    # the operand type and answer for every one of them.
+    funcops = {
+        "FloorDiv": (",", "py_floor"),
+        "Mod": (",", "py_mod"),
+        "MatMult": (",", "dace::gemm"),
+    }
 
     def _BinOp(self, t):
         # Operations that require a function call

dace/runtime/include/dace/math.h

@@ -212,9 +212,14 @@ DACE_CONSTEXPR __device__ __forceinline__ dace::float16 max(const T& a, const da
 // https://stackoverflow.com/a/39304947
 template <typename T, std::enable_if_t<std::is_integral<T>::value && std::is_signed<T>::value>* = nullptr>
 static DACE_CONSTEXPR DACE_HDFI T int_floor_ni(const T& numerator, const T& denominator) {
-    auto divresult = std::div(numerator, denominator);
-    T corr = (divresult.rem != 0 && ((divresult.rem < 0) != (denominator < 0)));
-    return (T)divresult.quot - corr;
+    // ``/`` and ``%`` rather than ``std::div``: the latter is host-only, and nvcc answers a call to
+    // it from device code with a warning rather than an error, then removes the guarded region that
+    // holds the call -- the kernel launches and stores nothing. C++11 pins ``/`` to truncation
+    // toward zero, so the correction below is exact, and this form is constexpr where std::div is not.
+    const T quotient = numerator / denominator;
+    const T remainder = numerator % denominator;
+    const T corr = (remainder != 0 && ((remainder < 0) != (denominator < 0)));
+    return quotient - corr;
 }
 template <typename T, std::enable_if_t<std::is_integral<T>::value && std::is_unsigned<T>::value>* = nullptr>
 static DACE_CONSTEXPR DACE_HDFI T int_floor_ni(const T& numerator, const T& denominator) {
@@ -233,6 +238,14 @@ template<typename T, std::enable_if_t<!std::is_integral<T>::value && std::is_flo
 static DACE_CONSTEXPR DACE_HDFI T py_floor(const T& numerator, const T& denominator) {
     return (T)std::floor(numerator / denominator);
 }
+// Mixed-operand-type overload (e.g. ``a // 7``, where ``a`` is int64_t and the literal ``7`` is
+// int): promote both operands to their common arithmetic type and delegate, since the single-type
+// templates above cannot deduce T from two different types.
+template<typename T1, typename T2, std::enable_if_t<!std::is_same<T1, T2>::value>* = nullptr>
+static DACE_CONSTEXPR DACE_HDFI auto py_floor(const T1& numerator, const T2& denominator) -> decltype(numerator + denominator) {
+    using T = decltype(numerator + denominator);
+    return py_floor<T>((T)numerator, (T)denominator);
+}
 template<typename T>
 static DACE_CONSTEXPR DACE_HDFI std::complex<T> py_floor(const std::complex<T>& numerator, const std::complex<T>& denominator) {
     std::complex<T> quotient = numerator / denominator;
@@ -261,6 +274,13 @@ static DACE_CONSTEXPR DACE_HDFI T py_mod(const T& numerator, const T& denominato
     return (T)(numerator - quotient * denominator);
 }
 
+// Mixed-operand-type overload, as for py_floor above.
+template<typename T1, typename T2, std::enable_if_t<!std::is_same<T1, T2>::value>* = nullptr>
+static DACE_CONSTEXPR DACE_HDFI auto py_mod(const T1& numerator, const T2& denominator) -> decltype(numerator + denominator) {
+    using T = decltype(numerator + denominator);
+    return py_mod<T>((T)numerator, (T)denominator);
+}
+
 // Computes C/C++ modulus (operator % and fmod)
 template<typename T, std::enable_if_t<std::is_integral<T>::value>* = nullptr>
 static DACE_CONSTEXPR DACE_HDFI T cpp_mod(const T& numerator, const T& denominator) {

tests/cppunparse_test.py

@@ -60,7 +60,7 @@ def test():
     auto result = 0;
     while (((i < woo) && (i > 0))) {
         for (auto j : range(i)) {
-            result += dace::math::pow(dace::math::ifloor(2 / 1), j);
+            result += dace::math::pow(py_floor(2, 1), j);
         }
     }
     return result;
@@ -88,7 +88,7 @@ def lfunc() -> None:
 
     # Operations (augmented assignment)
     success &= _test_py2cpp('l *= 3; l //= 8', """l *= 3;
-l = dace::math::ifloor(l / 8);""")
+l = py_floor(l, 8);""")
 
     success &= _test_pyexpr2cpp('a << 3', '(a << 3)')
 

tests/python_division_semantics_test.py

@@ -0,0 +1,85 @@
+# Copyright 2019-2026 ETH Zurich and the DaCe authors. All rights reserved.
+""" Tests that ``//`` and ``%`` in a tasklet follow Python (hence NumPy) semantics. """
+import itertools
+
+import numpy as np
+import pytest
+
+import dace
+
+N = dace.symbol('N', dtype=dace.int64)
+
+# Both signs, and divisors that do and do not divide evenly: the disagreement between Python and C
+# needs a nonzero remainder and operands of opposite sign.
+VALUES = (-32, -7, -3, -1, 1, 3, 7, 32)
+
+DTYPES = (('int32', dace.int32, np.int32), ('int64', dace.int64, np.int64), ('float32', dace.float32, np.float32),
+          ('float64', dace.float64, np.float64))
+
+TARGETS = (pytest.param('host'), pytest.param('device', marks=pytest.mark.gpu))
+
+
+def division_program(op: str, dtype):
+    """ ``out = a <op> b`` elementwise, as a map so that the GPU lowering is a real kernel. """
+    if op == '//':
+
+        @dace.program
+        def prog(a: dtype[N], b: dtype[N], out: dtype[N]):
+            for i in dace.map[0:N]:
+                out[i] = a[i] // b[i]
+    else:
+
+        @dace.program
+        def prog(a: dtype[N], b: dtype[N], out: dtype[N]):
+            for i in dace.map[0:N]:
+                out[i] = a[i] % b[i]
+
+    return prog
+
+
+def operand_pairs(nptype):
+    pairs = list(itertools.product(VALUES, VALUES))
+    return (np.array([x for x, _ in pairs], dtype=nptype), np.array([y for _, y in pairs], dtype=nptype))
+
+
+@pytest.mark.parametrize('op', ('//', '%'))
+@pytest.mark.parametrize('name,dtype,nptype', DTYPES)
+@pytest.mark.parametrize('target', TARGETS)
+def test_floor_division_and_modulo_agree_with_numpy(op, name, dtype, nptype, target):
+    """ Python rounds the quotient toward negative infinity, so the remainder takes the divisor's
+        sign: ``-32 // 7 == -5`` and ``-32 % 7 == 3``. C rounds toward zero and answers ``-4`` for
+        both.
+
+        The GPU half of the table is not redundant: the correction term is a branch, and that branch
+        used to hold a call to host-only ``std::div``. nvcc reports such a call with a warning rather
+        than an error and then removes the region around it, so the kernel launched, reported
+        success, and stored nothing.
+    """
+    a, b = operand_pairs(nptype)
+    expected = (a // b) if op == '//' else (a % b)
+
+    sdfg = division_program(op, dtype).to_sdfg()
+    if target == 'device':
+        sdfg.apply_gpu_transformations()
+    out = np.zeros(a.shape, dtype=nptype)
+    sdfg(a=a, b=b, out=out, N=a.shape[0])
+
+    disagreements = [(a[i], b[i], out[i], expected[i]) for i in range(a.shape[0]) if out[i] != expected[i]]
+    assert not disagreements, f'{op} on {name}: {disagreements[:4]}'
+
+
+@pytest.mark.parametrize('op,call', (('//', 'py_floor('), ('%', 'py_mod(')))
+def test_neither_operator_is_emitted_infix(op, call):
+    """ The numeric table above catches an infix ``%`` on integers, but not on floats, where the
+        emitted line does not compile and there is no number left to compare.
+    """
+    code = division_program(op, dace.int64).to_sdfg().generate_code()[0].clean_code
+    assert call in code, f'{op} lowered without {call}, so it lowered infix'
+
+
+if __name__ == '__main__':
+    for operator in ('//', '%'):
+        for dtype_name, dace_type, numpy_type in DTYPES:
+            test_floor_division_and_modulo_agree_with_numpy(operator, dtype_name, dace_type, numpy_type, 'host')
+    test_neither_operator_is_emitted_infix('//', 'py_floor(')
+    test_neither_operator_is_emitted_infix('%', 'py_mod(')