> computational anomaly in the double precision case
Computational anomaly is right!
The following C code (gcc 2.95.3-5 on Athlon, W2k) (see full test code
below):
double x,y,z;
x=1.0+1.1105e-16; FMEM(&x); /* Very near epsilon */
y=x-1.0; FMEM(&y);
z=1.0+(y*0.500001); FMEM(&z); /* Should be > 1.0 */
printf("%.17e %.17e %.17e\n",x,y,z-1);
produces
1.00000000000000022e+00 2.22044604925031308e-16
0.00000000000000000e+00
but the equivalent Lisp code (in both GCL 2.3.6 and Emacs Lisp 20.7.1
(i386-*-nt5.0.2195))
(setq x (+ 1.0 1.1105e-16))
(setq y (- x 1))
(setq z (+ 1 (* y .500001)))
(princ (format "%.17e %.17e %.17e\n" x y (- z 1))) [or "~,17E"]
gives
1.00000000000000020e+000 2.22044604925031310e-016
2.22044604925031310e-016
How is this possible? (If the IEEE rounding mode were set to round
down, then 1+(.51*y)-1 would also be 0.0, so it's not that....)
After all, both GCL and Emacs Lisp are built with gcc, and the Lisp
answer looks like the correct IEEE floating point result. And in
interpreted Lisp, the intermediate results must be going to memory. Or
am I missing something blindingly obvious?
-s
(in test.c)
main()
{
double x,y,
m51,y51,z51,
m500001,y500001,z500001;
x=1.0+1.1105e-16; FMEM(&x);
y=x-1.0; FMEM(&y);
printf("x = %.17e y = %.17e\n",x,y);
m51=0.51;
y51=y*m51; FMEM(&y51);
z51=1.0+y51; FMEM(&z51);
printf("m = %.17e z = %.17e\n",m51,z51-1);
m500001=0.500001;
y500001=y*m500001; FMEM(&y500001);
z500001=1.0+y500001; FMEM(&z500001);
printf("m = %.17e z = %.17e\n",m500001,z500001-1);
}
--------------
(in fmem.c)
FMEM(q) void *q; {;}
--------------
gcc test.c fmem.c -o test; ./test
x = 1.00000000000000022e+00 y = 2.22044604925031308e-16
m = 5.10000000000000009e-01 z = 2.22044604925031308e-16
m = 5.00001000000000029e-01 z = 0.00000000000000000e+00