Taking discussion out of code: You are right, result_type doesn't work for things like result_type(input-quantity, rtol-float).

I think in the end, perhaps one just has to write more carefully... Looking at this code, really the two different paths for finite and non-finite are not needed if one is a bit more careful, which means one can avoid the broadcasting altogether -- see mdhaber#1 (@mdhaber - I made this PR to your branch so that we can keep credit for both you and the earlier commit).

From that, perhaps it is useful to have a general helper routine that does asanyarray except passes through python scalars? I now do

    x, y, atol, rtol = (
        a if isinstance(a, (int, float, complex)) else asanyarray(a)
        for a in (a, b, atol, rtol))

Though it may be that we simply have to see how things pan out.

p.s. Just for completeness, the NEP 50 changes in astropy were not that large - astropy/astropy#15525 - with some code becoming more and other code becoming less obvious.

Thanks @mdhaber, @mhvk.

We did notice something related to this downstream in SciPy per: scipy/scipy#19526

I didn't open a NumPy issue, because it may actually be an improvement to sensitivity, can't tell just yet..

The failure is

E   Mismatched elements: 1 / 1 (100%)
E   Max absolute difference among violations: nan
E   Max relative difference among violations: nan
E    ACTUAL: array(-8.113438e+295-infj)
E    DESIRED: array(-8.113438e+295-infj)

I think this should still pass; the new failure is not an intended effect of the PR. It is probably failing because the magnitudes are not finite, so exact equality comparison is performed, but presumably the real parts are not exactly equal. What is surprising is that the new logic is similar to the old logic in this respect, so it's not jumping out at me why it didn't fail before. I'll take a closer look soon.

Update:
In the old logic:

        x = x * ones_like(cond)
        y = y * ones_like(cond)

(where x and y are np.asarray([-8.113438e+295-np.inf*1j]) and cond is a boolean array) turns x and y into complex NaNs.

import numpy as np
x = np.asarray([-8.113438e+295-np.inf*1j])
cond = np.asarray([True])
x * np.ones_like(cond)  # array([nan+nanj])

Then the NaN comparison logic kicked in, and the test passed. This was not exactly correct behavior either because the original inputs would be considered close even if the real parts were not close. So I think this case has exposed a bug in both the old and new code.

The obvious workaround to me is to perform real and imaginary part comparison separately, but maybe there is a more elegant solution.

@mhvk would you like to adapt your solution to work with complex numbers in which either the real or complex part is infinite?

I'm actually not sure what the expectation should be: the current test assumes either x, y are finite, or, if not, that they should be the same infinite (same sign). Following that logic, for complex, I guess we can that say that the infinite should point in the same direction (i.e., any of 8 possibilities). Still, for the case in question, of x=finite ± infinite j, I guess one would then logically have to ignore the finite values altogether. Is that what we want?

On the other hand, one can see isclose as a way to look for machine rounding errors, in which case the answer is clearer: just check real and imaginary separately if the other is infinite. But then, we do not do that for reals: the maximum number representable is not treated as the same as infinite, even if it is just a rounding error away. But it is not obvious to me one would ever want to do isclose(x.real, y.real) & isclose(x.imag, y.imag) - it is important that rtol goes with abs(y) which includes both real and imaginary.

Anyway, I think it makes most sense to check for "same infinity" (ignoring real for -8.113438e+295+0j-np.inf*1j), which would make the scipy tests pass again, but probably best if we first agree this is OK...

p.s. It is all a bit weird: if I do,

np.asarray([-8.113438e+295-np.inf*1j])
# array([nan-infj])

This actually is a python problem:

-8.113438e+295-np.inf*1j
(nan-infj)

So, it looks like one has to construct the test values explicitly assigning real and imaginary...

This is what the directional infinity check would look like:

    with errstate(invalid='ignore'), _no_nep50_warning():
        result = less_equal(abs(x-y), atol + rtol * abs(y))
        if not np.all(finite_y := isfinite(y)):
            result &= finite_y
            # We allow matches for infinities pointing in the same direction.
            if x.dtype.kind != "c" and y.dtype.kind != "c":
                result |= (x == y)
            else:
                # For complex, 8 ways to point in the same direction.
                x_real_inf = isinf(x.real)
                x_imag_inf = isinf(x.imag)
                match_real_inf = x_real_inf & (x.real == y.real)
                match_imag_inf = x_imag_inf & (x.imag == y.imag)
                neither_real_inf = ~(x_real_inf | isinf(y.real))
                neither_imag_inf = ~(x_imag_inf | isinf(y.imag))
                result |= match_real_inf & (match_imag_inf | neither_imag_inf)
                result |= match_imag_inf & neither_real_inf

            if equal_nan:
                result |= isnan(x) & isnan(y)