[Numpy-discussion] optimizing power() for complex and real cases

[David M. Cooke]

Tim Hochberg <tim.hochberg at cox.net> writes:
> David M. Cooke wrote:
>
>>[1] which brings up another point. Would using the 3-multiplication
>>version for complex multiplication be good? There might be some
>>effects with cancellation errors due to the extra subtractions...
>>
>>
> I'm inclined to leave this be for now. Both because I'm unsure of the
> rounding issues and because I'm not sure it would actually be faster.
> It has one less multiplication, but several more additions, so it
> would depend on the relative speed add/sub with multiplication and how
> things end up getting scheduled in the FP pipeline. At some point it's
> probably worth trying; if it turns out to be signifigantly faster we
> can think about rounding then. If it's not faster then no need to
> think.

I did some thinking, and looked up how to analyse it. 3M goes like
this:

xy = (a+bi)(c+di) = (ac-bd) + ((a+c)(b+d)-ac-bd)i

Consider x = y = t + i/t, for which x**2 = (t**2-1/t**2) + 2i, then

xy=x^2 = t*t - 1/t*1/t + ((t+1/t)(t+1/t) - t**2 - 1/t**2)i

Consider when t is large enough that (t+1/t)**2 = t**2 in floating
point; then Im fl(xy) will be -1/t**2, instead of 2.

So...let's leave it as is.

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