Every now and then, the behaviour of `sqrt`` for negative values is discussed on the mailing list. Hopefully, by summarising the latest thread] here, we don't need to do so again in the future. `

## Why does numpy.sqrt(-1) return nan? Why does it not yield 1j, which is surely the correct answer?

If you assume the domain of computation is the field of complex numbers, then yes, the above assumption is true -- that the square root of `-1`` is ``1j``. However, ``numpy`` supports many different data types, and in the context of these data types, the answer may vary. `

By default, `numpy`` assumes that you are working with real numbers (which is a reasonable assumption). What you are then asking ``numpy`` is: "What is the square root of ``-1`` (but please limit your answer to real numbers)?". Since the answer is ``1j``, numpy cannot provide a real value, so it rather yields ``nan``. `

On the other hand, if you specifically ask `numpy`` for the complex answer in one of several ways, it will be happy to oblige: `

```
In [49]: N.sqrt(-1+0j)
Out[49]: 1j
In [50]: N.sqrt(complex(-1))
Out[50]: 1j
In [50]: N.sqrt(N.asarray(-1,dtype=complex))
Out[50]: 1j
```

Simply put, `sqrt``'s output is of the same type as the input (``numpy`` does no checking). The only exception is for integers, in which case ``sqrt`` returns floats. `

## Why can't numpy just provide complex numbers by default?

Imagine having to calculate the square roots of a large array of numbers. If numpy calculates complex answers, that would mean that the complex answer would take up twice as much memory as the real answer. Clearly, this case is not a desirable default. Of course, `numpy`` could check whether negative numbers are present, and only *then* calculate complex answers, but again, this leads to inconsistent behaviour. Furthermore, if a user wants to operate on positive real numbers only, why burden his calculations with the extra time-consuming complexity? `

## How can I change the default behaviour?

`Numpy`` provides functions to automatically return complex numbers for negative real arguments. These reside in ``numpy.lib.scimath`` and can be imported using `

```
from numpy.lib.scimath import *
```

to replace `sqrt`` in the current namespace. `

However, it is better style to use

```
from numpy.lib import scimath as SM
SM.sqrt(-1)
```

Note that this doesn't actually change the default behavior. Still `numpy.sqrt(-1)`` will return ``nan``. It just makes the version of sqrt that checks for complex results more accessible via ``SM.sqrt`` instead of ``numpy.lib.scimath.sqrt``. `

Another alternative is to take advantage of the fact that any python file is a module, and make a mini-module that makes complex behavior the default. Then you just import that instead of numpy.

```
# File: cnumpy.py
from numpy import *
from numpy.lib.scimath import *
```

Then from your files, or from the interpreter, you use your `cnumpy`` module like so: `

```
>>> import cnumpy
>>> cnumpy.sqrt(-1)
1j
>>> a = cnumpy.ones((2,2)) * -1
>>> cnumpy.sqrt(a)
array([[ 0. +1.00000000e+00j, 0. +1.00000000e+00j],
[ 0. +1.00000000e+00j, 0. +1.00000000e+00j]])
```

In all respects it looks just like a new version of `numpy`` that happens to return ``1j`` for ``sqrt(-1)``. `

#### Changing error reporting behavior

The `numpy`` error handler can also be asked to complain if negative values are passed to ``sqrt``: `

```
In [1]: N.sqrt(-1)
Out[1]: nan
In [2]: saved_error_handler = N.seterr(invalid='raise') # or 'warn'
In [3]: N.sqrt(-1)
---------------------------------------------------------------------------
exceptions.FloatingPointError Traceback (most recent call last)
/home/stefan/work/prasa2006/<ipython console>
FloatingPointError: invalid encountered in sqrt
In [4]: N.seterr(**saved_error_handler)
Out[4]: {'over': 'ignore', 'divide': 'ignore', 'invalid': 'warn', 'under': 'ignore'}
In [5]: N.sqrt(-1)
Out[5]: nan
```

With Python 2.5 it should soon be possible to do

```
with errstate(invalid='raise'):
sqrt(-1) # raise exception
```

Alternatively, you can simply use complex arrays throughout.

## Do numpy and scipy behave differently?

Yes, `scipy`` exposes ``numpy.lib.scimath`` by default. `

```
>>> import scipy
>>> scipy.sqrt(-1)
1j
```