Should numpy.sqrt(-1) return 1j rather than nan?

Christopher Barker Chris.Barker at noaa.gov
Thu Oct 12 12:49:25 CDT 2006


As someone on this list (sorry to quick with the delete button) said:

"numpy is closer to the metal than MATLAB"

MATLAB aside, numpy is somewhat "close to the metal". I think this is 
clearly part of its design philosophy, and also something I personally 
like about out.

Getting back to MATLAB (at least version 4, the last one I used 
regularly) -- everything in MATLAB is a matrix of double precision 
Floating points. (or complex - that does get weird!). This makes things 
simpler for a lot of use, but restricts the power and flexibility 
substantially.

Everything in numpy is an array of arbitrary shape and many possible 
data types -- this is more to think about, but give you a lot of power 
and flexibility. In any case, that is what numpy is. Period.

Given that, arrays of different types behave differently -- that's 
something you'd darn well better understand pretty soon in your numpy 
career:

 >>> a = N.array((1,2,3), dtype=N.uint8)
 >>> a *= 200
 >>> a
array([200, 144,  88], dtype=uint8)

Oh my gosh! 2 * 200 isn't 144! what happened?

This isn't any different than sqrt(-1) resulting in NaN or an error. In 
fact, there are all sorts of places where upcasting is not happening 
automagically -- in fact, I think that's become more consistent in the 
new numpy. So, in numpy you need to choose the datatype appropriate for 
your problem at hand. I know I always know whether I want to work in the 
complex plane or not.

Another change recently made is to make doubles the default in more 
places -- I like that change.

So, given the entire philosophy and current status of how numpy works, 
the ONLY question that is legitimate here is:

Should the "default" datatype be complex or double? I vote double.

Travis Oliphant wrote:
> Now, it would be possible to give ufuncs a dtype keyword argument that 
> allowed you to specify which underlying loop was to be used for the 
> calculation. 

> sqrt(a, dtype=complex). 

I like this OK, and isn't is similar to the dtype argument in the 
accumulating methods? (sum, etc.). However, it's probably not worth the 
C-api change -- is it that different than calling:

sqrt(a.astype(complex)) ?

As for scipy, I heartily agree that scipy should be a superset of numpy, 
and NOT change the behavior of numpy functions and methods. We're taking 
the opportunity at this point for numpy to break backward compatibility 
with Numeric/numarray -- this is probably the best time to   do the same 
for scipy.

couldn't we at least introduce new names?

numpy.scimath.csqrt()

for instance?

-Chris


-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chris.Barker at noaa.gov

-------------------------------------------------------------------------
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642




More information about the Numpy-discussion mailing list