test_square test failure with numpy complex64

[asmeurer]

E                   AssertionError: out=(4+512128000j), but should be roughly x²=512128008j [square()]
E                     x=(16002+16002j)
E                   assert False
E                    +  where False = isclose(4.0, 0.0, 3.4028234663852886e+38)
E                    +    where 4.0 = (4+512128000j).real
E                    +    and   0.0 = 512128008j.real
E                   Falsifying example: test_square(
E                       x=Array(16002.+16002.j, dtype=complex64),
E                   )

I haven't investigated further other than to confirm that it does indeed do that

>>> np.square(np.asarray(16002.+16002.j, dtype=np.complex64))
(4+512128000j)

I don't know if this is a bug in NumPy or just an expected numerical inaccuracy that we shouldn't b testing. I suspect the latter. Note that the relative error is around 1e-8

>>> abs(512128008j - (4+512128000j))/abs(512128008j)
1.7464914572684647e-08

OTOH, PyTorch does give a purely imaginary answer as expected:

>>> torch.square(torch.asarray(16002.+16002.j, dtype=torch.complex64))
tensor(0.+5.1213e+08j)

[asmeurer]

Actually this failure comes up with all array libraries. It seems the test is just wrong for floating-point arithmetic.

[asmeurer]

A related error with numpy

array_api_tests/test_operators_and_elementwise_functions.py:1403: in test_square
    unary_assert_against_refimpl(
array_api_tests/test_operators_and_elementwise_functions.py:264: in unary_assert_against_refimpl
    assert isclose(scalar_o.real, expected.real, M), msg
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

a = -inf, b = nan, M = 1.7976931348623157e+308

    def isclose(
        a: float,
        b: float,
        M: float,
        *,
        rel_tol: float = 0.25,
        abs_tol: float = 1,
    ) -> bool:
        """Wraps math.isclose with very generous defaults.
    
        This is useful for many floating-point operations where the spec does not
        make accuracy requirements.
        """
        if math.isnan(a) or math.isnan(b):
>           raise ValueError(f"{a=} and {b=}, but input must be non-NaN")
E           ValueError: a=-inf and b=nan, but input must be non-NaN
E           Falsifying example: test_square(
E               x=array(1.+8.98846567e+307j),
E           )
E           
E           You can reproduce this example by temporarily adding @reproduce_failure('6.78.3', b'AXicY2ZkQAAQe/8HKAcADsUBtQ==') as a decorator on your test case

>>> import numpy as np
>>> np.square(np.array(1.+8.98846567e+307j))
<stdin>:1: RuntimeWarning: overflow encountered in square
(-inf+1.797693134e+308j)

[asmeurer]

The example in the OP varies between NumPy versions:

>>> import numpy as np
>>> np.square(np.asarray(16002.+16002.j, dtype=np.complex64))
(4+512128000j)
>>> np.__version__
'1.25.0.dev0+1503.g58224647e8'

>>> import numpy as np
>>> np.square(np.asarray(16002.+16002.j, dtype=np.complex64))
512128000j
>>> np.__version__
'1.24.2'

The spec for square doesn't say anything about how accurate the result should be, other than that the special cases should match for x*x (I don't think that applies to any of the examples here).

A big part of the problem is that the helpers are comparing the real and imaginary parts separately. For example, here

array-api-tests/array_api_tests/test_operators_and_elementwise_functions.py

Lines 264 to 265 in f82c7bc

	assert isclose(scalar_o.real, expected.real, M), msg
	assert isclose(scalar_o.imag, expected.imag, M), msg

For floating-point closeness, a complex number should be compared directly, for example, by using cmath.isclose, or by applying the isclose function directly to the complex number (like abs(x - y) or whatever). The relative tolerance of a complex floating-point computation applies to the whole complex number, meaning that if the imaginary part is much larger than the real part, as in this example, the answer is still relatively close.

Also, for square, I'm not sure if we should really test against x**2 at all. If we are going to do that, x*x would be better (although for NumPy they seem to be the same), and with very low tolerances.