[Matrix-SIG] Problem with integer arrays?

[dlr@postbox.ius.cs.cmu.edu]

Charles Waldman wrote:
>> >>> a = Numeric.array((0,0))
>> >>> a**2
>> Traceback (innermost last):
>>   File "<stdin>", line 1, in ?
>> OverflowError: math range error

In umathmodule.c, the function powll() has a range check:

> /* Overflow check: overflow will occur if log2(abs(x)) * n > nbits. */

On my machine, the math range error is happening when we take 
log2(0) = -inf.  I think the best fix is to put a test around the
overflow check so that it's not executed if x == 0.  I've modified my
version of umathmodule.c, and it seems to work fine (diffs appended).

Also, I noticed this strange behavior yesterday:

>>> import Numeric
>>> Numeric.__version__
'1.4'
>>> a = Numeric.ones(5)
>>> a
array([1, 1, 1, 1, 1])
>>> (-1) * a
array([-1, 1, -1, 1, -1])

I've made and mildly tested a change to umathmodule.c which I think
fixes this too (see below).

Of course, I'd welcome corrections or comments.
David LaRose

-------------------------------------
Note: This diff is for Numeric 1.4, from the LLNLPython5 distribution

diff umathmodule.c umathmodule.orig.c
261,271c261,263
< 
<       /******Added by dlr@cs.cmu.edu, 10/6/98**********/
<       if(x != 0) {   /* Range check unnecessary & dangerous if x == 0 */
<       /******End of Addition by dlr@cs.cmu.edu, 10/6/98**********/
<         logtwox = log10 (fabs ( (double) x))/log10 ( (double) 2.0);
<         if (logtwox * (double) n > (double) nbits)
<           PyErr_SetString(PyExc_ArithmeticError, "Integer overflow in power.");
<       /******Added by dlr@cs.cmu.edu, 10/6/98**********/
<       }  /* if(x != 0) */
<       /******End of Addition by dlr@cs.cmu.edu, 10/6/98**********/
< 
---
>       logtwox = log10 (fabs ( (double) x))/log10 ( (double) 2.0);
>       if (logtwox * (double) n > (double) nbits)
>               PyErr_SetString(PyExc_ArithmeticError, "Integer overflow in power.");
439,441d430
<     /********Added by dlr, 10/6/98*******/
<     s = 1;  /* reset sign flag */
<     /********End of addition by dlr, 10/6/98*******/
546,548d534
<     /********Added by dlr, 10/6/98*******/
<     s = 1;  /* reset sign flag */
<     /********End of addition by dlr, 10/6/98*******/
653,655d638
<     /********Added by dlr, 10/6/98*******/
<     s = 1;  /* reset sign flag */
<     /********End of addition by dlr, 10/6/98*******/