[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]

On 2007-06-25 20:01-0700 Alan W. Irwin wrote:

> I now think I have a reasonable working hypothesis for what is going on.
>
> These weird results where x/xscale is greater than 0.e0 while y=x/xscale is
> not (see my latest test.c code) may be related to an old gcc bug (see
> http://gcc.gnu.org/ml/gcc-bugs/2000-06/msg00113.html) that only appears for
> Intel hardware.  Intel hardware has an 80-bit floating-point word (64-bit
> mantissa, 15-bit exponent, and one bit for the sign, see
> http://webster.cs.ucr.edu/AoA/DOS/ch14/CH14-1.html#HEADING1-52). 15 bits is
> 4 bits more than the 11 bits of exponent in ieee 754 (double precision) so
> the underflow limit must be something like 10.d-1200. Thus the x/xscale
> value in the 80-bit register cannot possibly underflow (even for the extreme
> test case I used with double x=1.e-305, xscale=1.e+305;) while if you store
> y = x/xscale, then the value in y must fit into a 11-bit exponent and is
> guaranteed to underflow.

Hmm. That underflow limit for the 80-bit register should be approximately
(10^{-305})^16 = 10^{-1488400}, not 10^{-1200} as stated above.

Here is a test programme and results that are strongly consistent with the
hypothesis that when x is slightly greater than the 64-bit IEEE 754
underflow limit, then x^16 does not underflow the 80-bit register while x^32
does.

int main(void) {
    double x=1.e-305;
    double y;
    int greater;
    y = x*x;
    greater = y > 0.e0;
    printf("%10s %5.1d\n", "y=x^2", greater);
    greater = (x*x) > 0.e0;
    printf("%10s %5.1d\n","x^2", greater);
    greater = (x*x*x*x) > 0.e0;
    printf("%10s %5.1d\n","x^4", greater);
    greater = (x*x*x*x*x*x*x*x) > 0.e0;
    printf("%10s %5.1d\n","x^8", greater);
    greater = (x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x) > 0.e0;
    printf("%10s %5.1d\n","x^16", greater);
    greater = (x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x) > 0.e0;
    printf("%10s %5.1d\n","x^32", greater);
}

results:

      y=x^2     0
        x^2     1
        x^4     1
        x^8     1
       x^16     1
       x^32     0

The results, of course, are also consistent with the hypothesis that gcc has
failed to normalize the 80-bit register results to IEEE 754 64-bit before
doing comparisons.  To my mind that is a bug since people do test, for
example, the argument of log to be greater than zero to avoid -inf results
(which is how I found this bug in the first place for a result that was
underflowing in 64 bit, but not in 80 bit).  On Intel, Daniel has effectively
done the first two tests confirming my Intel results, and also shown on
PowerPC that the x^2 case yields the correct false (0) rather than the
incorrect true (1) result you get for Intel hardware.

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
__________________________