diff --git a/dace/codegen/cppunparse.py b/dace/codegen/cppunparse.py index f80a3f841c..7f95cdfda4 100644 --- a/dace/codegen/cppunparse.py +++ b/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 diff --git a/dace/runtime/include/dace/math.h b/dace/runtime/include/dace/math.h index 629db7734e..d0482378ce 100644 --- a/dace/runtime/include/dace/math.h +++ b/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 ::value && std::is_signed::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 ::value && std::is_unsigned::value>* = nullptr> static DACE_CONSTEXPR DACE_HDFI T int_floor_ni(const T& numerator, const T& denominator) { @@ -233,6 +238,14 @@ template::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::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)numerator, (T)denominator); +} template static DACE_CONSTEXPR DACE_HDFI std::complex py_floor(const std::complex& numerator, const std::complex& denominator) { std::complex 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::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)numerator, (T)denominator); +} + // Computes C/C++ modulus (operator % and fmod) template::value>* = nullptr> static DACE_CONSTEXPR DACE_HDFI T cpp_mod(const T& numerator, const T& denominator) { diff --git a/tests/cppunparse_test.py b/tests/cppunparse_test.py index ebe8548ebc..a0adc6e87b 100644 --- a/tests/cppunparse_test.py +++ b/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)') diff --git a/tests/python_division_semantics_test.py b/tests/python_division_semantics_test.py new file mode 100644 index 0000000000..b398c1660b --- /dev/null +++ b/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 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(')