diff --git a/graphblas/core/operator/udt_utils.py b/graphblas/core/operator/udt_utils.py index 0ce992cfc..6625483ce 100644 --- a/graphblas/core/operator/udt_utils.py +++ b/graphblas/core/operator/udt_utils.py @@ -128,8 +128,39 @@ def _compile_codegen(src, *, func_name, source_label, extra_ns=None): # (no ordering for min/max, no integer-mod for floordiv). _OPS_NOT_FOR_COMPLEX = frozenset({"min", "max", "floordiv"}) -# C operator equivalents for JIT code generation -_C_INFIX_OPS = {"+": "+", "-": "-", "*": "*", "/": "/", "//": "/"} +# ``_MIN`` spellings from , keyed by itemsize. Signed integer +# division traps on ``MIN / -1``; the generated C tests for it explicitly. +_C_INT_MIN = {1: "INT8_MIN", 2: "INT16_MIN", 4: "INT32_MIN", 8: "INT64_MIN"} + +# Same constants for the generated Python source, spelled ``-MAX - 1`` rather +# than as a direct literal that reads as the negation of an out-of-range value. +# MAINT 2026-07-30: Numba 0.65 types both spellings as ``Literal[int]``, so this +# can likely be simplified; re-check on the oldest supported Numba first. +_PY_INT_MIN = { + 1: "(-127 - 1)", + 2: "(-32767 - 1)", + 4: "(-2147483647 - 1)", + 8: "(-9223372036854775807 - 1)", +} + + +def _is_float(np_dtype): + return np_dtype.kind == "f" + + +def _is_signed_int(np_dtype): + return np_dtype.kind == "i" + + +def _is_int(np_dtype): + """True for signed and unsigned integers, but not bool. + + Bool is excluded from the integer division guards: ``_Bool`` has no range + to overflow and C converts a double to it by comparing against zero, so + both execution paths already agree without help. + """ + return np_dtype.kind in ("i", "u") + # Vanilla strips GxB callables but keeps GxB constants, so the bare # ``hasattr`` would lie; gate on the backend too. @@ -539,12 +570,64 @@ def _expr_binary(py_op, x_expr, y_expr, x_dtype, y_dtype): """Python-source builder; sibling of :func:`_c_expr_binary` for JIT C. ``x_dtype`` and ``y_dtype`` are each operand's numpy dtype at this - leaf, which ``_check_udt_pair`` allows to differ. Only ``min`` and - ``max`` consult them, to tell a floating-point leaf from an integer - one; the rest are type-agnostic. + leaf, which ``_check_udt_pair`` allows to differ. Four ops consult + them: + + - ``min`` / ``max`` follow SuiteSparse's own ``GrB_MIN_FP64``, which + is C99 ``fmin`` and ignores a NaN operand rather than ordering or + propagating it; see :func:`_minmax_expr`. Only floating-point + leaves carry the NaN and signed-zero rules, so the dtype decides + which of the two forms is emitted. + - ``//`` on two signed integers. Numba deliberately returns 0 for + ``INT_MIN // -1`` to dodge the SIGFPE that x86 raises on the + unrepresentable quotient, while numpy wraps to ``INT_MIN``. Since + ``a // -1`` is exactly ``-a``, routing that divisor through + negation (which wraps) reaches numpy's answer without the trap. + - ``/`` on two integers. The quotient is a float that then has to fit + back into the integer field; when it doesn't, the conversion is + undefined in both C and LLVM, and the two disagree in ways that + vary with the field's width. Ruling out the two ways a division of + two same-signedness integers can leave the field's range (a zero + divisor, and ``INT_MIN / -1``) is what makes the two paths agree. + Note this gives ``x / 0 == 0``, which agrees with + ``np.floor_divide``; ``np.true_divide`` produces an infinity whose + cast back to an integer is undefined in numpy too. + + Two range escapes are left, and both move the two paths together + rather than apart. A 64-bit field can leave the range through the + ``(double)`` conversion itself (``(2**63 - 1) / 1`` rounds up to + ``2**63``), and there the answer is whatever the hardware does: + measured saturating to ``INT64_MAX`` on arm64. Operands of mixed + signedness escape through a negative divisor, which the + ``INT_MIN`` guard, keyed on the left operand, doesn't see: + ``uint8(200) / int8(-1)`` is -200.0 before the cast back. That one + is unguarded rather than known-wrong (it lands on 56, which is what + ``np.float64(-200).astype(np.uint8)`` gives), and a mixed pair gets + no JIT kernel at all, so there is no second path for it to + disagree with. + - ``/`` where either side is complex. Numba's complex division raises + ``ZeroDivisionError`` unconditionally, outside the error model's + control, so a zero divisor left the element unwritten while the C + kernel returned numpy's infinities. Spelling out the zero case here + makes both paths match numpy. """ if py_op in _FUNC_BINARY_OPS: return _minmax_expr(py_op, x_expr, y_expr, x_dtype, y_dtype) + if py_op == "//" and _is_signed_int(x_dtype) and _is_signed_int(y_dtype): + return f"(-{x_expr} if {y_expr} == -1 else {x_expr} // {y_expr})" + if py_op == "/" and (x_dtype.kind == "c" or y_dtype.kind == "c"): + return ( + f"({x_expr} / {y_expr} if {y_expr} != 0 " + f"else complex({x_expr}.real / 0.0, {x_expr}.imag / 0.0))" + ) + if py_op == "/" and _is_int(x_dtype) and _is_int(y_dtype): + quotient = f"{x_expr} / {y_expr}" + if _is_signed_int(x_dtype): + type_min = _PY_INT_MIN[x_dtype.itemsize] + quotient = ( + f"({x_expr} if ({x_expr} == {type_min} and {y_expr} == -1) else {quotient})" + ) + return f"(0 if {y_expr} == 0 else {quotient})" return f"{x_expr} {py_op} {y_expr}" def _minmax_expr(py_op, x_expr, y_expr, x_dtype, y_dtype): @@ -571,7 +654,7 @@ def _minmax_expr(py_op, x_expr, y_expr, x_dtype, y_dtype): over an int record still produces float results. """ cmp_op = "<" if py_op == "min" else ">" - if x_dtype.kind != "f" and y_dtype.kind != "f": + if not _is_float(x_dtype) and not _is_float(y_dtype): return f"({x_expr} if {x_expr} {cmp_op} {y_expr} else {y_expr})" # ``y != y`` is the NaN test on the right operand: when it holds, the # left one wins whatever it is, which is fmin's "ignore the NaN" rule @@ -855,6 +938,42 @@ def compile_udt_unary_wrapper(op_name, py_op, dtype): # JIT C code generators below. +def _c_assign_binary(py_op, target, lhs, rhs, field_dtype): + """Return a C *statement* assigning ``lhs py_op rhs`` to ``target``. + + Almost every op is a single expression, but floating-point ``//`` needs + temporaries (see :func:`_c_float_floordiv_stmt`), so this is the entry + point the kernel builder uses rather than :func:`_c_expr_binary`. + """ + if py_op == "//" and field_dtype.kind == "f": + return _c_float_floordiv_stmt(target, lhs, rhs, field_dtype) + return f"{target} = {_c_expr_binary(py_op, lhs, rhs, field_dtype)} ;" + + +def _c_expr_binary(py_op, lhs, rhs, field_dtype): + """Return a C expression for a binary op: e.g., ``(x->a) + (y->a)``. + + ``field_dtype`` is the numpy dtype of the element being written. Four of + the ops need it because the C spelling that looks obvious disagrees with + what SuiteSparse's own operators compute: ``/`` is integer division in C + but true division in Python, ``//`` is neither, and ``min`` / ``max`` have + to ignore NaN the way ``fmin`` does. It is required rather than optional + so a new caller can't silently get the wrong kernel. + + Floating-point ``//`` has no single-expression form; go through + :func:`_c_assign_binary`. + """ + if py_op in _FUNC_BINARY_OPS: + return _c_minmax_expr(py_op, lhs, rhs, field_dtype) + if py_op == "/": + return _c_truediv_expr(lhs, rhs, field_dtype) + if py_op == "//": + return _c_int_floordiv_expr(lhs, rhs, field_dtype) + # Everything left (``+``, ``-``, ``*``) spells the same in C as in Python + # and needs no dtype-specific handling. + return f"({lhs}) {py_op} ({rhs})" + + def _c_minmax_expr(py_op, lhs, rhs, field_dtype): """Return a C expression for ``min`` / ``max`` matching SuiteSparse. @@ -863,7 +982,7 @@ def _c_minmax_expr(py_op, lhs, rhs, field_dtype): different things depending on the dtype it is typed for. Calling ``fmin`` itself is the way to be sure of that, tie-break on signed zeros included. ```` reaches the JIT kernel through ``GraphBLAS.h``, the same - route ``_c_floordiv_expr`` relies on for ``floor``. + route ``_c_float_floordiv_stmt`` relies on for ``floor`` and ``fmod``. ``float`` fields need ``fminf``: passing them to ``fmin`` would compute in double and round on the way back, which costs nothing in accuracy for @@ -873,12 +992,8 @@ def _c_minmax_expr(py_op, lhs, rhs, field_dtype): exact for them. ``_minmax_expr`` emits the matching Python source; the two must agree, because SuiteSparse picks between the JIT kernel and the Numba cfunc on its own. - - ``field_dtype`` is None only if a caller omits it. Falling back to the - comparison keeps that case compiling; ``_make_jit_c_definition`` always - passes the real leaf dtype. """ - if field_dtype is not None and field_dtype.kind == "f": + if _is_float(field_dtype): fn = "fmin" if py_op == "min" else "fmax" suffix = "f" if field_dtype.itemsize == 4 else "" return f"{fn}{suffix}(({lhs}), ({rhs}))" @@ -886,52 +1001,124 @@ def _c_minmax_expr(py_op, lhs, rhs, field_dtype): return f"(({lhs}) {cmp_op} ({rhs}) ? ({lhs}) : ({rhs}))" -def _c_expr_binary(py_op, lhs, rhs, field_dtype=None): - """Return a C expression for a binary op: e.g., ``(x->a) + (y->a)``. +def _c_truediv_expr(lhs, rhs, field_dtype): + """Return a C expression for Python-semantics true division. + + C ``/`` on two integers is integer division; Python's ``/`` always + divides in floating point and only then does the result land back in the + integer field. The difference is visible whenever the exact quotient + doesn't fit a double (``10**18 / 3`` is 333333333333333333 in C but + 333333333333333312 through float64) and it is what makes integer + division by zero trap: ``(double) 7 / 0`` is ``inf``, but ``7 / 0`` in + integers raises SIGFPE and takes the whole process down. + + The quotient then has to come back down into the field, and converting a + double that doesn't fit the destination is undefined in C and poison in + LLVM, so the two paths need not agree. The two range escapes that a + same-signedness division can take are ruled out instead, matching the + guards :func:`_expr_binary` emits for the cfunc. A 64-bit field can still + escape through the ``(double)`` conversion (``(2**63 - 1) / 1`` rounds up + to ``2**63``). Both paths land on the same hardware conversion there + rather than on defined behaviour, so they agree with each other but need + not agree across machines; measured saturating to ``INT64_MAX`` on arm64. + + A complex field spells the zero-divisor case out rather than leaning on + C99 Annex G. SuiteSparse's JIT compiles with ``-fcx-limited-range`` under + GCC, which replaces the Annex G division with the naive formula, so + ``z / 0`` came out ``nan+nanj`` on Linux while clang (no such flag) and + the cfunc gave numpy's infinities. Dividing the parts by ``0.0`` as reals + sidesteps the flag entirely and matches the cfunc's spelling (see + ``_expr_binary``). ``CMPLX`` rather than arithmetic on ``I``: under the + naive formula ``inf * I`` multiplies out to ``nan``, which is the exact + failure being avoided. ``rhs == 0`` on a ``_Complex`` operand compares + both parts, mirroring the cfunc's ``y != 0``. + """ + kind = field_dtype.kind + if kind == "c": + if field_dtype.itemsize == 8: + cmplx, creal, cimag, zero = "CMPLXF", "crealf", "cimagf", "0.0f" + else: + cmplx, creal, cimag, zero = "CMPLX", "creal", "cimag", "0.0" + return ( + f"(({rhs}) == 0 " + f"? {cmplx}({creal}({lhs}) / {zero}, {cimag}({lhs}) / {zero}) " + f": ({lhs}) / ({rhs}))" + ) + if kind == "f": + return f"({lhs}) / ({rhs})" + quotient = f"(double)({lhs}) / (double)({rhs})" + if kind == "b": + # C converts a double to _Bool by comparing against zero, so there is + # no range to leave and nothing to guard. + return f"({quotient})" + if kind == "u": + return f"(({rhs}) == 0 ? 0 : ({quotient}))" + type_min = _C_INT_MIN[field_dtype.itemsize] + return ( + f"(({rhs}) == 0 ? 0 : " + f"((({lhs}) == {type_min} && ({rhs}) == -1) ? {type_min} : ({quotient})))" + ) + - ``field_dtype`` is the numpy dtype of the *result* element. It is - consulted for ``floordiv`` (``//``), which needs Python ``//`` semantics - rather than C ``/``, and for ``min`` / ``max``, which need C99 ``fmin`` / - ``fmax`` on floating-point fields so NaN is ignored rather than ordered. - The remaining ops are type-agnostic at the C level. +def _c_int_floordiv_expr(lhs, rhs, field_dtype): + """Return a C expression for integer floor division. + + Three ways C ``/`` differs from what the cfunc computes: + + - It truncates toward zero; Python ``//`` floors. Subtract one when the + operands have different signs and the remainder is non-zero. + - It raises SIGFPE (process death, not an exception) when the divisor is + zero. Numba's numpy error model returns 0, as does ``np.floor_divide``. + - It raises SIGFPE on ``INT_MIN / -1``, whose true quotient is not + representable. numpy wraps to ``INT_MIN``, which is what the cfunc + reaches through negation (see ``_expr_binary``). """ - if py_op in _FUNC_BINARY_OPS: - return _c_minmax_expr(py_op, lhs, rhs, field_dtype) - if py_op == "//": - return _c_floordiv_expr(lhs, rhs, field_dtype) - c_op = _C_INFIX_OPS.get(py_op, py_op) - return f"({lhs}) {c_op} ({rhs})" + if field_dtype.kind in ("u", "b"): + # Both operands are non-negative, so truncation already floors. + return f"(({rhs}) == 0 ? 0 : (({lhs}) / ({rhs})))" + type_min = _C_INT_MIN[field_dtype.itemsize] + floored = f"({lhs}) / ({rhs}) - ((({lhs}) % ({rhs}) != 0) && ((({lhs}) < 0) != (({rhs}) < 0)))" + return ( + f"(({rhs}) == 0 ? 0 : " + f"((({lhs}) == {type_min} && ({rhs}) == -1) ? {type_min} : ({floored})))" + ) -def _c_floordiv_expr(lhs, rhs, field_dtype): - """Return a C expression for Python-semantics floor division. +def _c_float_floordiv_stmt(target, lhs, rhs, field_dtype): + """Return a C block computing floating-point ``//`` into ``target``. - Python ``//`` is floor (rounds toward negative infinity); C ``/`` is - trunc toward zero for ints and true division for floats. The two only - agree for non-negative integer operands; for everything else the JIT - path silently disagreed with the Numba cfunc path before this helper. + ``floor(a / b)`` is not floor division. numpy and CPython both compute + the quotient from the remainder (``(a - fmod(a, b)) / b``) and snap it, + which is more accurate and handles infinities differently: ``1.0 // 0.1`` + is 9.0 but ``floor(1.0 / 0.1)`` is 10.0, and ``inf // 2.0`` is NaN but + ``floor(inf / 2.0)`` is ``inf``. This mirrors numpy's ``npy_divmod`` + (which CPython's ``float_divmod`` matches), so the JIT kernel, the Numba + cfunc, and ``np.floor_divide`` all agree. - Float fields use ``floor()`` / ``floorf()`` from ````, which is - available in the JIT kernel via SuiteSparse's include chain - (``GraphBLAS.h`` -> ````). Signed integer fields use the - standard trunc-to-floor adjustment. Unsigned integers don't need - adjusting because both operands are non-negative. + The sign fix-up needs the remainder twice and the quotient twice, so + this emits a statement with temporaries rather than one expression. """ - if field_dtype is None: - # Caller didn't pass dtype info. The C ``/`` semantics match Python - # ``//`` for non-negative integer operands only. - return f"({lhs}) / ({rhs})" - kind = field_dtype.kind - if kind == "f": - if field_dtype.itemsize == 4: - return f"floorf((float)({lhs}) / (float)({rhs}))" - return f"floor((double)({lhs}) / (double)({rhs}))" - if kind in ("u", "b"): - return f"({lhs}) / ({rhs})" - # Signed integer: trunc-toward-zero is one greater than floor when the - # signs of ``a`` and ``b`` differ and the division has a non-zero - # remainder; subtract 1 in that case. - return f"(({lhs}) / ({rhs}) - ((({lhs}) % ({rhs}) != 0) && ((({lhs}) < 0) != (({rhs}) < 0))))" + is_f32 = field_dtype.itemsize == 4 + ctype = "float" if is_f32 else "double" + sfx = "f" if is_f32 else "" + half = "0.5f" if is_f32 else "0.5" + one = "1.0f" if is_f32 else "1.0" + zero = "0.0f" if is_f32 else "0.0" + return ( + f"{{ {ctype} gb_a = ({lhs}) ; {ctype} gb_b = ({rhs}) ; {ctype} gb_q ; " + # A zero divisor is the one case numpy answers straight from the + # division: +-inf, or NaN for 0/0. + f"if (gb_b == {zero}) {{ gb_q = gb_a / gb_b ; }} " + f"else {{ {ctype} gb_m = fmod{sfx} (gb_a, gb_b) ; " + f"{ctype} gb_d = (gb_a - gb_m) / gb_b ; " + # NaN compares false both ways, so a NaN remainder skips the + # adjustment and carries through to the quotient. + f"if (gb_m != {zero} && ((gb_b < {zero}) != (gb_m < {zero}))) {{ gb_d -= {one} ; }} " + f"if (gb_d != {zero}) {{ gb_q = floor{sfx} (gb_d) ; " + f"if (gb_d - gb_q > {half}) {{ gb_q += {one} ; }} }} " + f"else {{ gb_q = copysign{sfx} ({zero}, gb_a / gb_b) ; }} }} " + f"{target} = gb_q ; }}" + ) def _c_expr_unary(py_op, operand, field_dtype=None): @@ -1001,7 +1188,7 @@ def _make_jit_c_definition(op_name, py_op, dtype, arity): # Pass the leaf dtype to the binary expression builder so # type-sensitive ops (currently floordiv) can emit correct C. assigns = " ".join( - f"z->{c} = {_c_expr_binary(py_op, f'x->{c}', f'y->{c}', leaf_dtype)} ;" + _c_assign_binary(py_op, f"z->{c}", f"x->{c}", f"y->{c}", leaf_dtype) for _py, c, leaf_dtype in leaves ) else: @@ -1014,7 +1201,7 @@ def _make_jit_c_definition(op_name, py_op, dtype, arity): size = reduce(mul, shape) if arity == 2: assigns = " ".join( - f"z->v[{i}] = {_c_expr_binary(py_op, f'x->v[{i}]', f'y->v[{i}]', base_dtype)} ;" + _c_assign_binary(py_op, f"z->v[{i}]", f"x->v[{i}]", f"y->v[{i}]", base_dtype) for i in range(size) ) else: diff --git a/graphblas/tests/test_op.py b/graphblas/tests/test_op.py index 7987f0f54..c4c14b384 100644 --- a/graphblas/tests/test_op.py +++ b/graphblas/tests/test_op.py @@ -1625,6 +1625,234 @@ def test_udt_min_max_answer_what_the_builtin_dtype_answers(udt_op_path, np_dtype assert u.reduce(monoid.max[udt]).new().value[0] == 3.0, f"{udt_op_path} {data}" +@pytest.mark.skipif("not supports_udfs") +def test_udt_truediv_divides_in_floating_point(udt_op_path): + """``binary.truediv`` on integer fields must divide the way Python does. + + Regression: the JIT kernel emitted C ``/``, which is integer division. + ``10**18 / 3`` came out as 333333333333333333 under the JIT and + 333333333333333312 (float64, like numpy) through the cfunc, so the same + program gave different answers depending on whether a C compiler was + installed. + """ + udt = dtypes.register_anonymous( + np.dtype([("tdv_i", np.int64), ("tdv_j", np.int64)], align=True), "_TrueDivIntUDT" + ) + xs = [10**18, 10**18 + 1, 7, 22] + ys = [3, 3, 2, 7] + expected = (np.array(xs, np.int64) / np.array(ys, np.int64)).astype(np.int64) + assert expected[0] == 333333333333333312 # not 333333333333333333 + v, w = _udt_vectors(udt, xs, ys) + result = binary.truediv(v & w).new() + got = [result[i].new().value[0] for i in range(len(xs))] + assert got == list(expected), udt_op_path + + +@pytest.mark.skipif("not supports_udfs") +def test_udt_floordiv_matches_numpy_on_floats(udt_op_path): + """``binary.floordiv`` on float fields is not ``floor(a / b)``. + + Regression: the JIT kernel computed ``floor(a / b)``, which rounds + differently from the remainder-based algorithm numpy and CPython use and + treats infinities as ordinary values. ``1.0 // 0.1`` came out as 10.0 + instead of 9.0, and ``inf // 2.0`` as ``inf`` instead of NaN, while the + cfunc agreed with numpy all along. + """ + nan = float("nan") + inf = float("inf") + # ``floor(a / b)`` disagrees with numpy on the first four pairs: inf and + # -inf where numpy gives NaN, 10.0 rather than 9.0 for 1.0 // 0.1, and + # -0.0 rather than -1.0 for -2.0 // inf. + xs = [inf, -inf, 1.0, -2.0, 2.0, -2.0, 0.0, nan, -7.0, 7.0, -0.0, 7.5] + ys = [2.0, 2.0, 0.1, inf, 0.0, 0.0, 0.0, 2.0, 2.0, -2.0, 4.0, 2.5] + for np_dtype, name in ((np.float64, "_FloorDivF64UDT"), (np.float32, "_FloorDivF32UDT")): + udt = dtypes.register_anonymous( + np.dtype([("fdv_a", np_dtype), ("fdv_b", np_dtype)], align=True), name + ) + with np.errstate(divide="ignore", invalid="ignore"): + expected = np.floor_divide(np.array(xs, np_dtype), np.array(ys, np_dtype)) + v, w = _udt_vectors(udt, xs, ys) + result = binary.floordiv(v & w).new() + got = np.array([result[i].new().value[0] for i in range(len(xs))], np_dtype) + np.testing.assert_array_equal(got, expected, err_msg=f"{udt_op_path} {np_dtype.__name__}") + # ``assert_array_equal`` reads -0.0 and 0.0 as equal, so the sign of a + # zero quotient needs its own assertion. It is the whole job of the + # ``copysign`` branch the JIT kernel emits for an exact-zero result. + np.testing.assert_array_equal( + np.signbit(got), + np.signbit(expected), + err_msg=f"{udt_op_path} {np_dtype.__name__} sign of zero", + ) + + +@pytest.mark.skipif("not supports_udfs") +def test_udt_integer_division_by_zero_is_defined(udt_op_path): + """Integer division by zero must return a value rather than trap. + + The JIT kernel divided in integers, where a zero divisor is undefined + behaviour: on x86-64 ``idiv`` raises #DE, which is SIGFPE and process + death rather than an exception. AArch64's ``sdiv`` returns 0 and does not + trap, so this cannot be exhibited on an arm64 machine. The same trap + fires on ``INT_MIN / -1``, whose quotient is not representable. + + The values below are a choice, not a discovery. A quotient that doesn't + fit the field is an undefined conversion in both C and LLVM, and the two + answered differently per field width, so both paths now rule those cases + out: a zero divisor gives 0, as ``np.floor_divide`` does, and + ``INT_MIN // -1`` wraps to ``INT_MIN``, as numpy does. + """ + signed = dtypes.register_anonymous( + np.dtype([("dvz_a", np.int32), ("dvz_b", np.int8)], align=True), "_DivZeroSignedUDT" + ) + unsigned = dtypes.register_anonymous( + np.dtype([("dvz_c", np.uint32), ("dvz_d", np.uint64)], align=True), "_DivZeroUnsignedUDT" + ) + v, w = _udt_vectors(signed, [7, -7, 100, -128], [0, 0, 3, -1]) + for gb_op in (binary.truediv, binary.floordiv): + result = gb_op(v & w).new() + got = [result[i].new().value[0] for i in range(4)] + assert got[:2] == [0, 0], f"{udt_op_path} {gb_op.name}" + assert got[2] == 33, f"{udt_op_path} {gb_op.name}" # 100 / 3 truncates either way + # ``-128 // -1`` is the second trapping case; numpy wraps it to INT8_MIN. + result = binary.floordiv(v & w).new() + assert result[3].new().value[1] == np.iinfo(np.int8).min, udt_op_path + + v, w = _udt_vectors(unsigned, [7, 9, 100, 5], [0, 0, 3, 2]) + for gb_op in (binary.truediv, binary.floordiv): + result = gb_op(v & w).new() + got = [result[i].new().value[0] for i in range(4)] + assert got == [0, 0, 33, 2], f"{udt_op_path} {gb_op.name}" + + # Floor division still floors for signed operands of mixed sign. + v, w = _udt_vectors(signed, [-7, 7, -9, 11], [2, -2, 2, 3]) + result = binary.floordiv(v & w).new() + got = [result[i].new().value[0] for i in range(4)] + assert got == [-4, -4, -5, 3], udt_op_path + + +@pytest.mark.skipif("not supports_udfs") +@pytest.mark.skipif("not dtypes._supports_complex") +def test_udt_complex_truediv_by_zero(udt_op_path): + """``binary.truediv`` on a complex field survives a zero divisor. + + Numba's complex division raises ``ZeroDivisionError`` unconditionally, + outside the error model's control, so the cfunc left the element unwritten + while the JIT kernel returned numpy's infinities. Reading back the + abandoned element gave uninitialized memory, or the previous element's + answer, either of which looks like a plausible value. + """ + udt = dtypes.register_anonymous( + np.dtype([("cxz_a", np.complex128)], align=True), "_ComplexDivZeroUDT" + ) + xs = [3 + 4j, 0j, 1 + 1j, 2 - 2j] + ys = [0j, 0j, 2 + 0j, 0j] + v, w = _udt_vectors(udt, xs, ys) + result = binary.truediv(v & w).new() + got = np.array([result[i].new().value[0] for i in range(len(xs))]) + with np.errstate(divide="ignore", invalid="ignore"): + expected = np.array(xs) / np.array(ys) + np.testing.assert_array_equal(got, expected, err_msg=udt_op_path) + + +@pytest.mark.skipif("not supports_udfs") +# 136-byte UDT, which SS < 9 rejects; see test_udt_large_array. +@pytest.mark.skipif( + "ss_version_major < 9", + reason="SuiteSparse < 9 rejects a 136-byte UDT on builds without VLA support", +) +def test_udt_float_truediv_by_zero_is_infinite(udt_op_path): + """A zero divisor on a float field gives numpy's infinity, not a lost element. + + Unlike the integer case, nothing here is guarded: the generated code + divides and lets IEEE produce the infinity. That only holds because the + generated wrapper is compiled under Numba's numpy error model, which + nothing else in the suite pins down. + """ + udt = dtypes.register_anonymous(np.dtype((np.float64, (17,))), "_FloatDivZeroArr17") + xs = [1.0, -1.0, 0.0, 6.0] + ys = [0.0, 0.0, 0.0, 3.0] + v, w = _udt_vectors(udt, xs, ys) + result = binary.truediv(v & w).new() + got = np.array([result[i].new().value[0] for i in range(len(xs))]) + with np.errstate(divide="ignore", invalid="ignore"): + expected = np.array(xs) / np.array(ys) + np.testing.assert_array_equal(got, expected, err_msg=udt_op_path) + + +@pytest.mark.skipif("not supports_udfs") +@pytest.mark.slow +def test_udt_array_ops_match_record_ops(udt_op_path): + """The array-UDT codegen carries the same division and NaN fixes as records. + + Records and flat arrays go through separate branches in both the Numba + and the JIT C generators, so each fix has to land in both. + """ + nan = float("nan") + inf = float("inf") + float_udt = dtypes.register_anonymous(np.dtype((np.float64, (13,))), "_ArrOpsF64") + xs = [inf, 1.0, -2.0, nan, -7.0, 2.0] + ys = [2.0, 0.1, inf, 2.0, 2.0, 1.0] + v, w = _udt_vectors(float_udt, xs, ys) + # The reference computations touch inf and nan, and numpy raises the FP + # invalid flag for them on some platforms (Linux and Windows, via fmod) + # but not others; pyproject promotes the RuntimeWarning to an error. + with np.errstate(divide="ignore", invalid="ignore"): + expected_floordiv = np.floor_divide(np.array(xs), np.array(ys)) + expected_min = np.fmin(np.array(xs), np.array(ys)) + np.testing.assert_array_equal( + [binary.floordiv(v & w).new()[i].new().value[0] for i in range(len(xs))], + expected_floordiv, + err_msg=udt_op_path, + ) + # ``np.fmin``, not ``np.minimum``: ``binary.min`` is SuiteSparse's + # ``GrB_MIN_FP64``, which ignores a NaN operand rather than propagating it. + np.testing.assert_array_equal( + [binary.min(v & w).new()[i].new().value[0] for i in range(len(xs))], + expected_min, + err_msg=udt_op_path, + ) + + int_udt = dtypes.register_anonymous(np.dtype((np.int64, (6,))), "_ArrOpsI64") + ixs = [10**18, 7, -7, 100, -9, 5] + iys = [3, 0, 0, 3, 2, 2] + v, w = _udt_vectors(int_udt, ixs, iys) + result = binary.truediv(v & w).new() + got = [result[i].new().value[0] for i in range(len(ixs))] + assert got == [333333333333333312, 0, 0, 33, -4, 2], udt_op_path + result = binary.floordiv(v & w).new() + got = [result[i].new().value[0] for i in range(len(ixs))] + assert got == [333333333333333333, 0, 0, 33, -5, 2], udt_op_path + + +@pytest.mark.skipif("not supports_udfs") +@pytest.mark.slow +def test_udt_multidim_array_ops_match_numpy(udt_op_path): + """Built-in ops on a multi-dimensional array UDT agree with numpy on both paths. + + The JIT C typedef flattens any rank to ``double v [N]`` and the Numba + wrapper walks the same flat run, so a 2-D UDT covers codegen that the 1-D + cases reach only by accident of both being contiguous. + """ + udt = dtypes.register_anonymous(np.dtype((np.float64, (3, 2))), "_ArrOps2D") + xs = [1.0, -7.0, float("inf"), 2.0] + ys = [0.1, 2.0, 2.0, 0.0] + v, w = _udt_vectors(udt, xs, ys) + for gb_op, reference in ( + (binary.floordiv, np.floor_divide), + (binary.truediv, np.true_divide), + # ``fmin`` rather than ``minimum``: these inputs carry no NaN, so the + # two agree here, but ``binary.min`` is the NaN-ignoring one. + (binary.min, np.fmin), + ): + result = gb_op(v & w).new() + element = result[0].new().value + assert element.shape == (3, 2) + got = np.array([result[i].new().value[1, 1] for i in range(len(xs))]) + with np.errstate(divide="ignore", invalid="ignore"): + expected = reference(np.array(xs), np.array(ys)) + np.testing.assert_array_equal(got, expected, err_msg=f"{udt_op_path} {gb_op.name}") + + @pytest.mark.skipif("not supports_udfs") @pytest.mark.slow def test_udt_tuple_return_binaryop(record_udt):