numpy, overflow, inf, ieee, and rich comparison

[Tim Peters]

[Huaiyu Zhu]
> ...
>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

[Tim Roberts]
> Hmmm.  My Win98 box with Python 2.0b2 gives very different results:
>
> >>> exp(-740)
> 4.150151425066471e-322
> >>> exp(-744)
> 4.9406564584124654e-324
> >>> exp(-745)
> 0.0
> >>> exp(-746)
> 0.0
> >>> exp(-800)
> 0.0
>
> Interestoid #1 is that I don't get a range error; I assume something was
> fixed in 2.0.

No, Python inherits its exp results from the platform C implementation.
Huaiyu is simply suffering from not running Windows <wink>.

> Interestoid #2 is that my exp(-744) == your exp(-745).  I can't see
> that either result is correct; both my Casio calculator (using
> logarithms) and trusty old "bc", exp(-744) should be about 7.67e-324.
>
> I assume we've reached the lower limit of an IEEE 64-bit float?

Good eye!  You can get more insight by picking apart the bits:

>>> math.frexp(math.exp(-744))
(0.5, -1073)
>>> math.ldexp(0.5, -1073)
4.9406564584124654e-324
>>> math.ldexp(.5, -1074)
0.0
>>>

That is, the result is (exactly) 2.**-1074, which is the very limit of 754's
"denormal" range.  There's only 1 bit of precision left, and 4.94 *is* 7.67
to within one bit <ahem>.  The closest representable number on the other
side is 2.**-1073:

>>> math.ldexp(2, -1073)
1.9762625833649862e-323
>>>

Neither of {4.94e-324, 1.98e-323} is a reasonable approximation to
7.67e-324, but they're the best that can be done with a 754 double, and FWIW
the Windows result is the better one (probably not an accident, BTW:  while
MS people have shown no particular ability to sling floating-point, the
Pentium chip has much of this stuff built in, and the Intel people knew what
they were doing).

unfortunately-much-of-your-intel-fp-hardware-sits-unused-ly y'rs  - tim