[Discuss] Test for positive gives wacko results for underflowed values for the gcc compiler suite (including the gcc compiler, and g77 and gfortran compilers)

[Alan W. Irwin]

I have just run into something really wacko for some fortran code (see
corresponding C and python results below) I was
debugging.  Here is a test routine that exhibits the behaviour.
(The corresponding C and python routines are below.)

       real*8 x, xscale
       x = 1.d-200
       xscale = 1.d+250
       write(*,'(1p3d15.5)') x, xscale, x/xscale
       write(*,*) x/xscale.gt.0.d0
       write(*,*) x/xscale.gt.1.d-305
       write(*,*) log(x/xscale)
       end

Here is the strange result
I get from this code for both g77 and gfortran.

     1.000000000000000-200    9.999999999999999+249    0.000000000000000E+00
  T
  F
  -INF

The point is that x/xscale = 1.d-450 underflows (i.e., the above output
gives a zero result for x/xscale and -infinity for log(x/xscale)) as
expected. From http://steve.hollasch.net/cgindex/coding/ieeefloat.html it
appears the log (base 10) of the underflow limit ranges from -307.6
(normalized) to -323.2 (unnormalized) so x/xscale should always underflow
(by far). Nevertheless, the test on x/xscale.gt.0.d0 yields true rather than
the expected false while the test on x/xscale.gt.1.d-305 yields the expected
false.

I ran across this floating-point issue in complicated FreeEOS code where I
do one thing involving log(x/xscale) if x/scale is .gt. 0.d0, and another
thing not involving logs if x/xscale is zero.  The workaround is of course
to use a test like x/xscale.gt.1.d-307 (where 1.d-307 is a bit larger than
the above double-precision underflow limits), but I would still like to know
why the x/xscale.gt.0.d0 test is not giving false under the above
circumstances.

It turns out the same issue occurs in C as well. The test C code is

#include <stdio.h>
#include <math.h>
int main(void) {
    double x=1.e-200, xscale=1.e+250;
    double y;
    int greater;
    printf("%15.5e %15.5e %15.5e \n", x, xscale, x/xscale);
    greater = x/xscale > 0.e0;
    printf("%5i  \n", greater);
    greater = x/xscale > 1.e-305;
    printf("%5i  \n", greater);
    y = log(x/xscale);
    printf("%15.5e \n", y);

    return(0);
}


and the result is

    1.00000e-200    1.00000e+250     0.00000e+00
     1
     0
            -inf

However, python does not have this issue:

The code is

#!/usr/bin/python
from math import *
x = 6.317532735950019e-126
xscale = 1.0400000000000001e+200
print x, xscale, x/xscale
print x/xscale > 0.e0
print x/xscale > 1.e-305
print log(x/xscale)

and the result is

1e-200 1e+250 0.0
False
False
Traceback (most recent call last):
   File "test.py", line 7, in ?
     print log(x/xscale)
OverflowError: math range error

(note the two False results.)

Anybody have a clue about what is going on here between the bad results for
fortran and C on the one hand and good results for python on the other?

Is this a gcc (and g77 and gfortran) issue?  I would appreciate it if
somebody would run the above C test codes with different C
compilers/platforms.

Alan
__________________________
Alan W. Irwin

Astronomical research affiliation with Department of Physics and Astronomy,
University of Victoria (astrowww.phys.uvic.ca).

Programming affiliations with the FreeEOS equation-of-state implementation
for stellar interiors (freeeos.sf.net); PLplot scientific plotting software
package (plplot.org); the libLASi project (unifont.org/lasi); the Loads of
Linux Links project (loll.sf.net); and the Linux Brochure Project
(lbproject.sf.net).
__________________________

Linux-powered Science
__________________________