[Numpy-discussion] Re: indexing problem

[David M. Cooke]

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.

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|cookedm at physics.mcmaster.ca