BUG: fix reciprocal of 0+0j to be inf+nanj

[kurtamohler]

This change makes np.reciprocal(0+0j) return the same result as
1 / np.array(0+0j). See #17425

[kurtamohler]

I chose to mimic what np.divide is doing here:

numpy/numpy/core/src/umath/loops.c.src

Lines 2698 to 2702 in f4a4ddd

	if (in2r_abs == 0 && in2i_abs == 0) {
	/* divide by zero should yield a complex inf or nan */
	((@ftype@ *)op1)[0] = in1r/in2r_abs;
	((@ftype@ *)op1)[1] = in1i/in2i_abs;
	}

[eric-wieser]

I don't think this comment is correct - any result is valid as long as it contains at least one infinity

[kurtamohler]

Good point, I've changed the comment.

[seberg]

I wonder if we can compare to some other reciprocal implementation to make sure we are not missing anything. One thing to keep in mind, is that this should give a division warning I think. That could be achieved using npy_set_floatstatus_divbyzero() I think. But I am a bit curious if there is a way to spell it out differently, even if reordering below and testing if 1/d is an infinity.

[mattip]

I thought our policy was to not emit a warning where we can set a value instead.

[seberg]

Sorry, yes... I was being silly it is only an issue when a new NaN is created, but here an infinity is created and the NaN is just a funny side-product.

[seberg]

Hmmm, no wait:

In [6]: np.divide([1.], [0])                                                                                            
<ipython-input-6-7923c6b44821>:1: RuntimeWarning: divide by zero encountered in true_divide
  np.divide([1.], [0])

Further:

In [8]: np.divide([1.], [-0.])                                                                                          
<ipython-input-8-c274e82a14bd>:1: RuntimeWarning: divide by zero encountered in true_divide
  np.divide([1.], [-0.])
Out[8]: array([-inf])

reciprocal of -0. should probably be -inf. That said, this is of course already much better than before.

[kurtamohler]

Should I update the reciprocal of -0+0j and -0-0j to be -inf+nanj? If so, I suppose it would be a good idea to update divide as well, since it returns positive infinity in these cases:

>>> numpy.divide([1.], [-0.-0j])
array([inf+nanj])
>>> numpy.divide([1.], [-0.+0j])
array([inf+nanj])

[seberg]

I am confused :/. I am pretty sure it would be best to set the division flags (but maybe we can defer that). As for the -0., probably that doesn't really matter... Since the magnitude is infinite and the complex part is NaN, I am not certain that -inf would make much sense, since -0-0j, etc. are also possible causes for infinity.

[eric-wieser]

From what I remember, unlike the reals IEEE complex numbers have no semantic differences between infinities, so it doesn't matter which one we return.

[kurtamohler]

Since inf+nanj does not behave nicely in Python, would it be better to set the imaginary value to something other than nan? For instance:

>>> 1.0 / (float('inf') + float('nan')*1j)
(nan+nanj)
>>> 1.0 / (float('inf') + 0j)
0j

[seberg]

Sorry, that this stalled. I am still a bit confused on it, but it seems to me:

1. It doesn't matter what we return for the imaginary part. NaN may not always be ideal, but it is probably the safest bet. For example, if you calculate the np.angle you will get an invalid result and not something that might just be wrong.
2. 1 / 0j raises a ZeroDivisionError in python, so I think we should add a npy_set_floatstatus_divbyzero() here, or write the code in a way that ensures setting the flag.

[eric-wieser]

@kurtamohler, Your test is bad there, float('inf') + float('nan')*1j is not the same as complex(math.inf, math.nan) due to how multiplication works on complex.

[mattip]

ping. This has stalled.

[eric-wieser]

I thought our policy was to not emit a warning where we can set a value instead.

I don't remember seeing that conversation. If that is indeed our policy, can we document it somewhere? Either in the ufunc docs, or perhaps in np.errstate

[seberg]

I think the floating point error flags for division by zero, overflow, and "invalid" are normally always set when a new NaN or inf is created and never when one is propagated. (Even inf / inf sets it in NumPy, although not in Python since a new NaN is created out of inf.)

[seberg]

The PR seems pretty far along, but does not yet make sure that the correct warnings are given. The complex division code does this correctly (both divide by zero and invalid value).

This also needs tests (which could also just test against np.divide).

[kurtamohler]

Sorry, I forgot all about this. I will fix it up in the next couple days

[kurtamohler]

I've added the warning and updated the test to check it

[seberg]

It would be cool to extend the test to all complex dtypes, but the whole file doesn't do this, so I can ignore that part.

However, np.divide(1, 0j) gives both division by zero (read "a new inf is created") and invalid value (read "a new NaN is created"). The invalid value warning is missing and should also be tested.

But, besides that small addition, I think this is ready. I am not sure if there is a preference for a specific infinity, but I could not find one so I think it can wait until someone else comes around and wants it.

[InessaPawson]

@kurtamohler Let us know if you have any questions regarding extending the test.