[Numpy-discussion] Calculating roots with negative numbers

[Robert Kern]

On Tue, Jul 15, 2008 at 21:59, Matthias Hillenbrand
<mailanhilli at googlemail.com> wrote:
> Hello,
>
> I want to calculate the root of a numpy array with negative numbers.
> Here is an example:
>
> x = linspace(-10,10,100)
> h = zeros(100)
>
> h[where(abs(x) < 2)]  = sqrt( abs(x) -2 )
> h[where(2 <= abs(x))]  = 1j * sqrt( 2 - abs(x) )
>
> Unfortunately I get the following error:     Warning: invalid value
> encountered in sqrt
> and h contains some NaN-values

This scheme would work if you masked the right-hand-sides, too.

In [14]: from numpy import *
In [15]: x = linspace(-4, 4, 21)
In [16]: h = empty(21, dtype=complex)
In [17]: m = abs(x) > 2
In [18]: h[m] = sqrt(abs(x[m]) - 2)
In [19]: h[~m] = 1j*sqrt(2 - abs(x[~m]))
In [20]: h
Out[20]:
array([ 1.41421356+0.j        ,  1.26491106+0.j        ,
        1.09544512+0.j        ,  0.89442719+0.j        ,
        0.63245553+0.j        ,  0.00000000+0.j        ,
        0.00000000+0.63245553j,  0.00000000+0.89442719j,
        0.00000000+1.09544512j,  0.00000000+1.26491106j,
        0.00000000+1.41421356j,  0.00000000+1.26491106j,
        0.00000000+1.09544512j,  0.00000000+0.89442719j,
        0.00000000+0.63245553j,  0.00000000+0.j        ,
        0.63245553+0.j        ,  0.89442719+0.j        ,
        1.09544512+0.j        ,  1.26491106+0.j        ,
1.41421356+0.j        ])


> How can I make a correct case differentiation in combination with
> roots that might contain negative values?

The versions of functions in scimath will extend the domain to include
negatives and return the complex result. You don't need to do any
testing yourself.

In [21]: from numpy.lib import scimath

In [22]: scimath.sqrt(abs(x) - 2)
Out[22]:
array([ 1.41421356+0.j        ,  1.26491106+0.j        ,
        1.09544512+0.j        ,  0.89442719+0.j        ,
        0.63245553+0.j        ,  0.00000000+0.j        ,
        0.00000000+0.63245553j,  0.00000000+0.89442719j,
        0.00000000+1.09544512j,  0.00000000+1.26491106j,
        0.00000000+1.41421356j,  0.00000000+1.26491106j,
        0.00000000+1.09544512j,  0.00000000+0.89442719j,
        0.00000000+0.63245553j,  0.00000000+0.j        ,
        0.63245553+0.j        ,  0.89442719+0.j        ,
        1.09544512+0.j        ,  1.26491106+0.j        ,
1.41421356+0.j        ])

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco