ENH: support reduction loops in ufunc.reduceat - #32212

Broader suggestion, probably for a shared prep patch rather than this PR: the
[acc_0.., x, out_0..] index arithmetic (nout, nout + 1 + i,
2*nout + 1) is now hand-written at every site — here, the seed loops below,
reduce_loop:

Can I do this later? I don't feel super confident I can extract the right macros before seeing how the reduceat and accumulate use them? There is potential to replace loops with macros but I don't easily see the right design right now which is why I would prefer it to do it after. I understand that we are repeating these loops but I don't think they will diverge cause they are trivial in what they do.

This is UB for NaN, infinity, and values outside the representable range of intp. This is an internal test helper but it's at minimum bad style to intentionally add UB in C code like this. In principle compilers could also note that UB is happening here and generate nasal demons.

We don't have a C helper for this right now so you could either add one or add NPY_METH_REQUIRES_PYAPI and use e.g. PyLong_FromDouble and then PyArray_PyIntAsIntP to do this conversion. My AI model identified a few other places in the numpy implementation that have fishy casts like this. IMO that probably argues for landing a C helper and fixing these cases with the helper in a PR and then rebase this PR on top of the PR with the fix.

There is an _arange_safe_ceil_to_intp helper but it's buggy so some of these are slipping through that:

1. arange’s length conversion — ctors.c:3089. Its check permits ivalue <= (double)NPY_MAX_INTP. With 64-bit intp, that upper bound rounds to 2⁶³, which is already outside the integer range. The subsequent cast therefore has UB at that boundary. This corrects my earlier characterization: _arange_safe_ceil_to_intp has checks, but they are incomplete.

2. casting="same_value" validation — lowlevel_strided_loops.c.src:989. The integer maximum rounds upward during comparison here too, allowing an invalid cast inside the round-trip check. This produces an observable bug on this machine:

a = np.array([float(2**63)])
a.astype(np.int64, casting="same_value")
# array([9223372036854775807]), with RuntimeWarning

The input is exactly 9223372036854775808. It should reject the conversion. Converting the incorrect result back to double rounds upward again, concealing the change from the round-trip check. The analogous uint64 boundary also passes incorrectly.

3. Timedelta multiplication and division by doubles — loops.c.src:1136, also lines 1157 and 1178. These check npy_isfinite(result), but finite values can still exceed the int64 range. Examples include:

np.timedelta64(1, "ns") * 1e20
np.timedelta64(1, "ns") / 1e-20

Both reach an out-of-range cast. Here they produce the maximum timedelta with a warning; that result is not guaranteed by C.

I confirmed these casts with Clang’s float-cast-overflow sanitizer using isolated functions extracted from the source, alongside the Python reproductions above.

So there is a concrete reuse case for a small, allocation-free C helper that returns conversion status, letting each caller choose its error handling. An intp variant would cover arange and the PR fixture; timedelta needs an int64 variant even on 32-bit platforms. The casting machinery would need additional destination types. Centralizing the boundary calculation would avoid repeating the deceptively unsafe “compare against integer maximum converted to double” pattern.

Would something like this suffice here or do you want to go for the separate PR with a more general helper?

diff --git a/numpy/_core/src/umath/_reduction_loop_tests.c b/numpy/_core/src/umath/_reduction_loop_tests.c
index 33228b7430..8c7bbce279 100644
--- a/numpy/_core/src/umath/_reduction_loop_tests.c
+++ b/numpy/_core/src/umath/_reduction_loop_tests.c
@@ -226,6 +226,21 @@ object_minimummaximum_get_reduction_loop(
  * The second output is the maximum converted to its dtype, so reduce and
  * reduceat seed it by casting the first float64 element of the reduction.
  */
+/* Convert to intp, or set a ValueError if `val` is NaN or out of range. */
+static int
+double_to_intp(double val, npy_intp *out)
+{
+    /* Both bounds are powers of two, so they are exact as doubles. */
+    if (!(val >= (double)NPY_MIN_INTP && val < -(double)NPY_MIN_INTP)) {
+        PyErr_SetString(PyExc_ValueError,
+                        "maximum is NaN or out of the intp range");
+        return -1;
+    }
+    *out = (npy_intp)val;
+    return 0;
+}
+
+
 static int
 double_minimum_intp_maximum_loop(PyArrayMethod_Context *NPY_UNUSED(context),
         char *const data[], npy_intp const dimensions[],
@@ -238,8 +253,10 @@ double_minimum_intp_maximum_loop(PyArrayMethod_Context *NPY_UNUSED(context),
     for (npy_intp i = 0; i < n; i++) {
         double a = *(double *)in1;
         double b = *(double *)in2;
+        if (double_to_intp(nan_max(a, b), (npy_intp *)out2) < 0) {
+            return -1;
+        }
         *(double *)out1 = nan_min(a, b);
-        *(npy_intp *)out2 = (npy_intp)nan_max(a, b);
         in1 += strides[0]; in2 += strides[1];
         out1 += strides[2]; out2 += strides[3];
     }
@@ -259,7 +276,10 @@ double_minimum_intp_maximum_reduce_loop(

     for (npy_intp i = 0; i < n; i++) {
         double val = *(double *)x;
-        npy_intp ival = (npy_intp)val;
+        npy_intp ival;
+        if (double_to_intp(val, &ival) < 0) {
+            return -1;
+        }
         npy_intp cur_max = *(npy_intp *)acc_max;
         *(double *)out_min = nan_min(*(double *)acc_min, val);
         *(npy_intp *)out_max = cur_max >= ival ? cur_max : ival;
@@ -281,7 +301,7 @@ minimum_intp_maximum_get_reduction_loop(
 {
     *out_loop = &double_minimum_intp_maximum_reduce_loop;
     *out_transferdata = NULL;
-    *flags = NPY_METH_NO_FLOATINGPOINT_ERRORS;
+    *flags = NPY_METH_REQUIRES_PYAPI;
     return 0;
 }

@@ -560,7 +580,7 @@ _reduction_loop_tests_exec(PyObject *m)
     if (add_mixed_minimummaximum(m, "minimum_intp_maximum", NPY_INTP,
             &double_minimum_intp_maximum_loop,
             &minimum_intp_maximum_get_reduction_loop,
-            NPY_METH_NO_FLOATINGPOINT_ERRORS) < 0) {
+            NPY_METH_REQUIRES_PYAPI) < 0) {
         return -1;
     }
     if (add_mixed_minimummaximum(m, "minimum_object_maximum", NPY_OBJECT,
@@ -584,9 +604,10 @@ static struct PyModuleDef_Slot _reduction_loop_tests_slots[] = {
 #endif
 #if PY_VERSION_HEX >= 0x030d00f0  // Python 3.13+
     /*
-     * No GIL is needed: the inner loops never touch Python objects,
-     * apart from the object ones, which declare NPY_METH_REQUIRES_PYAPI
-     * so the iterator does not release the GIL around them.
+     * No GIL is needed: the inner loops never touch the Python API,
+     * apart from the object and intp ones, which declare
+     * NPY_METH_REQUIRES_PYAPI so the iterator does not release the GIL
+     * around them.
      */
     {Py_mod_gil, Py_MOD_GIL_NOT_USED},
 #endif

Let's fix this in a followup.