Robert Kern wrote:

Stefan van der Walt wrote:

...

Hi all

Two quick questions regarding unintuitive numpy behaviour:

Why is the square root of -1 not equal to the square root of -1+0j?

In [5]: N.sqrt(-1.) Out[5]: nan

In [6]: N.sqrt(-1.+0j) Out[6]: 1j

It is frequently the case that the argument being passed to sqrt() is expected to be non-negative and all of their code strictly deals with numbers in the real domain. If the argument happens to be negative, then it is a sign of a bug earlier in the code or a floating point instability. Returning nan gives the programmer the opportunity for sqrt() to complain loudly and expose bugs instead of silently upcasting to a complex type. Programmers who *do* want to work in the complex domain can easily perform the cast explicitly.

...

Is there an easier way of dividing two scalars than using divide?

In [9]: N.divide(1.,0) Out[9]: inf

x/y ?

...

(also

In [8]: N.divide(1,0) Out[8]: 0

should probably ruturn inf / nan?)

inf and nan are floating point values. The definition of int division used when both arguments to divide() are ints also yields ints, not float

This relates to the discussion that Travis and I we're having about error handling last week. The current defaults for handling errors is to ignore them all. This is for speed reasons, although our discussion may have alleviated some of these. The numarray default was to ignore underflow, but warn for the rest; this seemed to work well in practice. However, this example points in another possible direction.... Travis mentioned that checking the various error conditions in integer operations was painful and slowed things down since there wasn't machine support for it. My current opinion is that we should just punt on overflow and let integers overflow silently. That's what bit twiddlers want anyway and it'll be somewhere between difficult and impossible to do a good job. I don't think invalid and underflow apply to integers, so that leaves divide. I think me preference here would be for int divide to raise by default. That would require that there by five error classes, shown here with my preferred defaults: divide_by_zero="warn", overflow="warn", underflow="ignore", invalid="warn" int_divide_by_zero="raise" The first four apply to floating point (and complex) operations, while the last applies to integer operations. The separation of warnings into two classes also helps avoid the expectation that we should be doing something useful about integer overflow. I don't *think* this should be too difficult; just stick a int_divide_by_zero flag on some thread_local variable and set it to true when there's been a divide by zero, checking on the way out of the ufunc machinery. I haven't tried it though, so it may be much harder than I envision. In any event , the current divide by zero checking seems to be a bit broken. I took a quick look at the code and it's not obvious why, (unless my optimizer is eliding the error generation code?). This is the behaviour I see under windows compiled using VC7:

...

...

one = np.array(1) zero = np.array(0) one/zero 0 np.seterr(divide='raise') one/zero # Should raise an error 0 (one*1.0 / zero) # Works for floats though?! Traceback (most recent call last): File "<stdin>", line 1, in ? FloatingPointError: divide by zero encountered in divide

Regards, -tim

Stefan van der Walt wrote:

On Wed, Apr 12, 2006 at 01:14:54AM -0500, Robert Kern wrote:

...

Stefan van der Walt wrote:

...

Why is the square root of -1 not equal to the square root of -1+0j?

In [5]: N.sqrt(-1.) Out[5]: nan

In [6]: N.sqrt(-1.+0j) Out[6]: 1j

It is frequently the case that the argument being passed to sqrt() is expected to be non-negative and all of their code strictly deals with numbers in the real domain. If the argument happens to be negative, then it is a sign of a bug earlier in the code or a floating point instability. Returning nan gives the programmer the opportunity for sqrt() to complain loudly and expose bugs instead of silently upcasting to a complex type. Programmers who *do* want to work in the complex domain can easily perform the cast explicitly.

The current docstring (specified in generate_umath.py) states

y = sqrt(x) square-root elementwise.

It would help a lot if it could explain the above constraint, e.g.

y = sqrt(x) square-root elementwise. If x is real (and not complex), the domain is restricted to x>0.

I'll get around to it sometime. In the meantime, please make a ticket: http://projects.scipy.org/scipy/numpy/newticket

...

...

In [9]: N.divide(1.,0) Out[9]: inf

x/y ?

On my system, x/y (for x=0., y=1) throws a ZeroDivisionError. Are the two divisions supposed to behave the same?

Not exactly, no. Specifically, the error handling is, by design, more flexible with numpy than regular float objects. If you want that flexibility, then you need to use numpy scalars or ufuncs. -- Robert Kern robert.kern@gmail.com "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