x86 extended precision hazards: dtoa.c

[David G. Hough at validgh]

As part of the linux-233 project (http://www.validgh.com/linux-233 is updated
from time to time with progress notes)
I've wasted a day trying to debug f2c 
failures on x86-linux that are not present on sparc-linux.   I had been
using netlib fp/dtoa.c and g_fmt.c for f2c's in order to get
reproducible-across-platforms constant base conversion, but the x86-linux
failures went away when I quit using dtoa.c.

My suspicion
was a problem with fp stack overflow or underflow producing NaNs, which can
be caused by inconsistent declarations somewhere in the f2c source code.
And indeed I was able to work around the first failure like this:

#if defined(__linux__) && defined(__i386__)

/* some kind of expression evaluation problem with GCC 2.7.2.2 -
        x86 stack overflows? and ds gets Nan value */

{ double t1,t2,t3,t4,t5;
t1=(d2-1.5);  
t4 = i;    
t2 = t1 * 0.289529654602168;
t5 = t4 * 0.301029995663981;
t3 = t2 + 0.1760912590558;
ds = t3 + t5; 
if (ds != ds) fprintf(stderr," ds %g d2 %g i %d \n",ds,d2,i);
}
#else
        ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;#endif   


but that just postponed the failure.
Then I noticed similar failures with my hacked-up version of f2c running
on SPARC with Apogee compilers.   The symptom is a termination with a message
like "Compiler error line 23864 of : Constant expression yields NaN."
on plain vanilla source code.

On a hunch I went to netlib to see if dtoa and g_fmt had been revised
lately, and sure enough they have.    And to get to the point of this
posting, the first prominent change is a new comment in dtoa.c:

/* On a machine with IEEE extended-precision registers, it is
 * necessary to specify double-precision (53-bit) rounding precision
 * before invoking strtod or dtoa.  If the machine uses (the equivalent
 * of) Intel 80x87 arithmetic, the call
 *      _control87(PC_53, MCW_PC);
 * does this with many compilers.  Whether this or another call is
 * appropriate depends on the compiler; for this to work, it may be
 * necessary to #include "float.h" or another system-dependent header
 * file.
 */

An all-too-familiar warning.    The point relative to the ongoing Java
discussions in this mailing list is not that C compilers and C programs
can't be written to handle IEEE 754 extended precision correctly and
consistently, but that after 15 years they still usually aren't.    In the case
of dtoa.c it probably isn't economic to figure out a way to preserve the
correct operation of the program in the face of unpredictable extended-
precision evaluation over various compilers, releases, and
optimization levels, and so it's easier just to add a warning.

Rather than invent another language in which the same old problems arose,
the Java inventors tried to do something different.

As for my f2c failures, I will see if they go away when I install the revised
dtoa.c and g_fmt.c.