Numpy - square root of -1 leaves a small real part

Question

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.

Answer 1

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).