BUG: Ensure that scalar binops prioritize __array_ufunc__

If array-ufunc is implemented, we must call always use it for all
operators (that seems to be the promise).

If __array_function__ is defined we are in the clear w.r.t. recursion
because the object is either an array (can be unpacked, but already
checked earlier now also), or it cannot call the ufunc without unpacking
itself (otherwise it would cause recursion).

There is an oddity about `__array_wrap__`.  Rather than trying to
do odd things to deal with it, I added a comment explaining why it
doens't matter (roughly: don't use our scalar priority if you want
to be sure to get a chance).

Author: seberg

numpy/_core/src/multiarray/scalartypes.c.src

@@ -194,15 +194,30 @@ find_binary_operation_path(
             PyLong_CheckExact(other) ||
             PyFloat_CheckExact(other) ||
             PyComplex_CheckExact(other) ||
-            PyBool_Check(other)) {
+            PyBool_Check(other) ||
+            PyArray_Check(other)) {
         /*
          * The other operand is ready for the operation already.  Must pass on
          * on float/long/complex mainly for weak promotion (NEP 50).
          */
-        Py_INCREF(other);
-        *other_op = other;
+        *other_op = Py_NewRef(other);
         return 0;
     }
+    /*
+     * If other has __array_ufunc__ always use ufunc.  If array-ufunc was None
+     * we already deferred.  And any custom object with array-ufunc cannot call
+     * our ufuncs without preventing recursion.
+     * It may be nice to avoid double lookup in `BINOP_GIVE_UP_IF_NEEDED`.
+     */
+    PyObject *attr = PyArray_LookupSpecial(other, npy_interned_str.array_ufunc);
+    if (attr != NULL) {
+        Py_DECREF(attr);
+        *other_op = Py_NewRef(other);
+        return 0;
+    }
+    else if (PyErr_Occurred()) {
+        PyErr_Clear(); /* TODO[gh-14801]: propagate crashes during attribute access? */
+    }
 
     /*
      * Now check `other`.  We want to know whether it is an object scalar
@@ -216,7 +231,13 @@ find_binary_operation_path(
     }
 
     if (!was_scalar || PyArray_DESCR(arr)->type_num != NPY_OBJECT) {
-        /* The array is OK for usage and we can simply forward it
+        /* 
+         * The array is OK for usage and we can simply forward it.  There
+         * is a theoretical subtlety here:  If the other object implements
+         * `__array_wrap__`, we may ignore that.  However, this only matters
+         * if the other object has the identical `__array_priority__` and
+         * additionally already deferred back to us.
+         * (`obj + scalar` and `scalar + obj` are not symmetric.)
          *
          * NOTE: Future NumPy may need to distinguish scalars here, one option
          *       could be marking the array.

numpy/_core/tests/test_multiarray.py

@@ -4025,6 +4025,18 @@ class LowPriority(np.ndarray):
         assert res.shape == (3,)
         assert res[0] == 'result'
 
+    @pytest.mark.parametrize("scalar", [
+            np.longdouble(1), np.timedelta64(120, 'm')])
+    @pytest.mark.parametrize("op", [operator.add, operator.xor])
+    def test_scalar_binop_guarantees_ufunc(self, scalar, op):
+        # Test that __array_ufunc__ will always cause ufunc use even when
+        # we have to protect some other calls from recursing (see gh-26904).
+        class SomeClass:
+            def __array_ufunc__(self, ufunc, method, *inputs, **kw):
+                return "result"
+
+        assert SomeClass() + np.longdouble(1) == "result"
+        assert np.longdouble(1) + SomeClass() == "result"
 
     def test_ufunc_override_normalize_signature(self):
         # gh-5674

numpy/_core/tests/test_scalarmath.py

@@ -27,7 +27,7 @@
 floating_types = np.floating.__subclasses__()
 complex_floating_types = np.complexfloating.__subclasses__()
 
-objecty_things = [object(), None]
+objecty_things = [object(), None, np.array(None, dtype=object)]
 
 binary_operators_for_scalars = [
     operator.lt, operator.le, operator.eq, operator.ne, operator.ge,