[Numpy-discussion] Floating point question

[Anne Archibald]

On 02/03/2009, Gideon Simpson <simpson at math.toronto.edu> wrote:
> I recently discovered that for 8 byte floating point numbers, my
>  fortran compilers (gfortran 4.2 and ifort 11.0) on an OS X core 2 duo
>  machine believe the  smallest number 2.220507...E-308.  I presume that
>  my C compilers have similar results.
>
>  I then discovered that the smallest floating point number in python
>  2.5 is 4.9065...E-324.  I have been using numpy to generate data,
>  saving it with savetxt, and then reading it in as ASCII into my
>  fortran code.  Recently, it crapped out on something because it didn't
>  like reading it a number that small, though it is apparently perfectly
>  acceptable to python.
>
>  My two questions are:
>
>  1.  What is the best way to handle this?  Is it just to add a filter
>  of the form
>
>         u = u * ( np.abs(u) > 2.3 e-308)
>
>  2.  What gives?  What's the origin of this (perceived) inconsistency
>  in floating points across languages within the same platform?
>
>  I recognize that this isn't specific to Scipy/Numpy, but thought
>  someone here might have the answer.

What's happening is that numbers like 1e-310 are represented by
"denormals". If we use base 10 to explain, suppose that floating point
numbers could only have five digits of mantissa and two digits of
exponent:

1.3000 e 00
2.1000 e-35
1.0000 e-99

Now what to do if you divide that last number in half? You can write it as:

0.5000 e-99

But this is a bit of an anomaly: unlike all normal floating-point
numbers, it has a leading zero, and there are only four digits of
information in the mantissa.  In binary it's even more of an anomaly,
since in binary you can take advantage of the fact that all normal
mantissas start with one by not bothering to store the one. But it
turns out that implementing denormals is useful to provide graceful
degradation as numbers underflow, so it's in IEEE. It appears that
some FORTRAN implementations cannot handle denormals, which is giving
you trouble. It's usually fairly safe to simply replace all denormals
by zero.

Anne