[Numpy-discussion] Re: indexing problem

[Tim Hochberg]

David M. Cooke wrote:

>Tim Hochberg <tim.hochberg at cox.net> writes:
>
>  
>
>>Ryan Krauss wrote:
>>
>>    
>>
>>>At the risk of sounding silly, can you explain to me in simple terms
>>>why s**2 is less accurate than s*s.  I can sort of intuitively
>>>appreciate that that would be true, but but might like just a little
>>>more detail.
>>>
>>>
>>>      
>>>
>>I don't know that it has to be *less* accurate, although it's unlikely
>>to be more accurate since s*s should be nearly as accurate as you get
>>with floating point. Multiplying two complex numbers in numpy is done
>>in the most straightforward way imaginable:
>>
>>   result.real    = z1.real*z2.real - z1.imag*z2.imag
>>   result.image = z1.real*z2.imag + z1.imag*z2.real
>>
>>The individual results lose very little precision and the overall
>>result will be nearly exact to within the limits of floating point.
>>
>>On the other hand, s**2 is being calculated by a completely different
>>route. Something that will look like:
>>
>>   result = pow(s, 2.0)
>>
>>Pow is some general function that computes the value of s to any
>>power. As such it's a lot more complicated than the above simple
>>expression. I don't think that there's any reason in principle that
>>pow(s,2) couldn't be as accurate as s*s, but there is a tradeoff
>>between accuracy, speed and simplicity of implementation.
>>    
>>
>
>On a close tangent, I had a patch at one point for Numeric (never
>committed) that did pow(s, 2.0) (= s**2) actually as s*s at the C level (no
>pow), which helped a lot in speed (currently, s**2 is slower than s*s).
>
>I should have another look at that. The difference is speed is pretty
>bad: for an array of 100 complex elements, s**2 is 68.4 usec/loop as
>opposed to s*s with 4.13 usec/loop on my machine.
>  
>
Python's complex object also special cases integer powers. Which is why 
you won't see the inaccuracy that started this thread using basic 
complex objects.

However, I'm not convinced this is a good idea for numpy. This would 
introduce a discontinuity in a**b that could cause problems in some 
cases. If, for instance, one were running an iterative solver of some 
sort (something I've been known to do), and b was a free variable, it 
could get stuck at b = 2 since things would go nonmonotonic there. I 
would recomend that we just prominently document that x*x is faster and 
more accurate than x**2 and that people should use x*x where that's a 
concern.

-tim