About the last two points, the problem is that we are dividing by an integer scalar? So there is a cast necessary and the promotion isn't quite trivial unfortunately (also because it is a scalar, so the legacy code kicks in, which inspects the value :/, trying to push that out, but its slow progress).
So we cannot easily avoid the "Resolver" because of the scalar integer here (or I couldn't think of a way).

Now, I admit, I think some of those paths that need to fall back got slower when scalars were involved, maybe that is part of what you see now. (I remember there was one path that got slower, while most others got faster, but I forgot which it was)

Te overlap check, it is behind a if (op[nin] == NULL) { check so shouldn't get called if out is not passed explicitly? But, np.mean does pass out explicitly. In fact, out is identical to the first input there (which is OK, the overlap detection should say that they are the same).
Now, what might be possible is to do an identity check for these situations in PyArray_EQUIVALENTLY_ITERABLE_OVERLAP_OK (or similar), since this is a common situation for in-place operators!

A fun way to improve the __array_function__ __array_ufunc__ might be to do rcount = int(rcount) to make it a python integer.

It is good to look at this closer, but mean is a pure Python function (you probably already found it in numpy/core/_methods.py, and to really knock away the last bit of overheads, we might have to reconsider that?!

Interesting, when where= was introduced, it may be that we switched from getting 0 to getting np.intp(0) out from the counting function. That might well be the big slow-down here, since as I had said, returning int might hit a few faster paths (i.e. avoid the _array_ufunc__ stuff).

Just to confirm that:

diff --git a/numpy/core/_methods.py b/numpy/core/_methods.py
index a239e2c87..b6ec15714 100644
--- a/numpy/core/_methods.py
+++ b/numpy/core/_methods.py
@@ -71,7 +71,7 @@ def _count_reduce_items(arr, axis, keepdims=False, where=True):
             axis = tuple(range(arr.ndim))
         elif not isinstance(axis, tuple):
             axis = (axis,)
-        items = nt.intp(1)
+        items = 1
         for ax in axis:
             items *= arr.shape[mu.normalize_axis_index(ax, arr.ndim)]
     else:

speeds things up quite a lot in your example. I do not know if scalar math is involved here (or how much), so my scalar-math update may help. But I suspect this also helps for the overrides etc.

With the PRs linking to this issue merged the difference is smaller. np.mean is a python function, so some overhead is to be expected.

There are two more items that can be addressed to reduce the difference:

• PERF: Drastically lower the overhead of __array_function__ #20790 by @seberg
• Avoiding the legacy resolving. See PERF: Use python integer on _count_reduce_items #21465 (comment)

Closing this issue

Now that branching happened, I can put you on the trace that our generic objects should probably not have a custom tp_alloc at all ;).
There is a tiny chance that this affects downstream subclasses, but I doubt it, since user dtypes shouldn't really be able to require NEEDS_INIT for anything reasonable.

Now that branching happened, I can put you on the trace that our generic objects should probably not have a custom tp_alloc at all ;). There is a tiny chance that this affects downstream subclasses, but I doubt it, since user dtypes shouldn't really be able to require NEEDS_INIT for anything reasonable.

@seberg Not sure I understand this.

The np.generic scalar provides a tp_alloc which is inherited by all scalars. This does some unnecessary initialization (only void really needs to initialize).
Removing that tp_alloc() will mean that np.float64, etc. use the default Python allocation directly, which probably removes a tiny bit of overhead.

Do you mean we can modify PyGenericArrType_Type.tp_alloc = gentype_alloc here https://github.com/numpy/numpy/blob/main/numpy/core/src/multiarray/scalartypes.c.src#L3874 and set if to _PyObject_New?

The difference between gentype_alloc and _PyObject_New is allocation of one item more, a memset(obj, 0, size) and a special case for type->tp_itemsize == 0.
Is there are information why those differences are there?

@seberg

Hmmm, the generic version is PyType_GenericAlloc I think (i.e. what is used if you just delete that line of code).
Seems that does exactly the same stuff, except that it checks for whether it needs storage for the GC.

I thought I had timed it to make things faster if I remove the custom version. Maybe a fluke, or maybe reusing that version is actually good (because it is used all over the place).

EDIT: I suppose in principle, we could not remove it, but just cut it down a lot (no need to memset or itemsize check for our scalars really)