BUG: invalid value encountered in power resulting in NaNs

[jecampagne]

Describe the issue:

NaN values in exponentiation of y as np.array either y**val or np.power(y,val) depending on val values but strange. I'm using fresh version from conda install 1.21.2. May be someone can tell me if he/she reproduces the problem

Reproduce the code example:

import numpy as np
y = np.array([-8.18254833e-09, -9.87463911e-09, -1.19166418e-08, -1.43809156e-08,
 -1.73547831e-08, -2.09436246e-08, -2.52746120e-08, -3.05012158e-08,
 -3.68086427e-08, -4.44203992e-08, -5.36062108e-08, -6.46915807e-08,
 -7.80693235e-08, -9.42134850e-08, -1.13696140e-07, -1.37207665e-07,
 -1.65581201e-07, -1.99822175e-07, -2.41143930e-07, -2.91010720e-07,
 -3.51189595e-07, -4.23813019e-07, -5.11454432e-07, -6.17219445e-07,
 -7.44855884e-07, -8.98886599e-07, -1.08476973e-06, -1.30909213e-06,
 -1.57980275e-06, -1.90649433e-06, -2.30074333e-06, -2.77652011e-06,
 -3.35068404e-06, -4.04358084e-06, -4.87976360e-06, -5.88886280e-06,
 -7.10663629e-06, -8.57623638e-06, -1.03497389e-05, -1.24899888e-05,
 -1.50728266e-05, -1.81897762e-05, -2.19512880e-05, -2.64906528e-05,
 -3.19687249e-05, -3.85796219e-05, -4.65576038e-05, -5.61853736e-05,
 -6.78040954e-05, -8.18254833e-05])

y**4.0  # np.power(y,4.0)  gives correct values
y**3.0000000000000001   # gives correct values
y**3.000000000000001 # leeds to Nan

array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])

Error message:

RuntimeWarning: invalid value encountered in power
  y**3.000000000000001

NumPy/Python version information:

1.21.2 3.8.12 (default, Oct 12 2021, 13:49:34)
[GCC 7.5.0]

[PeterTillema]

I can reproduce it, but the problem is that 3.0000000000000001 is too small to be represented by NumPy types, so it automatically gets rounded to 3.0. 3.000000000000001 is large enough so the actual value is indeed 3.000000000000001, which is an invalid power value.

See this:

>>> a = np.array([3.0000000000000001])
>>> a.tobytes().hex()
'0000000000000840'
>>> a = np.array([3.000000000000001])
>>> a.tobytes().hex()
'0200000000000840'

[jecampagne]

Hello @PeterTillema In fact you get also pbs with 3.1 as exponent. But I realize that all the y-values are negatives, so I was expected either complex values or a prescription like sign(y) |y|^val.

[WarrenWeckesser]

Cross-reference (so it can also be closed when we close this issue): #12118

Edit: The issue of getting nan when negative values are raised to a non-integer float comes up now and then. Here are a few similar issues:

#14913
#14273
#12319
#12182
#12118

Edit 2: Added another issue (still open at the moment) to the above list.

[melissawm]

Is this a bug, or could adding a note to the effect of #14273 (comment) to the docs help?

[seberg]

There is no bug. We should not change the default of returning float when the input is float (and should keep scalars aligned to arrays). There are the, maybe somewhat controversial, np.emath.… functions that try to get around the issue, they do this by inspecting if the data/result(?) is "almost real" though, I think.

But in the end, you should probably make sure the input is complex if you want a complex output, or use dtype=complex.

[WarrenWeckesser]

A comment about this behavior in the docstrings of np.power and np.float_power is a good idea. Updating the docs won't necessarily cut down on reports like this in the future (since I suspect many new users confused by the difference between (-1.5)**2.5 and np.array(-1.5)**2.5 won't think to look in the docstring of np.power) but it will allow us to respond quickly with a pointer to the docs, instead of explaining it again in the new issue or looking through the old issues to provide a link to a good explanation.