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

[Murray Strome]

Alan W. Irwin wrote:
> 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
> __________________________
I don't really know the answer to this one, but just as a guess, is it 
possible that Python would automatically use a larger word size for 
calculating these things (e.g. 64 bit or 128 bit floating point as 
opposed to 32bit)? Might that be the difference?

Murray