# [Numpy-discussion] Re: indexing problem

Tim Hochberg tim.hochberg at cox.net
Mon Feb 13 21:18:12 CST 2006

```David M. Cooke wrote:

>Gary Ruben <gruben at bigpond.net.au> writes:
>
>
>
>>Tim Hochberg wrote:
>><snip>
>>
>>
>>>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 don't quite understand the problem here. Tim says Python special
>>cases integer powers but then talks about the problem when b is a
>>floating type. I think special casing x**2 and maybe even x**3 when
>>the power is an integer is still a good idea.
>>
>>
>
>Well, what I had done with Numeric did special case x**0, x**1,
>x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
>exponent was a scalar (so x**y where y was an array wouldn't be). I
>think this is very useful, as I don't want to microoptimize my code to
>x*x instead of x**2. The reason for just scalar exponents was so
>choosing how to do the power was lifted out of the inner loop. With
>that, x**2 was as fast as x*x.
>
>
This is getting harder to object to since, try as I might I can't get
a**b to go nonmontonic in the vicinity of b==2. I run out of floating
point resolution before the slight shift due to special casing at 2
results in nonmonoticity. I suspect that I could manage it with enough
work, but it would require some unlikely function of a**b. I'm not sure
if I'm really on board with this, but let me float a slightly modified
proposal anyway:

1. numpy.power stays as it is now. That way in the rare case that
someone runs into trouble they can drop back to power. Alternatively
there could be rawpower and power where rawpower has the current
behaviour. While the name rawpower sounds cool/cheesy, power is used
infrequently enough that I doubt it matters whether it has these special
case optimazations.

2, Don't distinguish between scalars and arrays -- that just makes
things harder to explain.

3. Python itself special cases all integral powers between -100 and
100. Beg/borrow/steal their code. This makes it easier to explain since
all smallish integer powers are just automagically faster.

4. Is the performance advantage of special casing a**0.5 signifigant?
If so use the above trick to special case all half integral  and
integral powers between -N and N. Since sqrt probably chews up some time
the cutoff. The cutoff probably shifts somewhat if we're optimizing half
integral as well as integral powers. Perhaps N would be 32 or 64.

The net result of this is that a**b would be computed using a
combination of repeated multiplication and sqrt for  real integral and
half integral values of b between -N and N. That seems simpler to
explain and somewhat more useful as well.

It sounds like a fun project although I'm not certain yet that it's a
good idea.

-tim

```