On Fri, 18 Oct 2002, Stavros Macrakis wrote:
> > computational anomaly in the double precision case
>
> Computational anomaly is right!
>
> [omitted]
You're seeing the fact that 1) the x86 fpu normally runs in extended
precision, 2) the x86 fpu has multiple registers so some calculations
are done there without storing back into memory, and 3) C compilers are
allowed to optimize certain operations involving constants.
I modified your test code for Linux where I can easily set the fpu mode.
If I do
gcc -o try try.c
./try
then I see what you see. If I do
./try d
then I see what lisp does. Either lisp is not maintaining variables on
the fpu stack or it is setting the fpu control word.
Generally running the fpu in extended precision is a good idea but there
are sometimes problems.
The Linux C library has the ISP C99 %a format spec so you can see that
the differences in the numbers are small.
BTW, if you change the definition of VOLATILE below and compile with
optimization, you'll see what else the compiler is allowed to do.
10
===================================================================
/* try.c */
#include <stdio.h>
#include <fpu_control.h>
void set_double_precision (void)
{
fpu_control_t mode;
fputs ("Setting double precision\n", stderr);
_FPU_GETCW (mode);
mode &= ~_FPU_EXTENDED;
mode |= _FPU_DOUBLE;
_FPU_SETCW (mode);
return;
}
void set_single_precision (void)
{
fpu_control_t mode;
fputs ("Setting single precision\n", stderr);
_FPU_GETCW (mode);
mode &= ~_FPU_EXTENDED;
mode |= _FPU_DOUBLE;
_FPU_SETCW (mode);
return;
}
void set_extended_precision (void)
{
fpu_control_t mode;
fputs ("Setting extended precision\n", stderr);
_FPU_GETCW (mode);
mode |= _FPU_EXTENDED;
_FPU_SETCW (mode);
return;
}
#if 1
#define VOLATILE volatile
#else
#define VOLATILE
#endif
int main (int argc, char *argv[])
{
double VOLATILE x;
double VOLATILE y;
double VOLATILE m51;
double VOLATILE y51;
double VOLATILE z51;
double VOLATILE m500001;
double VOLATILE y500001;
double VOLATILE z500001;
if (argc > 1)
{
unsigned c = (unsigned char) argv[1][0];
if (c == 's' || c == 'S')
set_single_precision ();
else if (c == 'd' || c == 'D')
set_double_precision ();
else if (c == 'e' || c == 'E' || c == 'x' || c == 'X')
set_extended_precision ();
else
fprintf (stderr, "unknown argument: \"%s\"\n", argv[1]);
}
x = 1.0 + 1.1105e-16;
y = x - 1.0;
printf ("x = %.17e (%a), y = %.17e (%a)\n", x, x, y, y);
m51 = 0.51;
y51 = y * m51;
z51 = 1.0 + y51;
printf ("m = %.17e (%a), z = %.17e (%a)\n", m51, m51, z51 - 1,
z51 - 1);
m500001 = 0.500001;
y500001 = y * m500001;
z500001 = 1.0 + y500001;
printf ("m = %.17e (%a), z = %.17e (%a)\n", m500001, m500001,
z500001 - 1, z500001 - 1);
return 0;
}