Invalid corner cases (resulting in nan+nanj) in complex division about cpython

[eendebakpt]

I can confirm the behavior described above. Some similar cases:

In [24]: 1/(cmath.inf + cmath.infj)
Out[24]: (nan+nanj)

In [25]: cmath.inf * 1j
Out[25]: (nan+infj)

But numpy shows the same behaviour as current cpython main

In [47]: import cmath
    ...: import numpy as np
    ...: print(np.__version__)
    ...: x = np.array([np.inf], dtype=np.complex128)
    ...: y = np.array([cmath.infj])
    ...: print(1. / (x + y))
1.26.4
[nan+nanj]

C:\Users\eendebakpt\AppData\Local\Temp\ipykernel_10856\2212785139.py:6: RuntimeWarning: invalid value encountered in divide
  print(1. / (x + y))

[serhiy-storchaka]

I agree. If the numerator is finite and the denominator is infinite, the result should be zero. If the numerator is also infinite or contains a NaN, the result should contain a NaN (but not always both component should be NaN). Do you mind to create a PR? Pay attention to the sign of zeros, in some cases it can be clearly defined.

[skirpichev]

Some similar cases

Please note that this issue was restricted to true division and (nan+nanj) output.

1/(cmath.inf + cmath.infj)

That's similar to my example: CPython's complex arithmetic has no mixed-mode rules (e.g. for float x complex), so it's an equivalent of (1+0j)/(cmath.inf + cmath.infj). Should be 0j.

cmath.inf * 1j

Seems to be valid: inf * 1j == (inf+0j)*(0+1j)=(inf*0-0*1)+(inf*1+0*0)j=(nan+infj). You might expect infj output, per C11 Annex G, but this will imply mentioned above mixed-mode rules for complex arithmetic (and special imaginary data type, e.g. for 1j).

But we have some cases, when CPython math give (nan+nanj) for multiplication and C code - something different. If this worth an issue - it should be a separate one.

an example

>>> complex('(inf+1j)')*complex('(nan+1j)')
(nan+nanj)

#include <complex.h>
#include <stdio.h>
#include <math.h>

int
main()
{
    complex double z;
    z = CMPLX(INFINITY, 1) * CMPLX(NAN, 1);
    printf("%f%+fi\n", creal(z), cimag(z));
    return 0;
}

$ cc -std=c11 a.c -lm && ./a.out
-nan+infi

But numpy shows the same behaviour as current cpython main

I guess, implementation copied from CPython.

Do you mind to create a PR?

@serhiy-storchaka, yes, see above diff.

[skirpichev]

FYI, PR #119457 is ready for review.