[Gcl-devel] Re: [Maxima] Re: Epsilon calculation

[Raymond Toy]

>>>>> "Stavros" == Stavros Macrakis <address@hidden> writes:

    >> If that extra bit weren't there, the 80-bit number will round-to-even
    >> and therefore round down, back to 1, which isn't what is wanted.

    Stavros> It's been a few years since I read the spec, but as I remember it, 
the
    Stavros> IEEE sticky bit takes care of situations like this.  If anything
    Stavros> non-zero gets shifted off the right of the mantissa, the sticky bit
    Stavros> remains turned on to ensure correct rounding.

Yes, this is correct.  The 80-bit number is correctly rounded.
However, when you store that 80-bit number away, there's an additional
rounding from that 80-bit number to a 64-bit number.  Since the 53rd
bit (54?) is zero, it rounds down and you get 1d0, precisely.

    >> I think this is about the best one can do on an x86 architecture.

    Stavros> Testing in Emacs Lisp and GCL shows that the machine correctly
    Stavros> implements the sticky bit, so what is happening in C?

I think to be sure you'd have to look at the disassembly of the
relevant code in GCL and C.  GCL may be storing intermediate results
to memory but C may not.  Maybe if you compile the C code with -O0 (no
optimization at all), you'll get the answer you want?

I also note the Clisp on x86 now uses that number for it's
double-float epsilon.

Ray