eriknw · Nov 13, 2026 · Nov 3, 2026 · Nov 4, 2026 · Nov 4, 2026 · Nov 8, 2026
diff --git a/graphblas/core/matrix.py b/graphblas/core/matrix.py
@@ -459,7 +459,8 @@ def isequal(self, other, *, check_dtype=False, **opts):
     def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False, **opts):
         """Check for approximate equality (including same size and same structure).
 
-        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``.
+        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``,
+        except that, as in ``math.isclose``, an infinity is close only to itself.
 
         Parameters
         ----------

diff --git a/graphblas/core/operator/binary.py b/graphblas/core/operator/binary.py
@@ -339,11 +339,51 @@ def _rpow(x, y):
 
 def _isclose(rel_tol=1e-7, abs_tol=0.0):
     def inner(x, y):  # pragma: no cover (numba)
-        return x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol)
+        if x == y:
+            return True
+        # As in math.isclose, an infinity is close only to itself
+        if abs(x) == np.inf or abs(y) == np.inf:
+            return False
+        return _abs_diff(x, y) <= max(rel_tol * max(_magnitude(x), _magnitude(y)), abs_tol)
 
     return inner
 
 
+# Integer arithmetic in Numba wraps: x - y in int64 or uint64 (Numba widens smaller
+# unsigned types to uint64, so 5 - 6 is 2**64 - 1), and abs at a signed type's
+# minimum. Python ints do not, so for them these are plain Python; the overloads
+# below compile them for Numba's integers so that they do not wrap either.
+def _abs_diff(x, y):
+    return abs(x - y)
+
+
+def _magnitude(x):
+    return abs(x)
+
+
+if _has_numba:
+    from numba.extending import overload
+
+    @overload(_abs_diff)
+    def _abs_diff_overload(x, y):
+        if isinstance(x, numba.types.Integer) and x == y:
+
+            def abs_diff(x, y):  # pragma: no cover (numba)
+                # The difference fits uint64, where both convert modulo 2**64
+                if x >= y:
+                    return np.uint64(x) - np.uint64(y)
+                return np.uint64(y) - np.uint64(x)
+
+            return abs_diff
+        return lambda x, y: abs(x - y)
+
+    @overload(_magnitude)
+    def _magnitude_overload(x):
+        if isinstance(x, numba.types.Integer):
+            return lambda x: abs(np.float64(x))
+        return lambda x: abs(x)  # noqa: PLW0108 (an overload returns a Python function)
+
+
 _MAX_INT64 = np.iinfo(np.int64).max
 
 
@@ -843,11 +883,16 @@ def register_new(
         >>> dir(gb.binary)
         [..., 'max_zero', ...]
 
-        This is how ``gb.binary.isclose`` is defined:
+        ``gb.binary.isclose`` is defined like this, except that it computes ``x - y``
+        and ``abs(x)`` for integers without wrapping:
 
         >>> def user_isclose(rel_tol=1e-7, abs_tol=0.0):
         >>>     def inner(x, y):
-        >>>         return x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol)
+        >>>         if x == y:
+        >>>             return True
+        >>>         if abs(x) == np.inf or abs(y) == np.inf:
+        >>>             return False
+        >>>         return abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol)
         >>>     return inner
         >>> gb.binary.register_new("user_isclose", user_isclose, parameterized=True)
 

diff --git a/graphblas/core/operator/unary.py b/graphblas/core/operator/unary.py
@@ -319,7 +319,8 @@ def register_anonymous(cls, func, name=None, *, parameterized=False, is_udt=Fals
             can be used as: ``A.ewise_mult(B, gb.binary.isclose(rel_tol=1e-5))``.
             When creating a parameterized user-defined operator, the ``func``
             parameter must be a callable that *returns* a function that will
-            then get compiled. See the ``user_isclose`` example below.
+            then get compiled. See the ``user_isclose`` example in
+            ``BinaryOp.register_new``.
         is_udt : bool, default False
             Whether the operator is intended to operate on user-defined types.
             If True, then the function will not be automatically compiled for
@@ -357,7 +358,8 @@ def register_new(cls, name, func, *, parameterized=False, is_udt=False, lazy=Fal
             can be used as: ``A.ewise_mult(B, gb.binary.isclose(rel_tol=1e-5))``.
             When creating a parameterized user-defined operator, the ``func``
             parameter must be a callable that *returns* a function that will
-            then get compiled. See the ``user_isclose`` example below.
+            then get compiled. See the ``user_isclose`` example in
+            ``BinaryOp.register_new``.
         is_udt : bool, default False
             Whether the operator is intended to operate on user-defined types.
             If True, then the function will not be automatically compiled for

diff --git a/graphblas/core/scalar.py b/graphblas/core/scalar.py
@@ -17,6 +17,7 @@
     _literal_type,
 )
 from .operator import get_typed_op
+from .operator.binary import _isclose
 from .udt import (
     _element_bytes,
     _element_shapes_differ,
@@ -34,7 +35,7 @@
 if _supports_udfs:
     from ..binary import isclose
 else:
-    from .operator.binary import _isclose as isclose
+    isclose = _isclose
 
 ffi_new = ffi.new
 
@@ -271,7 +272,8 @@ def isequal(self, other, *, check_dtype=False):
     def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False):
         """Check for approximate equality (including whether the value is missing).
 
-        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``.
+        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``,
+        except that, as in ``math.isclose``, an infinity is close only to itself.
 
         Parameters
         ----------
@@ -318,16 +320,19 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False):
         # We can't yet call a UDF on a scalar as part of the spec, so let's do it ourselves
         isclose_func = isclose(rel_tol, abs_tol)
         if not _has_numba:
-            # Check if types are compatible
-            get_typed_op(
+            # Check if types are compatible, and compare in the type they unify to, as
+            # with Numba (INT64 with UINT64 is in FP64). Python ints do not wrap in x - y.
+            eq = get_typed_op(
                 binary.eq,
                 self.dtype,
                 other.dtype,
                 is_left_scalar=True,
                 is_right_scalar=True,
                 kind="binary",
             )
-            return isclose_func(self.value, other.value)
+            x = eq.type.np_type.type(self.value).item()
+            y = eq.type2.np_type.type(other.value).item()
+            return isclose_func(x, y)
         isclose_func = get_typed_op(
             isclose_func,
             self.dtype,
@@ -336,7 +341,11 @@ def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False):
             is_right_scalar=True,
             kind="binary",
         )
-        return isclose_func._numba_func(self.value, other.value)
+        # Pass the values in the op's input types, as ewise_mult does; Numba types a
+        # Python int as int64, which a UINT64 value need not fit
+        x = isclose_func.type.np_type.type(self.value)
+        y = isclose_func.type2.np_type.type(other.value)
+        return isclose_func._numba_func(x, y)
 
     def clear(self):
         """In-place operation which clears the value in the Scalar.

diff --git a/graphblas/core/ss/matrix.py b/graphblas/core/ss/matrix.py
@@ -7,7 +7,13 @@
 
 from ... import binary, monoid
 from ...dtypes import _INDEX, BOOL, INT64, UINT64, lookup_dtype
-from ...exceptions import InvalidValue, _error_code_lookup, check_status, check_status_carg
+from ...exceptions import (
+    IndexOutOfBound,
+    InvalidValue,
+    _error_code_lookup,
+    check_status,
+    check_status_carg,
+)
 from .. import NULL, _has_numba, ffi, lib
 from ..base import call
 from ..operator import get_typed_op
@@ -863,7 +869,9 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m
             return rv
 
         if method == "export":
-            mhandle = ffi_new("GrB_Matrix*", parent._carg)
+            # GraphBLAS frees the Matrix as it hands the arrays over, and sets this
+            # handle to NULL, so an error after that cannot leave it dangling
+            mhandle = parent.gb_obj
             type_ = ffi_new("GrB_Type*")
             nrows_ = ffi_new("GrB_Index*")
             ncols_ = ffi_new("GrB_Index*")
@@ -1170,8 +1178,6 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m
         rv["is_iso"] = is_iso
         rv["format"] = format
         rv["values"] = values
-        if method == "export":
-            parent.gb_obj[0] = NULL
         if parent.dtype._is_udt:
             rv["dtype"] = parent.dtype
         return rv
@@ -3171,7 +3177,7 @@ def _import_coor(
             raise ValueError(f"Invalid format: {format!r}.  Must be None or 'coor'.")
         if not sorted_rows:
             raise ValueError("sorted_rows must be True when importing 'coor' format")
-        indptr = indices_to_indptr(rows, nrows + 1)
+        indptr = _coo_indptr(rows, nrows, "rows", "nrows")
         return cls._import_csr(
             nrows=nrows,
             ncols=ncols,
@@ -3344,7 +3350,7 @@ def _import_cooc(
             raise ValueError(f"Invalid format: {format!r}.  Must be None or 'cooc'.")
         if not sorted_cols:
             raise ValueError("sorted_cols must be True when importing 'cooc' format")
-        indptr = indices_to_indptr(cols, ncols + 1)
+        indptr = _coo_indptr(cols, ncols, "cols", "ncols")
         return cls._import_csc(
             nrows=nrows,
             ncols=ncols,
@@ -4558,16 +4564,36 @@ def issorted(arr):  # pragma: no cover (numba)
 
 @njit
 def indices_to_indptr(indices, size):  # pragma: no cover (numba)
-    """Calculate the indptr for e.g. CSR from sorted COO rows."""
+    """Calculate the indptr, e.g. for CSR, from uint64 COO rows, each below ``size - 1``.
+
+    Each row's count goes at ``indptr[row + 1]``, then a running sum makes the offsets,
+    so every empty row, wherever it is, gets an empty range.
+    """
     indptr = np.zeros(size, dtype=indices.dtype)
-    index = np.uint64(0)
     one = np.uint64(1)
     for i in range(indices.size):
-        row = indices[i]
-        if row != index:
-            indptr[index + one] = i
-            index = row
-    indptr[index + one] = indices.size
+        indptr[indices[i] + one] += one
+    for i in range(1, size):
+        indptr[i] += indptr[i - 1]
+    return indptr
+
+
+def _coo_indptr(indices, size, name, size_name):
+    """Return the CSR (or CSC) indptr, of length ``size + 1``, for sorted COO rows (or columns)."""
+    # Read 2-D indices as a flat buffer, as other imports do. The kernel would index
+    # a 2-D array by its rows, past the end of its first axis.
+    indices = ints_to_numpy_buffer(indices, np.uint64, name=name).ravel()
+    # indices_to_indptr writes at each index + 1, with no bounds check under Numba
+    if indices.size > 0 and (biggest := indices.max()) >= size:
+        raise IndexOutOfBound(f"{name} must be less than {size_name}, {size}; got {biggest}")
+    if _has_numba:
+        return indices_to_indptr(indices, size + 1)
+    # Without Numba the kernel's loops run in Python, once per entry and once per row.
+    # Every index is below size, so it reads the same as int64; counts are never negative.
+    indptr = np.empty(size + 1, np.uint64)
+    indptr[0] = 0
+    counts = np.bincount(indices.view(np.int64), minlength=size)
+    np.cumsum(counts.view(np.uint64), out=indptr[1:])
     return indptr
 
 

diff --git a/graphblas/core/ss/vector.py b/graphblas/core/ss/vector.py
@@ -515,7 +515,9 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m
 
         size = parent._size
         if method == "export":
-            vhandle = ffi_new("GrB_Vector*", parent._carg)
+            # GraphBLAS frees the Vector as it hands the arrays over, and sets this
+            # handle to NULL, so an error after that cannot leave it dangling
+            vhandle = parent.gb_obj
             type_ = ffi_new("GrB_Type*")
             size_ = ffi_new("GrB_Index*")
             args = (type_, size_)
@@ -642,8 +644,6 @@ def _export(self, format=None, *, sort=False, give_ownership=False, raw=False, m
             format=format,
             values=values,
         )
-        if method == "export":
-            parent.gb_obj[0] = NULL
         if parent.dtype._is_udt:
             rv["dtype"] = parent.dtype
         return rv

diff --git a/graphblas/core/vector.py b/graphblas/core/vector.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 from .. import backend, binary, monoid, select, semiring, unary
-from ..dtypes import _INDEX, FP64, INT64, lookup_dtype, unify
+from ..dtypes import _INDEX, FP64, INT64, UINT64, lookup_dtype, unify
 from ..exceptions import DimensionMismatch, NoValue, check_status
 from . import _supports_udfs, automethods, ffi, lib, utils
 from .base import BaseExpression, BaseType, call
@@ -110,8 +110,13 @@ def _select_mask(updater, obj, mask):
         updater << obj.dup(mask=mask)
 
 
+def _abs_type(dtype):
+    return dtype if dtype.np_type.kind == "c" else FP64
+
+
 def _isclose_recipe(self, other, rel_tol, abs_tol, **opts):
-    #  x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol)
+    #  x == y or abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol), except that
+    #  an infinity is close only to itself
     isequal = self.ewise_mult(other, binary.eq).new(bool, name="isclose", **opts)
     if isequal._nvals != self._nvals:
         return False
@@ -124,23 +129,38 @@ def _isclose_recipe(self, other, rel_tol, abs_tol, **opts):
     # So we can use structural mask below
     isequal(**opts) << select.value(isequal == True)  # noqa: E712
 
-    # abs(x)
-    x = self.apply(unary.abs).new(FP64, mask=~isequal.S, **opts)
-    # abs(y)
-    y = other.apply(unary.abs).new(FP64, mask=~isequal.S, **opts)
+    # abs(x) and abs(y), in the type x - y is in below, or in FP64 for a real type:
+    # abs wraps at an integer type's minimum
+    abs_type = _abs_type(unify(self.dtype, other.dtype))
+    x = self.apply(unary.abs[abs_type]).new(FP64, mask=~isequal.S, **opts)
+    y = other.apply(unary.abs[abs_type]).new(FP64, mask=~isequal.S, **opts)
     # max(abs(x), abs(y))
     x(**opts) << x.ewise_mult(y, binary.max)
     max_x_y = x
+    # As in math.isclose, an infinity is close only to itself, and those are equal
+    if max_x_y.ndim == 1:
+        biggest = max_x_y.reduce(monoid.max, allow_empty=False).new(**opts).value
+    else:
+        biggest = max_x_y.reduce_scalar(monoid.max, allow_empty=False).new(**opts).value
+    if biggest == np.inf:
+        return False
     # rel_tol * max(abs(x), abs(y))
     max_x_y(**opts) << max_x_y.apply(binary.times, rel_tol)
     # max(rel_tol * max(abs(x), abs(y)), abs_tol)
     max_x_y(**opts) << max_x_y.apply(binary.max, abs_tol)
 
-    # x - y
-    y(~isequal.S, replace=True, **opts) << self.ewise_mult(other, binary.minus)
-    abs_x_y = y
     # abs(x - y)
-    abs_x_y(**opts) << abs_x_y.apply(unary.abs)
+    if unify(self.dtype, other.dtype).np_type.kind in "iu":
+        # x - y wraps in an integer type, and the difference may not fit it. It fits
+        # UINT64, where max(x, y) and min(x, y) convert modulo 2**64.
+        hi = self.ewise_mult(other, binary.max).new(UINT64, mask=~isequal.S, **opts)
+        lo = self.ewise_mult(other, binary.min).new(UINT64, mask=~isequal.S, **opts)
+        y(~isequal.S, replace=True, **opts) << hi.ewise_mult(lo, binary.minus)
+    else:
+        # abs before storing in y, whose FP64 would keep only a complex number's real part
+        diff = self.ewise_mult(other, binary.minus).new(mask=~isequal.S, **opts)
+        y(~isequal.S, replace=True, **opts) << diff.apply(unary.abs)
+    abs_x_y = y
 
     # abs(x - y) <= max(rel_tol * max(abs(x), abs(y)), abs_tol)
     isequal(**opts) << abs_x_y.ewise_mult(max_x_y, binary.le)
@@ -398,7 +418,8 @@ def isequal(self, other, *, check_dtype=False, **opts):
     def isclose(self, other, *, rel_tol=1e-7, abs_tol=0.0, check_dtype=False, **opts):
         """Check for approximate equality (including same size and same structure).
 
-        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``.
+        Equivalent to: ``abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)``,
+        except that, as in ``math.isclose``, an infinity is close only to itself.
 
         Parameters
         ----------

diff --git a/graphblas/select/__init__.py b/graphblas/select/__init__.py
@@ -116,7 +116,7 @@ def _match_expr(parent, expr):
         and len(args) == 2
         and op.parent in _binary_to_select
         and args[1]._is_scalar
-        and (args[0] is parent or getattr(args[0], "T", None) is parent)
+        and args[0] is parent
         and op.type in _binary_to_select[op.parent]
     ):
         selectop = _binary_to_select[op.parent][op.type]