ENH: Annotate the arithmetic operations of `ndarray` and `generic` (#17273)

* ENH: Added annotations for arithmetic-based magic methods

* TST: Added arithmetic tests

* TST: Moved a number of tests to `arithmetic.py`

* ENH: Ensure that objects annotated as `number` support arithmetic operations

* MAINT: Arithmetic operations on 0d arrays return scalars

* MAINT: Clarify the type of generics returned by `ufunc.__call__`

* TST: Added more arithmetic tests

* MAINT: Use `_CharLike` when both `str` and `bytes` are accepted

* MAINT: Change the `timedelta64` baseclass to `generic`

* MAINT: Add aliases for common scalar unions

* MAINT: Update the defition of `_NumberLike`

* MAINT: Replace `_NumberLike` with `_ComplexLike` in the `complexfloating` annotations

* MAINT: Move the callback protocols to their own module

* MAINT: Make `typing._callback` available at runtime

* DOC: Provide further clarification about callback protocols

* MAINT: Replace `_callback` with `_callable`

Addresses https://github.com/numpy/numpy/pull/17273#discussion_r485821346

The use of `__call__`-defining protocols is not limited to callbacks. The module name name & docstring now reflects this.

* MAINT: Removed `__add__` from `str_` and `bytes_`

Most `np.bytes_` / `np.str_` methods return their builtin `bytes` / `str` counterpart.
This includes addition.

* MAINT: Fix the return type of boolean division

Addresses https://github.com/numpy/numpy/pull/17273#discussion_r486271220

Dividing a `np.bool_` by an integer (or vice versa) always returns `float64`

* MAINT: Renamed all `_<x>Arithmetic` protocols to `_<x>Op

Addresses https://github.com/numpy/numpy/pull/17273#discussion_r486272745

* TST: Add tests for boolean division

* ENH: Make `np.number` generic w.r.t. its precision

* ENH,WIP: Add a mypy plugin for casting `np.number` instances to appropiate subclasses

* Revert "ENH,WIP: Add a mypy plugin for casting `np.number` instances to appropiate subclasses"

This reverts commit c526fb619d20902bfd77709c8983c7a7d5477c95.

* Revert "ENH: Make `np.number` generic w.r.t. its precision"

This reverts commit dbf20183cf7ff71e379cd1a165d07e1a1d643135.

* MAINT: Narow the definition of `_ComplexLike`

Addresses https://github.com/numpy/numpy/pull/17273#discussion_r490440238

* MAINT: Refined the return type of `unint + int` ops

`unsignedinteger + signedinteger` generally returns a `signedinteger` subclass.
The exception to this is `uint64 + signedinteger`, which returns `float64`.

Addresses https://github.com/numpy/numpy/pull/17273#discussion_r490442023

* MAINT: Use `_IntLike` and `_FloatLike` in the definition of `_ComplexLike`

9 files changed, 811 insertions, 89 deletions

diff --git a/numpy/__init__.pyi b/numpy/__init__.pyi
index 30d15ed12..e712801eb 100644
--- a/numpy/__init__.pyi
+++ b/numpy/__init__.pyi
@@ -5,6 +5,18 @@ from abc import abstractmethod
 
 from numpy.core._internal import _ctypes
 from numpy.typing import ArrayLike, DtypeLike, _Shape, _ShapeLike
+from numpy.typing._callable import (
+    _BoolOp,
+    _BoolSub,
+    _BoolTrueDiv,
+    _TD64Div,
+    _IntTrueDiv,
+    _UnsignedIntOp,
+    _SignedIntOp,
+    _FloatOp,
+    _ComplexOp,
+    _NumberOp,
+)
 
 from typing import (
     Any,
@@ -646,23 +658,10 @@ class _ArrayOrScalarCommon(
     def __ne__(self, other): ...
     def __gt__(self, other): ...
     def __ge__(self, other): ...
-    def __add__(self, other): ...
-    def __radd__(self, other): ...
-    def __sub__(self, other): ...
-    def __rsub__(self, other): ...
-    def __mul__(self, other): ...
-    def __rmul__(self, other): ...
-    def __truediv__(self, other): ...
-    def __rtruediv__(self, other): ...
-    def __floordiv__(self, other): ...
-    def __rfloordiv__(self, other): ...
     def __mod__(self, other): ...
     def __rmod__(self, other): ...
     def __divmod__(self, other): ...
     def __rdivmod__(self, other): ...
-    # NumPy's __pow__ doesn't handle a third argument
-    def __pow__(self, other): ...
-    def __rpow__(self, other): ...
     def __lshift__(self, other): ...
     def __rlshift__(self, other): ...
     def __rshift__(self, other): ...
@@ -834,14 +833,26 @@ class ndarray(_ArrayOrScalarCommon, Iterable, Sized, Container):
     def __matmul__(self, other): ...
     def __imatmul__(self, other): ...
     def __rmatmul__(self, other): ...
+    def __add__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __radd__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __sub__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __rsub__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __mul__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __rmul__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __floordiv__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __rfloordiv__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __pow__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __rpow__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __truediv__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
+    def __rtruediv__(self, other: ArrayLike) -> Union[ndarray, generic]: ...
     # `np.generic` does not support inplace operations
-    def __iadd__(self, other): ...
-    def __isub__(self, other): ...
-    def __imul__(self, other): ...
-    def __itruediv__(self, other): ...
-    def __ifloordiv__(self, other): ...
+    def __iadd__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
+    def __isub__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
+    def __imul__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
+    def __itruediv__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
+    def __ifloordiv__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
+    def __ipow__(self: _ArraySelf, other: ArrayLike) -> _ArraySelf: ...
     def __imod__(self, other): ...
-    def __ipow__(self, other): ...
     def __ilshift__(self, other): ...
     def __irshift__(self, other): ...
     def __iand__(self, other): ...
@@ -857,6 +868,11 @@ class ndarray(_ArrayOrScalarCommon, Iterable, Sized, Container):
 # See https://github.com/numpy/numpy-stubs/pull/80 for more details.
 
 _CharLike = Union[str, bytes]
+_BoolLike = Union[bool, bool_]
+_IntLike = Union[int, integer]
+_FloatLike = Union[_IntLike, float, floating]
+_ComplexLike = Union[_FloatLike, complex, complexfloating]
+_NumberLike = Union[int, float, complex, number, bool_]
 
 class generic(_ArrayOrScalarCommon):
     @abstractmethod
@@ -869,6 +885,19 @@ class number(generic):  # type: ignore
     def real(self: _ArraySelf) -> _ArraySelf: ...
     @property
     def imag(self: _ArraySelf) -> _ArraySelf: ...
+    # Ensure that objects annotated as `number` support arithmetic operations
+    __add__: _NumberOp
+    __radd__: _NumberOp
+    __sub__: _NumberOp
+    __rsub__: _NumberOp
+    __mul__: _NumberOp
+    __rmul__: _NumberOp
+    __floordiv__: _NumberOp
+    __rfloordiv__: _NumberOp
+    __pow__: _NumberOp
+    __rpow__: _NumberOp
+    __truediv__: _NumberOp
+    __rtruediv__: _NumberOp
 
 class bool_(generic):
     def __init__(self, __value: object = ...) -> None: ...
@@ -876,6 +905,18 @@ class bool_(generic):
     def real(self: _ArraySelf) -> _ArraySelf: ...
     @property
     def imag(self: _ArraySelf) -> _ArraySelf: ...
+    __add__: _BoolOp[bool_]
+    __radd__: _BoolOp[bool_]
+    __sub__: _BoolSub
+    __rsub__: _BoolSub
+    __mul__: _BoolOp[bool_]
+    __rmul__: _BoolOp[bool_]
+    __floordiv__: _BoolOp[int8]
+    __rfloordiv__: _BoolOp[int8]
+    __pow__: _BoolOp[int8]
+    __rpow__: _BoolOp[int8]
+    __truediv__: _BoolTrueDiv
+    __rtruediv__: _BoolTrueDiv
 
 class object_(generic):
     def __init__(self, __value: object = ...) -> None: ...
@@ -892,10 +933,18 @@ class datetime64(generic):
         __format: Union[_CharLike, Tuple[_CharLike, _IntLike]] = ...,
     ) -> None: ...
     @overload
-    def __init__(self, __value: int, __format: Union[_CharLike, Tuple[_CharLike, _IntLike]]) -> None: ...
-    def __add__(self, other: Union[timedelta64, int]) -> datetime64: ...
-    def __sub__(self, other: Union[timedelta64, datetime64, int]) -> timedelta64: ...
-    def __rsub__(self, other: Union[datetime64, int]) -> timedelta64: ...
+    def __init__(
+        self,
+        __value: int,
+        __format: Union[_CharLike, Tuple[_CharLike, _IntLike]]
+    ) -> None: ...
+    def __add__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> datetime64: ...
+    def __radd__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> datetime64: ...
+    @overload
+    def __sub__(self, other: datetime64) -> timedelta64: ...
+    @overload
+    def __sub__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> datetime64: ...
+    def __rsub__(self, other: datetime64) -> timedelta64: ...
 
 # Support for `__index__` was added in python 3.8 (bpo-20092)
 if sys.version_info >= (3, 8):
@@ -911,8 +960,20 @@ class integer(number):  # type: ignore
     # NOTE: `__index__` is technically defined in the bottom-most
     # sub-classes (`int64`, `uint32`, etc)
     def __index__(self) -> int: ...
-
-class signedinteger(integer): ...  # type: ignore
+    __truediv__: _IntTrueDiv
+    __rtruediv__: _IntTrueDiv
+
+class signedinteger(integer):  # type: ignore
+    __add__: _SignedIntOp
+    __radd__: _SignedIntOp
+    __sub__: _SignedIntOp
+    __rsub__: _SignedIntOp
+    __mul__: _SignedIntOp
+    __rmul__: _SignedIntOp
+    __floordiv__: _SignedIntOp
+    __rfloordiv__: _SignedIntOp
+    __pow__: _SignedIntOp
+    __rpow__: _SignedIntOp
 
 class int8(signedinteger):
     def __init__(self, __value: _IntValue = ...) -> None: ...
@@ -926,24 +987,36 @@ class int32(signedinteger):
 class int64(signedinteger):
     def __init__(self, __value: _IntValue = ...) -> None: ...
 
-class timedelta64(signedinteger):
+class timedelta64(generic):
     def __init__(
         self,
         __value: Union[None, int, _CharLike, dt.timedelta, timedelta64] = ...,
         __format: Union[_CharLike, Tuple[_CharLike, _IntLike]] = ...,
     ) -> None: ...
-    @overload
-    def __add__(self, other: Union[timedelta64, int]) -> timedelta64: ...
-    @overload
-    def __add__(self, other: datetime64) -> datetime64: ...
-    def __sub__(self, other: Union[timedelta64, int]) -> timedelta64: ...
-    @overload
-    def __truediv__(self, other: timedelta64) -> float: ...
-    @overload
-    def __truediv__(self, other: float) -> timedelta64: ...
+    def __add__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> timedelta64: ...
+    def __radd__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> timedelta64: ...
+    def __sub__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> timedelta64: ...
+    def __rsub__(self, other: Union[timedelta64, _IntLike, _BoolLike]) -> timedelta64: ...
+    def __mul__(self, other: Union[_FloatLike, _BoolLike]) -> timedelta64: ...
+    def __rmul__(self, other: Union[_FloatLike, _BoolLike]) -> timedelta64: ...
+    __truediv__: _TD64Div[float64]
+    __floordiv__: _TD64Div[signedinteger]
+    def __rtruediv__(self, other: timedelta64) -> float64: ...
+    def __rfloordiv__(self, other: timedelta64) -> signedinteger: ...
     def __mod__(self, other: timedelta64) -> timedelta64: ...
 
-class unsignedinteger(integer): ...  # type: ignore
+class unsignedinteger(integer):  # type: ignore
+    # NOTE: `uint64 + signedinteger -> float64`
+    __add__: _UnsignedIntOp
+    __radd__: _UnsignedIntOp
+    __sub__: _UnsignedIntOp
+    __rsub__: _UnsignedIntOp
+    __mul__: _UnsignedIntOp
+    __rmul__: _UnsignedIntOp
+    __floordiv__: _UnsignedIntOp
+    __rfloordiv__: _UnsignedIntOp
+    __pow__: _UnsignedIntOp
+    __rpow__: _UnsignedIntOp
 
 class uint8(unsignedinteger):
     def __init__(self, __value: _IntValue = ...) -> None: ...
@@ -958,7 +1031,20 @@ class uint64(unsignedinteger):
     def __init__(self, __value: _IntValue = ...) -> None: ...
 
 class inexact(number): ...  # type: ignore
-class floating(inexact): ...  # type: ignore
+
+class floating(inexact):  # type: ignore
+    __add__: _FloatOp
+    __radd__: _FloatOp
+    __sub__: _FloatOp
+    __rsub__: _FloatOp
+    __mul__: _FloatOp
+    __rmul__: _FloatOp
+    __truediv__: _FloatOp
+    __rtruediv__: _FloatOp
+    __floordiv__: _FloatOp
+    __rfloordiv__: _FloatOp
+    __pow__: _FloatOp
+    __rpow__: _FloatOp
 
 _FloatType = TypeVar('_FloatType', bound=floating)
 
@@ -977,6 +1063,18 @@ class complexfloating(inexact, Generic[_FloatType]):  # type: ignore
     @property
     def imag(self) -> _FloatType: ...  # type: ignore[override]
     def __abs__(self) -> _FloatType: ...  # type: ignore[override]
+    __add__: _ComplexOp
+    __radd__: _ComplexOp
+    __sub__: _ComplexOp
+    __rsub__: _ComplexOp
+    __mul__: _ComplexOp
+    __rmul__: _ComplexOp
+    __truediv__: _ComplexOp
+    __rtruediv__: _ComplexOp
+    __floordiv__: _ComplexOp
+    __rfloordiv__: _ComplexOp
+    __pow__: _ComplexOp
+    __rpow__: _ComplexOp
 
 class complex64(complexfloating[float32]):
     def __init__(self, __value: _ComplexValue = ...) -> None: ...
@@ -987,7 +1085,7 @@ class complex128(complexfloating[float64], complex):
 class flexible(generic): ...  # type: ignore
 
 class void(flexible):
-    def __init__(self, __value: Union[int, integer, bool_, bytes]): ...
+    def __init__(self, __value: Union[_IntLike, _BoolLike, bytes]): ...
     @property
     def real(self: _ArraySelf) -> _ArraySelf: ...
     @property
@@ -995,6 +1093,9 @@ class void(flexible):
 
 class character(flexible): ...  # type: ignore
 
+# NOTE: Most `np.bytes_` / `np.str_` methods return their
+# builtin `bytes` / `str` counterpart
+
 class bytes_(character, bytes):
     @overload
     def __init__(self, __value: object = ...) -> None: ...
@@ -1396,7 +1497,3 @@ def sctype2char(sctype: object) -> str: ...
 def find_common_type(
     array_types: Sequence[DtypeLike], scalar_types: Sequence[DtypeLike]
 ) -> dtype: ...
-
-_NumberLike = Union[int, float, complex, number, bool_]
-_IntLike = Union[int, integer]
-_BoolLike = Union[bool, bool_]
diff --git a/numpy/typing/_callable.py b/numpy/typing/_callable.py
new file mode 100644
index 000000000..5e14b708f
--- /dev/null
+++ b/numpy/typing/_callable.py
@@ -0,0 +1,136 @@
+"""
+A module with various ``typing.Protocol`` subclasses that implement
+the ``__call__`` magic method.
+
+See the `Mypy documentation`_ on protocols for more details.
+
+.. _`Mypy documentation`: https://mypy.readthedocs.io/en/stable/protocols.html#callback-protocols
+
+"""
+
+import sys
+from typing import Union, TypeVar, overload, Any
+
+from numpy import (
+    _BoolLike,
+    _IntLike,
+    _FloatLike,
+    _ComplexLike,
+    _NumberLike,
+    generic,
+    bool_,
+    timedelta64,
+    number,
+    integer,
+    unsignedinteger,
+    signedinteger,
+    int32,
+    int64,
+    floating,
+    float32,
+    float64,
+    complexfloating,
+    complex128,
+)
+
+if sys.version_info >= (3, 8):
+    from typing import Protocol
+    HAVE_PROTOCOL = True
+else:
+    try:
+        from typing_extensions import Protocol
+    except ImportError:
+        HAVE_PROTOCOL = False
+    else:
+        HAVE_PROTOCOL = True
+
+if HAVE_PROTOCOL:
+    _NumberType = TypeVar("_NumberType", bound=number)
+    _NumberType_co = TypeVar("_NumberType_co", covariant=True, bound=number)
+    _GenericType_co = TypeVar("_GenericType_co", covariant=True, bound=generic)
+
+    class _BoolOp(Protocol[_GenericType_co]):
+        @overload
+        def __call__(self, __other: _BoolLike) -> _GenericType_co: ...
+        @overload  # platform dependent
+        def __call__(self, __other: int) -> Union[int32, int64]: ...
+        @overload
+        def __call__(self, __other: float) -> float64: ...
+        @overload
+        def __call__(self, __other: complex) -> complex128: ...
+        @overload
+        def __call__(self, __other: _NumberType) -> _NumberType: ...
+
+    class _BoolSub(Protocol):
+        # Note that `__other: bool_` is absent here
+        @overload  # platform dependent
+        def __call__(self, __other: int) -> Union[int32, int64]: ...
+        @overload
+        def __call__(self, __other: float) -> float64: ...
+        @overload
+        def __call__(self, __other: complex) -> complex128: ...
+        @overload
+        def __call__(self, __other: _NumberType) -> _NumberType: ...
+
+    class _BoolTrueDiv(Protocol):
+        @overload
+        def __call__(self, __other: Union[float, _IntLike, _BoolLike]) -> float64: ...
+        @overload
+        def __call__(self, __other: complex) -> complex128: ...
+        @overload
+        def __call__(self, __other: _NumberType) -> _NumberType: ...
+
+    class _TD64Div(Protocol[_NumberType_co]):
+        @overload
+        def __call__(self, __other: timedelta64) -> _NumberType_co: ...
+        @overload
+        def __call__(self, __other: _FloatLike) -> timedelta64: ...
+
+    class _IntTrueDiv(Protocol):
+        @overload
+        def __call__(self, __other: Union[_IntLike, float]) -> floating: ...
+        @overload
+        def __call__(self, __other: complex) -> complexfloating[floating]: ...
+
+    class _UnsignedIntOp(Protocol):
+        # NOTE: `uint64 + signedinteger -> float64`
+        @overload
+        def __call__(self, __other: Union[bool, unsignedinteger]) -> unsignedinteger: ...
+        @overload
+        def __call__(self, __other: Union[int, signedinteger]) -> Union[signedinteger, float64]: ...
+        @overload
+        def __call__(self, __other: float) -> floating: ...
+        @overload
+        def __call__(self, __other: complex) -> complexfloating[floating]: ...
+
+    class _SignedIntOp(Protocol):
+        @overload
+        def __call__(self, __other: Union[int, signedinteger]) -> signedinteger: ...
+        @overload
+        def __call__(self, __other: float) -> floating: ...
+        @overload
+        def __call__(self, __other: complex) -> complexfloating[floating]: ...
+
+    class _FloatOp(Protocol):
+        @overload
+        def __call__(self, __other: _FloatLike) -> floating: ...
+        @overload
+        def __call__(self, __other: complex) -> complexfloating[floating]: ...
+
+    class _ComplexOp(Protocol):
+        def __call__(self, __other: _ComplexLike) -> complexfloating[floating]: ...
+
+    class _NumberOp(Protocol):
+        def __call__(self, __other: _NumberLike) -> number: ...
+
+else:
+    _BoolOp = Any
+    _BoolSub = Any
+    _BoolTrueDiv = Any
+    _TD64Div = Any
+    _IntTrueDiv = Any
+    _UnsignedIntOp = Any
+    _SignedIntOp = Any
+    _FloatOp = Any
+    _ComplexOp = Any
+    _NumberOp = Any
diff --git a/numpy/typing/tests/data/fail/arithmetic.py b/numpy/typing/tests/data/fail/arithmetic.py
new file mode 100644
index 000000000..169e104f9
--- /dev/null
+++ b/numpy/typing/tests/data/fail/arithmetic.py
@@ -0,0 +1,19 @@
+import numpy as np
+
+b_ = np.bool_()
+dt = np.datetime64(0, "D")
+td = np.timedelta64(0, "D")
+
+b_ - b_  # E: No overload variant
+
+dt + dt  # E: Unsupported operand types
+td - dt  # E: Unsupported operand types
+td % 1  # E: Unsupported operand types
+td / dt  # E: No overload
+
+# NOTE: The 1 tests below currently don't work due to the broad
+# (i.e. untyped) signature of `.__mod__()`.
+# TODO: Revisit this once annotations are added to the
+# `_ArrayOrScalarCommon` magic methods.
+
+# td % dt  # E: Unsupported operand types
diff --git a/numpy/typing/tests/data/fail/scalars.py b/numpy/typing/tests/data/fail/scalars.py
index 47c031163..13bb45483 100644
--- a/numpy/typing/tests/data/fail/scalars.py
+++ b/numpy/typing/tests/data/fail/scalars.py
@@ -28,22 +28,6 @@ np.complex64(1, 2)  # E: Too many arguments
 
 np.datetime64(0)  # E: non-matching overload
 
-dt_64 = np.datetime64(0, "D")
-td_64 = np.timedelta64(1, "h")
-
-dt_64 + dt_64  # E: Unsupported operand types
-td_64 - dt_64  # E: Unsupported operand types
-td_64 % 1  # E: Unsupported operand types
-
-# NOTE: The 2 tests below currently don't work due to the broad
-# (i.e. untyped) signature of `generic.__truediv__()` and `.__mod__()`.
-# TODO: Revisit this once annotations are added to the
-# `_ArrayOrScalarCommon` magic methods.
-
-# td_64 / dt_64  # E: No overload
-# td_64 % dt_64  # E: Unsupported operand types
-
-
 class A:
     def __float__(self):
         return 1.0
diff --git a/numpy/typing/tests/data/pass/arithmetic.py b/numpy/typing/tests/data/pass/arithmetic.py
new file mode 100644
index 000000000..f26eab879
--- /dev/null
+++ b/numpy/typing/tests/data/pass/arithmetic.py
@@ -0,0 +1,257 @@
+import numpy as np
+
+c16 = np.complex128(1)
+f8 = np.float64(1)
+i8 = np.int64(1)
+u8 = np.uint64(1)
+
+c8 = np.complex64(1)
+f4 = np.float32(1)
+i4 = np.int32(1)
+u4 = np.uint32(1)
+
+dt = np.datetime64(1, "D")
+td = np.timedelta64(1, "D")
+
+b_ = np.bool_(1)
+
+b = bool(1)
+c = complex(1)
+f = float(1)
+i = int(1)
+
+AR = np.ones(1, dtype=np.float64)
+AR.setflags(write=False)
+
+# Time structures
+
+dt + td
+dt + i
+dt + i4
+dt + i8
+dt - dt
+dt - i
+dt - i4
+dt - i8
+
+td + td
+td + i
+td + i4
+td + i8
+td - td
+td - i
+td - i4
+td - i8
+td / f
+td / f4
+td / f8
+td / td
+td // td
+td % td
+
+
+# boolean
+
+b_ / b
+b_ / b_
+b_ / i
+b_ / i8
+b_ / i4
+b_ / u8
+b_ / u4
+b_ / f
+b_ / f8
+b_ / f4
+b_ / c
+b_ / c16
+b_ / c8
+
+b / b_
+b_ / b_
+i / b_
+i8 / b_
+i4 / b_
+u8 / b_
+u4 / b_
+f / b_
+f8 / b_
+f4 / b_
+c / b_
+c16 / b_
+c8 / b_
+
+# Complex
+
+c16 + c16
+c16 + f8
+c16 + i8
+c16 + c8
+c16 + f4
+c16 + i4
+c16 + b_
+c16 + b
+c16 + c
+c16 + f
+c16 + i
+c16 + AR
+
+c16 + c16
+f8 + c16
+i8 + c16
+c8 + c16
+f4 + c16
+i4 + c16
+b_ + c16
+b + c16
+c + c16
+f + c16
+i + c16
+AR + c16
+
+c8 + c16
+c8 + f8
+c8 + i8
+c8 + c8
+c8 + f4
+c8 + i4
+c8 + b_
+c8 + b
+c8 + c
+c8 + f
+c8 + i
+c8 + AR
+
+c16 + c8
+f8 + c8
+i8 + c8
+c8 + c8
+f4 + c8
+i4 + c8
+b_ + c8
+b + c8
+c + c8
+f + c8
+i + c8
+AR + c8
+
+# Float
+
+f8 + f8
+f8 + i8
+f8 + f4
+f8 + i4
+f8 + b_
+f8 + b
+f8 + c
+f8 + f
+f8 + i
+f8 + AR
+
+f8 + f8
+i8 + f8
+f4 + f8
+i4 + f8
+b_ + f8
+b + f8
+c + f8
+f + f8
+i + f8
+AR + f8
+
+f4 + f8
+f4 + i8
+f4 + f4
+f4 + i4
+f4 + b_
+f4 + b
+f4 + c
+f4 + f
+f4 + i
+f4 + AR
+
+f8 + f4
+i8 + f4
+f4 + f4
+i4 + f4
+b_ + f4
+b + f4
+c + f4
+f + f4
+i + f4
+AR + f4
+
+# Int
+
+i8 + i8
+i8 + u8
+i8 + i4
+i8 + u4
+i8 + b_
+i8 + b
+i8 + c
+i8 + f
+i8 + i
+i8 + AR
+
+u8 + u8
+u8 + i4
+u8 + u4
+u8 + b_
+u8 + b
+u8 + c
+u8 + f
+u8 + i
+u8 + AR
+
+i8 + i8
+u8 + i8
+i4 + i8
+u4 + i8
+b_ + i8
+b + i8
+c + i8
+f + i8
+i + i8
+AR + i8
+
+u8 + u8
+i4 + u8
+u4 + u8
+b_ + u8
+b + u8
+c + u8
+f + u8
+i + u8
+AR + u8
+
+i4 + i8
+i4 + i4
+i4 + i
+i4 + b_
+i4 + b
+i4 + AR
+
+u4 + i8
+u4 + i4
+u4 + u8
+u4 + u4
+u4 + i
+u4 + b_
+u4 + b
+u4 + AR
+
+i8 + i4
+i4 + i4
+i + i4
+b_ + i4
+b + i4
+AR + i4
+
+i8 + u4
+i4 + u4
+u8 + u4
+u4 + u4
+b_ + u4
+b + u4
+i + u4
+AR + u4
diff --git a/numpy/typing/tests/data/pass/scalars.py b/numpy/typing/tests/data/pass/scalars.py
index c02e1ed36..49ddb8ed9 100644
--- a/numpy/typing/tests/data/pass/scalars.py
+++ b/numpy/typing/tests/data/pass/scalars.py
@@ -108,19 +108,6 @@ np.timedelta64(dt.timedelta(2))
 np.timedelta64(None)
 np.timedelta64(None, "D")
 
-dt_64 = np.datetime64(0, "D")
-td_64 = np.timedelta64(1, "h")
-
-dt_64 + td_64
-dt_64 - dt_64
-dt_64 - td_64
-
-td_64 + td_64
-td_64 - td_64
-td_64 / 1.0
-td_64 / td_64
-td_64 % td_64
-
 np.void(1)
 np.void(np.int64(1))
 np.void(True)
diff --git a/numpy/typing/tests/data/pass/ufuncs.py b/numpy/typing/tests/data/pass/ufuncs.py
index 82172952a..ad4d483d4 100644
--- a/numpy/typing/tests/data/pass/ufuncs.py
+++ b/numpy/typing/tests/data/pass/ufuncs.py
@@ -6,7 +6,10 @@ np.sin(1, out=np.empty(1))
 np.matmul(np.ones((2, 2, 2)), np.ones((2, 2, 2)), axes=[(0, 1), (0, 1), (0, 1)])
 np.sin(1, signature="D")
 np.sin(1, extobj=[16, 1, lambda: None])
-np.sin(1) + np.sin(1)
+# NOTE: `np.generic` subclasses are not guaranteed to support addition;
+# re-enable this we can infer the exact return type of `np.sin(...)`.
+#
+# np.sin(1) + np.sin(1)
 np.sin.types[0]
 np.sin.__name__
 
diff --git a/numpy/typing/tests/data/reveal/arithmetic.py b/numpy/typing/tests/data/reveal/arithmetic.py
new file mode 100644
index 000000000..b8c457aaf
--- /dev/null
+++ b/numpy/typing/tests/data/reveal/arithmetic.py
@@ -0,0 +1,256 @@
+import numpy as np
+
+c16 = np.complex128()
+f8 = np.float64()
+i8 = np.int64()
+u8 = np.uint64()
+
+c8 = np.complex64()
+f4 = np.float32()
+i4 = np.int32()
+u4 = np.uint32()
+
+dt = np.datetime64(0, "D")
+td = np.timedelta64(0, "D")
+
+b_ = np.bool_()
+
+b = bool()
+c = complex()
+f = float()
+i = int()
+
+AR = np.array([0], dtype=np.float64)
+AR.setflags(write=False)
+
+# Time structures
+
+reveal_type(dt + td)  # E: numpy.datetime64
+reveal_type(dt + i)  # E: numpy.datetime64
+reveal_type(dt + i4)  # E: numpy.datetime64
+reveal_type(dt + i8)  # E: numpy.datetime64
+reveal_type(dt - dt)  # E: numpy.timedelta64
+reveal_type(dt - i)  # E: numpy.datetime64
+reveal_type(dt - i4)  # E: numpy.datetime64
+reveal_type(dt - i8)  # E: numpy.datetime64
+
+reveal_type(td + td)  # E: numpy.timedelta64
+reveal_type(td + i)  # E: numpy.timedelta64
+reveal_type(td + i4)  # E: numpy.timedelta64
+reveal_type(td + i8)  # E: numpy.timedelta64
+reveal_type(td - td)  # E: numpy.timedelta64
+reveal_type(td - i)  # E: numpy.timedelta64
+reveal_type(td - i4)  # E: numpy.timedelta64
+reveal_type(td - i8)  # E: numpy.timedelta64
+reveal_type(td / f)  # E: numpy.timedelta64
+reveal_type(td / f4)  # E: numpy.timedelta64
+reveal_type(td / f8)  # E: numpy.timedelta64
+reveal_type(td / td)  # E: float64
+reveal_type(td // td)  # E: signedinteger
+reveal_type(td % td)  # E: numpy.timedelta64
+
+# boolean
+
+reveal_type(b_ / b)  # E: float64
+reveal_type(b_ / b_)  # E: float64
+reveal_type(b_ / i)  # E: float64
+reveal_type(b_ / i8)  # E: float64
+reveal_type(b_ / i4)  # E: float64
+reveal_type(b_ / u8)  # E: float64
+reveal_type(b_ / u4)  # E: float64
+reveal_type(b_ / f)  # E: float64
+reveal_type(b_ / f8)  # E: float64
+reveal_type(b_ / f4)  # E: float32
+reveal_type(b_ / c)  # E: complex128
+reveal_type(b_ / c16)  # E: complex128
+reveal_type(b_ / c8)  # E: complex64
+
+reveal_type(b / b_)  # E: float64
+reveal_type(b_ / b_)  # E: float64
+reveal_type(i / b_)  # E: float64
+reveal_type(i8 / b_)  # E: float64
+reveal_type(i4 / b_)  # E: float64
+reveal_type(u8 / b_)  # E: float64
+reveal_type(u4 / b_)  # E: float64
+reveal_type(f / b_)  # E: float64
+reveal_type(f8 / b_)  # E: float64
+reveal_type(f4 / b_)  # E: float32
+reveal_type(c / b_)  # E: complex128
+reveal_type(c16 / b_)  # E: complex128
+reveal_type(c8 / b_)  # E: complex64
+
+# Complex
+
+reveal_type(c16 + c16)  # E: complexfloating
+reveal_type(c16 + f8)  # E: complexfloating
+reveal_type(c16 + i8)  # E: complexfloating
+reveal_type(c16 + c8)  # E: complexfloating
+reveal_type(c16 + f4)  # E: complexfloating
+reveal_type(c16 + i4)  # E: complexfloating
+reveal_type(c16 + b_)  # E: complex128
+reveal_type(c16 + b)  # E: complexfloating
+reveal_type(c16 + c)  # E: complexfloating
+reveal_type(c16 + f)  # E: complexfloating
+reveal_type(c16 + i)  # E: complexfloating
+reveal_type(c16 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(c16 + c16)  # E: complexfloating
+reveal_type(f8 + c16)  # E: complexfloating
+reveal_type(i8 + c16)  # E: complexfloating
+reveal_type(c8 + c16)  # E: complexfloating
+reveal_type(f4 + c16)  # E: complexfloating
+reveal_type(i4 + c16)  # E: complexfloating
+reveal_type(b_ + c16)  # E: complex128
+reveal_type(b + c16)  # E: complexfloating
+reveal_type(c + c16)  # E: complexfloating
+reveal_type(f + c16)  # E: complexfloating
+reveal_type(i + c16)  # E: complexfloating
+reveal_type(AR + c16)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(c8 + c16)  # E: complexfloating
+reveal_type(c8 + f8)  # E: complexfloating
+reveal_type(c8 + i8)  # E: complexfloating
+reveal_type(c8 + c8)  # E: complexfloating
+reveal_type(c8 + f4)  # E: complexfloating
+reveal_type(c8 + i4)  # E: complexfloating
+reveal_type(c8 + b_)  # E: complex64
+reveal_type(c8 + b)  # E: complexfloating
+reveal_type(c8 + c)  # E: complexfloating
+reveal_type(c8 + f)  # E: complexfloating
+reveal_type(c8 + i)  # E: complexfloating
+reveal_type(c8 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(c16 + c8)  # E: complexfloating
+reveal_type(f8 + c8)  # E: complexfloating
+reveal_type(i8 + c8)  # E: complexfloating
+reveal_type(c8 + c8)  # E: complexfloating
+reveal_type(f4 + c8)  # E: complexfloating
+reveal_type(i4 + c8)  # E: complexfloating
+reveal_type(b_ + c8)  # E: complex64
+reveal_type(b + c8)  # E: complexfloating
+reveal_type(c + c8)  # E: complexfloating
+reveal_type(f + c8)  # E: complexfloating
+reveal_type(i + c8)  # E: complexfloating
+reveal_type(AR + c8)  # E: Union[numpy.ndarray, numpy.generic]
+
+# Float
+
+reveal_type(f8 + f8)  # E: floating
+reveal_type(f8 + i8)  # E: floating
+reveal_type(f8 + f4)  # E: floating
+reveal_type(f8 + i4)  # E: floating
+reveal_type(f8 + b_)  # E: float64
+reveal_type(f8 + b)  # E: floating
+reveal_type(f8 + c)  # E: complexfloating
+reveal_type(f8 + f)  # E: floating
+reveal_type(f8 + i)  # E: floating
+reveal_type(f8 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(f8 + f8)  # E: floating
+reveal_type(i8 + f8)  # E: floating
+reveal_type(f4 + f8)  # E: floating
+reveal_type(i4 + f8)  # E: floating
+reveal_type(b_ + f8)  # E: float64
+reveal_type(b + f8)  # E: floating
+reveal_type(c + f8)  # E: complexfloating
+reveal_type(f + f8)  # E: floating
+reveal_type(i + f8)  # E: floating
+reveal_type(AR + f8)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(f4 + f8)  # E: floating
+reveal_type(f4 + i8)  # E: floating
+reveal_type(f4 + f4)  # E: floating
+reveal_type(f4 + i4)  # E: floating
+reveal_type(f4 + b_)  # E: float32
+reveal_type(f4 + b)  # E: floating
+reveal_type(f4 + c)  # E: complexfloating
+reveal_type(f4 + f)  # E: floating
+reveal_type(f4 + i)  # E: floating
+reveal_type(f4 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(f8 + f4)  # E: floating
+reveal_type(i8 + f4)  # E: floating
+reveal_type(f4 + f4)  # E: floating
+reveal_type(i4 + f4)  # E: floating
+reveal_type(b_ + f4)  # E: float32
+reveal_type(b + f4)  # E: floating
+reveal_type(c + f4)  # E: complexfloating
+reveal_type(f + f4)  # E: floating
+reveal_type(i + f4)  # E: floating
+reveal_type(AR + f4)  # E: Union[numpy.ndarray, numpy.generic]
+
+# Int
+
+reveal_type(i8 + i8)  # E: signedinteger
+reveal_type(i8 + u8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(i8 + i4)  # E: signedinteger
+reveal_type(i8 + u4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(i8 + b_)  # E: int64
+reveal_type(i8 + b)  # E: signedinteger
+reveal_type(i8 + c)  # E: complexfloating
+reveal_type(i8 + f)  # E: floating
+reveal_type(i8 + i)  # E: signedinteger
+reveal_type(i8 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(u8 + u8)  # E: unsignedinteger
+reveal_type(u8 + i4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u8 + u4)  # E: unsignedinteger
+reveal_type(u8 + b_)  # E: uint64
+reveal_type(u8 + b)  # E: unsignedinteger
+reveal_type(u8 + c)  # E: complexfloating
+reveal_type(u8 + f)  # E: floating
+reveal_type(u8 + i)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u8 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(i8 + i8)  # E: signedinteger
+reveal_type(u8 + i8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(i4 + i8)  # E: signedinteger
+reveal_type(u4 + i8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(b_ + i8)  # E: int64
+reveal_type(b + i8)  # E: signedinteger
+reveal_type(c + i8)  # E: complexfloating
+reveal_type(f + i8)  # E: floating
+reveal_type(i + i8)  # E: signedinteger
+reveal_type(AR + i8)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(u8 + u8)  # E: unsignedinteger
+reveal_type(i4 + u8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u4 + u8)  # E: unsignedinteger
+reveal_type(b_ + u8)  # E: uint64
+reveal_type(b + u8)  # E: unsignedinteger
+reveal_type(c + u8)  # E: complexfloating
+reveal_type(f + u8)  # E: floating
+reveal_type(i + u8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(AR + u8)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(i4 + i8)  # E: signedinteger
+reveal_type(i4 + i4)  # E: signedinteger
+reveal_type(i4 + i)  # E: signedinteger
+reveal_type(i4 + b_)  # E: int32
+reveal_type(i4 + b)  # E: signedinteger
+reveal_type(i4 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(u4 + i8)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u4 + i4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u4 + u8)  # E: unsignedinteger
+reveal_type(u4 + u4)  # E: unsignedinteger
+reveal_type(u4 + i)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u4 + b_)  # E: uint32
+reveal_type(u4 + b)  # E: unsignedinteger
+reveal_type(u4 + AR)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(i8 + i4)  # E: signedinteger
+reveal_type(i4 + i4)  # E: signedinteger
+reveal_type(i + i4)  # E: signedinteger
+reveal_type(b_ + i4)  # E: int32
+reveal_type(b + i4)  # E: signedinteger
+reveal_type(AR + i4)  # E: Union[numpy.ndarray, numpy.generic]
+
+reveal_type(i8 + u4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(i4 + u4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(u8 + u4)  # E: unsignedinteger
+reveal_type(u4 + u4)  # E: unsignedinteger
+reveal_type(b_ + u4)  # E: uint32
+reveal_type(b + u4)  # E: unsignedinteger
+reveal_type(i + u4)  # E: Union[numpy.signedinteger, numpy.float64]
+reveal_type(AR + u4)  # E: Union[numpy.ndarray, numpy.generic]
diff --git a/numpy/typing/tests/data/reveal/scalars.py b/numpy/typing/tests/data/reveal/scalars.py
index 882fe9612..ec3713b0f 100644
--- a/numpy/typing/tests/data/reveal/scalars.py
+++ b/numpy/typing/tests/data/reveal/scalars.py
@@ -12,22 +12,5 @@ reveal_type(x.itemsize)  # E: int
 reveal_type(x.shape)  # E: tuple[builtins.int]
 reveal_type(x.strides)  # E: tuple[builtins.int]
 
-# Time structures
-dt = np.datetime64(0, "D")
-td = np.timedelta64(0, "D")
-
-reveal_type(dt + td)  # E: numpy.datetime64
-reveal_type(dt + 1)  # E: numpy.datetime64
-reveal_type(dt - dt)  # E: numpy.timedelta64
-reveal_type(dt - 1)  # E: numpy.timedelta64
-
-reveal_type(td + td)  # E: numpy.timedelta64
-reveal_type(td + 1)  # E: numpy.timedelta64
-reveal_type(td - td)  # E: numpy.timedelta64
-reveal_type(td - 1)  # E: numpy.timedelta64
-reveal_type(td / 1.0)  # E: numpy.timedelta64
-reveal_type(td / td)  # E: float
-reveal_type(td % td)  # E: numpy.timedelta64
-
 reveal_type(np.complex64().real)  # E: numpy.float32
 reveal_type(np.complex128().imag)  # E: numpy.float64