numpy, overflow, inf, ieee, and rich comparison

[Huaiyu Zhu]

I'm stumbling over a seemingly trivial problem.  Following are some
descriptions of the problem, various attemps at solving it, and some
comments. 

We know that exp(-1000) is approximately zero to machine precision.
However, on my machine

>>> from math import *
>>> exp(-745)
4.9406564584124654e-324
>>> exp(-746)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
OverflowError: math range error

Now, for this simple example, it is easy to add an if statement.  But for a
Numeric.array, which has the same problem, an if statement obviously does
not work, in the sense that it defeats all the speed gains.

So the first question is: Can this be avoided?  My understanding is that
this is the issue of the underlying C library and cannot be avoided unless
we are willing to rewrite it.

The second question is: Can I change all the elements of the array that's
smaller than -745 into -inf?  It seems that inf only works under IEEE
arithmetic.  But since all the OSes I'm going to use have it, I'll try it
anyway.

Now the third question is, how do I do that?  The obvious way is

x = less(x,-745)*-inf + less(x,745)*greater(x,-745)*x + greater(x,745)*inf

But this does not work, because x*inf generates math range error if x
contains zero.  Interestingly, for plain float, 0*inf==nan, which seems more
reasonable, although it still does not solve the problem.

At this moment it seems that the only way available to me is to search the
whole array and replace the out of bound elements one by one.  This is dozens
of times less efficient.

This incidentally brings to a further observation concerning lack of the ?:
operator in Python.  If we had that, and with some rich comparison
operators, we could write the above as

x = (x<-745) ? -inf : (x>745) ? inf : x

This operator could be written in a C extension module that only goes over
the indices once.

So, all my ramblings aside, maybe someone can point out the things I've
overlooked?  For my needs now just something working on Linux would be
enough, but if it works generally I'll put it in MatPy as well.

Huaiyu