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

[skirpichev]

Bug report

Bug description:

Reproducer:

>>> (1+1j)/(cmath.inf+0j)  # ok
0j
>>> (1+1j)/cmath.infj  # ok
-0j
>>> (1+1j)/(cmath.inf+cmath.infj)  # should be 0j, isn't?
(nan+nanj)

c.f. MPC:

>>> gmpy2.mpc('(1+1j)')/gmpy2.mpc('(inf+infj)')
mpc('0.0+0.0j')

It seems, there are not so much such cases (i.e. when the numerator is finite and the denominator has infinite and infinite/nan components, or if the numerator has infinite components and the denominator has zeros). Following patch (mostly literally taken from the C11 Annex G.5.1 _Cdivd()'s example, except for the denom == 0.0 case) should fix the issue:

diff --git a/Objects/complexobject.c b/Objects/complexobject.c
index b3d57b8337..9041cbb8ae 100644
--- a/Objects/complexobject.c
+++ b/Objects/complexobject.c
@@ -122,6 +122,27 @@ _Py_c_quot(Py_complex a, Py_complex b)
         /* At least one of b.real or b.imag is a NaN */
         r.real = r.imag = Py_NAN;
     }
+
+    /* Recover infinities and zeros that computed as nan+nanj */
+    if (isnan(r.real) && isnan(r.imag)) {
+        if ((isinf(a.real) || isinf(a.imag))
+            && isfinite(b.real) && isfinite(b.imag))
+        {
+            const double x = copysign(isinf(a.real) ? 1.0 : 0.0, a.real);
+            const double y = copysign(isinf(a.imag) ? 1.0 : 0.0, a.imag);
+            r.real = Py_INFINITY * (x*b.real + y*b.imag);
+            r.imag = Py_INFINITY * (y*b.real - x*b.imag);
+        }
+        else if ((fmax(abs_breal, abs_bimag) == Py_INFINITY)
+                 && isfinite(a.real) && isfinite(a.imag))
+        {
+            const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);
+            const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);
+            r.real = 0.0 * (a.real*x + a.imag*y);
+            r.imag = 0.0 * (a.imag*x - a.real*y);
+        }
+    }
+
     return r;
 }
 

Perhaps, we could also clear ending remark in _Py_c_quot() comment:

cpython/Objects/complexobject.c

Lines 91 to 92 in d065edf

	* University). As usual, though, we're still ignoring all IEEE
	* endcases.

I'll prepare a PR, if this issue will be accepted.

CPython versions tested on:

CPython main branch

Operating systems tested on:

No response

Linked PRs

• gh-119372: recover inf's and zeros in _Py_c_quot #119457

[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.