Repository navigation
BUG: Invalid corner cases (resulting in nan+nanj) in complex division #26560
Description
Describe the issue:
It seems, that numpy uses Smith's algorithm to do complex true division. In some cases (see below) it produces (nan+nanj), but meaningful components could be recovered. See e.g. C11 Annex G.5.2, _Cdivd()'s example.
FYI: see also CPython's issue python/cpython#119372
Reproduce the code example:
>>> 1/np.complex128('inf')
0j
>>> 1/np.complex128('infj')
-0j
>>> 1/np.complex128('(inf+infj)') # oops, should be 0j
<stdin>:1: RuntimeWarning: invalid value encountered in scalar divide
(nan+nanj)Error message:
No response
Python and NumPy Versions:
>>> np.__version__
'1.26.4'Runtime Environment:
No response
Context for the issue:
No response
Activity
mdickinson commentedon May 29, 2024
For a bit of context the PR that @skirpichev mentions, python/cpython#119372, is a currently open PR against core Python that (if merged) will change the semantics of complex division for Python
complexobjects in these corner cases. NumPy (and in particular the NumPy test suite, if it tests for compatibility betweencomplexandnp.complex128) seems like the most likely 3rd party package to be broken by such a change, and I wanted to make sure that the NumPy developers were aware of the pending change and had an opportunity to push back if they deem the risk of breakage to be too high.seberg commentedon May 29, 2024
Some libraries may test against NumPy and copy us bug-for-bug, but those should just keep up.
NumPy tries to follow Python in many things, not sure I would say necessarily strictly for this kind of edge case (unlike
divmod). Seems like this one is pretty clear cut an improvement if the C standard defines it (and I doubt performance matters either).I believe we have actual a similar or even identical PR to fix
reciprocalfor this corner case.In any case, long story short: I think you can change this and NumPy should also fix it in either case. Neither does it sound like a burning issue or an issue if NumPy "diverges".
mdickinson commentedon May 29, 2024
Thanks, @seberg!
skirpichev commentedon May 30, 2024
JFR:
numpy/numpy/_core/src/npymath/npy_math_complex.c.src
Lines 105 to 107 in 012e63d
0+0jto beinf+nanj #17449 It's outdated, but I doubt it's identical or even a similar. In fact, mentioned CPython PR doesn't touch that condition of the Smith's algorithm (python just raises a ZeroDivisionError here).I'll work on a patch.
ikrommyd commentedon Sep 24, 2026
The reason for this problem is that NumPy does not perform the NaN recoveries for complex multiplication and division as defined in the Annex G of the C99 (and newer) standards. I had opened a PR that fixed this (#30806) but we have decided to not do this now. It adds a noticeable (few %) overhead in the nominal case (no NaNs to correct) and up to 50% overhead in the cases where there are a decent amount NaNs to correct. An additional argument is that no other array libraries do this and cupy and jax for example do not do these Annex G recoveries.
Edit: jax actually does implement Annex G recoveries but only for division.