[Numpy-discussion] RE: Python 2.2 seriously crippled for numerical computation?

[Huaiyu Zhu]

On Tue, 05 Mar 2002 19:50:10 -0500, Tim Peters <tim.one at comcast.net> wrote:
>
>>    Since both 1e-200**2 and 1e200**2 produce the same errno all the time,
>>    but Python still raises OverflowError for 1e200**2 when linked with
>>    -lieee, there must be a separate mechanism at work.
>
> ... What differs is when
>platform pow(x, y) does not set errno.  In that case, Python synthesizes
>errno=ERANGE if the pow() result equals +- platform HUGE_VAL.

Thanks!  Thanks!  Thanks!  That's exactly what I was looking for.  The
OverflowError was produced by the line commented out below:

	errno = 0;
	PyFPE_START_PROTECT("pow", return NULL)
	ix = pow(iv, iw);
	PyFPE_END_PROTECT(ix)
	  /*Py_SET_ERANGE_IF_OVERFLOW(ix); */
	if (errno != 0) {
		/* XXX could it be another type of error? */
		PyErr_SetFromErrno(PyExc_OverflowError);
		return NULL;
	}

With this little change, I got the following answer consistantly:

>>> 1e200**2
inf
>>> (1e200*1j)**2
(-inf+nanj)

instead of the inexplicable:

>>> 1e200**2
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
OverflowError: (34, 'Numerical result out of range')
>>> (1e200*1j)**2
(-inf+nanj)

>
>>  What is that and how can I override it?
>
>Sorry, you can't override it.

I just did.  So I guess I must have been communicating terribly.  But
anyway, all is fine and the sun is shining now.  :-)

Why do I want this?  I believe it is impossible for computation to carry its
own numerical analysis without prohibitive cost (think interval computation
in the 1970s).  If the user knows his numerical analysis, you don't want to
short circuit him, and if he does not, it's hopeless anyway.  Overflow error
should only be raised when the return value is out of bound AND not
represented by well defined out of bound values such as inf and nan.

I'm not going to write any lengthy rationale for this, as my itch has been
scratched.  The following example illustrate my point (Consider: How do you
know that y=x**2 is not used in exp(-y)?  And what should happen if this is
inside a matrix where you don't want try/except for each element?  Answer:
you can't.  Numerical analysis in general can only be performed globally.)

>>> from math import exp
>>> def f(x): return exp(-x**2)
... 
>>> f(1)
0.36787944117144233
>>> f(5)
1.3887943864964021e-11
>>> f(1e100)
0.0
>>> f(1e200)
0.0

versus

>>> from math import exp
>>> def f(x): return exp(-x**2)
... 
>>> f(1)
0.36787944117144233
>>> f(5)
1.3887943864964021e-11
>>> f(1e100)
0.0
>>> f(1e200)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "<stdin>", line 1, in f
OverflowError: (34, 'Numerical result out of range')


Hauiyu