13

I cannot understand the following output. I would expect Numpy to return -10 (or an approximation). Why is it a complex number?

print((-1000)**(1/3.))

Numpy answer

(5+8.660254037844384j)

Numpy official tutorial says the answer is nan. You can find it in the middle of this tutorial.

Glorfindel
  • 21,988
  • 13
  • 81
  • 109
user1700890
  • 7,144
  • 18
  • 87
  • 183
  • 7
    Why do you think `numpy` is involved in this process? Also, have a look at `(5+8.660254037844384j)**3`! – jonrsharpe Jul 05 '15 at 13:59
  • 3
    For the answer to the mathematical half of this question, see e.g. http://math.stackexchange.com/questions/25528/cubic-root-of-negative-numbers. There is more than one cubic root of `-1000`! – jonrsharpe Jul 05 '15 at 14:12
  • Thank you jonrsharpe. I completely forgot about roots of unity. How can I force python to return -10? – user1700890 Jul 05 '15 at 14:33
  • I'm not sure you can. – jonrsharpe Jul 05 '15 at 14:34
  • 2
    @user1700890: You could write your own `cbrt` function, using something like: `def cbrt(x): return copysign(abs(x)**(1/3.), x)`. If you import `copysign` from `numpy` instead of `math`, this definition should work for arrays as well as floats. – Mark Dickinson Jul 05 '15 at 14:36

1 Answers1

7

You are exponentiating a regular Python scalar rather than a numpy array.

Try this:

import numpy as np

print(np.array(-1000) ** (1. / 3))
# nan

The difference is that numpy does not automatically promote the result to a complex type, whereas a Python 3 scalar gets promoted to a complex value (in Python 2.7 you would just get a ValueError).

As explained in the link @jonrsharpe gave above, negative numbers have multiple cube roots. To get the root you are looking for, you could do something like this:

x = -1000
print(np.copysign(np.abs(x) ** (1. / 3), x))
# -10.0

Update 1

Mark Dickinson is absolutely right about the underlying cause of the problem - 1. / 3 is not exactly the same as a third because of rounding error, so x ** (1. / 3) is not quite the same thing as the cube root of x.

A better solution would be to use scipy.special.cbrt, which computes the 'exact' cube root rather than x ** (1./3):

from scipy.special import cbrt

print(cbrt(-1000))
# -10.0

Update 2

It's also worth noting that versions of numpy >= 0.10.0 will have a new np.cbrt function based on the C99 cbrt function.

ali_m
  • 71,714
  • 23
  • 223
  • 298
  • Thank you, ali_m! Why would not Numpy return -10? – user1700890 Jul 05 '15 at 14:35
  • *All* real (non-zero) numbers have a pair of complex cube roots in addition to a real cube root, not just the negative ones. – chepner Jul 05 '15 at 23:14
  • 2
    @user1700890: Note that this isn't computing the cube root: it's computing `-1000` to the power `0.333333333333333314829616256247390992939472198486328125`, which isn't the same thing at all. Having the power operation guess that you actually wanted a cube root and defer to `cbrt` in this situation would break continuity and in general be a horrible thing to do. (Should it also guess for 5th roots? 7th roots? Using what criterion?) IEEE 754 *does* define a `rootn` function, which behaves the way you'd expect for odd `n`, giving finite results for negative finite inputs. – Mark Dickinson Jul 06 '15 at 07:28
  • @ali_m: Nice find with `scipy.special.cbrt`. It would be nice to see Python's math module grow either a `cbrt` or a `rootn` function for this situation. – Mark Dickinson Jul 06 '15 at 17:42
  • @MarkDickinson Yes, I was hoping there would be a generic `rootn` function somewhere in the Python standard libraries (or at least in numpy). The [scipy source for `cbrt`](https://github.com/scipy/scipy/blob/d1203c280a8661536d64f17e6a90d7ef000302d5/scipy/special/cephes/cbrt.c) ultimately comes from cephes, which unfortunately lacks a `rootn`. – ali_m Jul 06 '15 at 18:59