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

To avoid the weird result I have documented with my simple test.c code, the
numbers taken out of the 80-bit register should always be properly
normalized to 64-bits before any comparisons are done, and it appears gcc is
skipping this step for comparisons like x/xscale > 0.e0.  Note, 80-bit
floating-point division should take a very long time compared to a
normalization operation so skipping the normalization step cannot be
justified on the basis of efficiency. If the skipped-normalization
hypothesis is correct, then this is a fundamental bug in gcc on Intel
hardware.

One might be tempted to think the -ffloat-store gcc compilation option is
relevant to the problem.  Here is what info gcc says about that:

This option prevents undesirable excess precision on machines such as the
68000 where the floating registers (of the 68881) keep more precision than a
double' is supposed to have.  Similarly for the x86 architecture.  For most
programs, the excess precision does only good, but a few programs rely on
the precise definition of IEEE floating point.  Use -ffloat-store' for such
programs, after modifying them to store all pertinent intermediate
computations into variables.

That last sentence may be the killer.  What I am complaining about is the
case where x/xscale is not stored in an intermediate variable, and of course
I am complaining about the range rather than the precision.  Anyhow, the
-ffloat-store gcc option made no difference to my test case.

If the skipped normalization hypothesis is correct, gcc (or any compiler) on
non-Intel platforms without the 80-bit registers should give correct results
for my test case. Does anybody have test results for such platforms? Also,
from Peter's test it looks like the windows C compiler skips the
normalization step on Intel hardware, and it would be interesting to see if
the Intel compiler also has the same issue.

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
__________________________