Re: [Numpy-discussion] integer power in scalarmath: how to handle overflow?

[David M. C. <co...@ph...>]

On Tue, 6 Jun 2006 10:44:20 +0200
David Douard <dav...@lo...> wrote:

> On Mon, Jun 05, 2006 at 05:10:23PM -0400, David M. Cooke wrote:
> > I just ran into the fact that the power function for integer types
> > isn't handled in scalarmath yet. I'm going to add it, but I'm
> > wondering what people want when power overflows the integer type?
> > 
> > Taking the concrete example of a = uint8(3), b = uint8(10), then
> > should a ** b return
> > 
> > 1) the maximum integer for the type (255 here)
> > 2) 0
> > 3) upcast to the largest type that will hold it (but what if it's
> >    larger than our largest type? Return a Python long?)
> > 4) do the power using "long" like Python does, then downcast it to
> > the type (that would return 169 for the above example)
> > 5) something else?
> > 
> > I'm leaning towards #3; if you do a ** 10, you get the right
> > answer (59049 as an int64scalar), because 'a' is upcasted to
> > int64scalar, since '10', a Python int, is converted to that type.
> > Otherwise, I would choose #1.
> 
> I agree, #1 seems the better solution for me.
> 
> BTW, I'm quite new on this list, and I don't know is this has already
> been discussed (I guess I has): why does uint_n arithmetics are done
> in the  Z/(2**n)Z field (not sure about the maths correctness here)?
> I mean:
> >>> a = numpy.uint8(10)
> >>> a*a
> 100
> >>> a*a*a # I'd like to have 255 here
> 232
> >>> 1000%256
> 232
> 

History, and efficiency. Detecting integer overflow in C portably requires
doing a division afterwards, or splitting the multiplication up into parts
that won't overflow, so you can see if the sum would. Both of those
options are pretty slow compared with multiplication.

Now, mind you, our scalar types *do* check for overflow: they use a
larger integer type for the result (or by splitting it up for the
largest type). So you can check for overflow by setting the overflow
handler:

>>> seterr(over='raise')
{'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore',
'under': 'ignore'}
>>> int16(32000) * int16(3)
Traceback (most recent call last):
      File "<stdin>", line 1, in ?
FloatingPointError: overflow encountered in short_scalars

Note that the integer array types don't check, though (huh, maybe they
should). It's easy enough to use the multiply routine for the power, so
you'll get overflow checking for free.

> It would be really a nice feature to be able (by the mean of a numpy
> flag or so) to have bound limited uint operations (especially when
> doing image processing...).

If you want to supply a patch ... :-)

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...