[Numpy-discussion] Re: indexing problem

[Tim Hochberg]

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