14

Perhaps this is an algorithmic issue, but the following piece of code

numpy.power((-1+0j),0.5)

produces the following output

(6.1230317691118863e-17+1j)

Analogous expressions e.g. numpy.power(complex(-1),.5) yield the same result, however - numpy.sqrt(complex(-1)) yields the expected result of 1j. Clearly the result should have no real part, so I am missing something crucial or do I need to report this to numpy dev's.

In case anyone asks, no I can not round away the real part (I need full precision for this calculation) and yes I need to use the power function.

Eric O. Lebigot
  • 91,433
  • 48
  • 218
  • 260
crasic
  • 1,918
  • 2
  • 18
  • 29
  • 4
    What do you mean "I need full precision for this calculation"? The value you have given is smaller than machine epsilon on my computer: `np.finfo(np.complex).eps == 2.2204460492503131e-16`. – Katriel Jun 09 '11 at 10:00
  • @katrielalex What I mean is that when there shouldn't be any real part, there shouldn't be any real part. Perhaps precision isn't the right word, but for example, this result is later plugged into the complex erf and the results later down the line are sensitive to the existence or not of a real component (this script is for calculations in quantum physics). I know that the precision of say, root(2) is limited by a larger number than this, but in this case, the results are rather exact if that makes any sense. – crasic Jun 09 '11 at 20:32
  • 1
    Precision of most complex floating point functions (e.g. in the standard C library, where `power = cpow` in this case comes from) is often only specified in terms of |error|/|value| metric. In that metric, the above result is as good as it gets in floating point. If you require more, you'll have to put in the extra knowledge by hand. – pv. Jun 10 '11 at 12:20

3 Answers3

15

What happens is that the square root of -1 is calculated as exp(i phase/2), where the phase (of -1) is approximately π. In fact,

>>> import cmath, math
>>> z = -1+0j
>>> cmath.phase(z)
3.141592653589793
>>> math.cos(_/2)
6.123233995736766e-17

This shows that the phase of -1 is π only up to a few 1e-17; the phase divided by 2 is also only approximately π/2, and its cosine is only approximately 0, hence your result (the real part of your result is this cosine).

The problem comes ultimately from the fact that there is only a fixed, finite number of floating point numbers. The number π is not in the list of floating point numbers, and can therefore only be represented approximately. π/2 cannot be exactly represented either, so that the real part of the square root of -1 is the cosine of the floating point approximation of π/2 (hence a cosine that differs from 0).

So, Python's approximate value for numpy.power(complex(-1), .5) is ultimately due to a limitation of floating point numbers, and is likely to be found in many languages.

What you observe is connected to this floating point limitation, through the implementation of the power of a number. In your example, the square root is calculated by evaluating the the module and the argument of your complex number (essentially via the log function, which returns log(module) + i phase). On the other hand, cmath.sqrt(-1) gives exactly 1j because it uses a different method, and does not suffer from the floating point approximation problem of (-1+0j)**0.5 (as suggested by TonyK).

Eric O. Lebigot
  • 91,433
  • 48
  • 218
  • 260
9

It's a side effect of the implementation of numpy.power() for complex numbers. The stdlib exhibits the same issue.

>>> numpy.power(-1+0j, 0.5)
(6.123233995736766e-17+1j)
>>> cmath.exp(cmath.log(-1)/2)
(6.123233995736766e-17+1j)
Ignacio Vazquez-Abrams
  • 776,304
  • 153
  • 1,341
  • 1,358
  • So is sqrt a separate function? I somehow always assumed it called power – crasic Jun 09 '11 at 08:53
  • 14
    In general, the complex `power` function has to evaluate (real) logarithms and trigonometric functions. But the `sqrt` function only needs to solve a couple of quadratic equations. So a good implementation will provide a custom `sqrt` function, which is faster and more accurate than `power`. – TonyK Jun 09 '11 at 08:59
5

Try numpy.real_if_close().

See: Real if close

roger314
  • 51
  • 1