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

[Raymond Toy]

>>>>> "Camm" == Camm Maguire <address@hidden> writes:

    Camm> "Stavros Macrakis" <address@hidden> writes:
    >> Camm,
    >> 
    >> The new Epsilon calculations look excellent.
    >> 
    >> But I am puzzled by the value 1.1107651257113995E-16 (call that Q765)
    >> for double-float-epsilon.  When running in Lisp, the correct value is in
    >> fact 1.1102230246251568d-16 (call that Q223).  That is, 1+Q223-1 is
    >> non-zero.
    >> 
    >> I have implemented your algorithm both in C and in Lisp.  In Lisp, it
    >> gives the correct result (Q223).  But in C, it gives Q765, UNLESS you
    >> don't force intermediate results to memory.  But surely in Lisp, all the
    >> intermediate results are in fact going to memory!  And the binary
    >> representation of Q765 looks bizarre.
    >> 
    >> It is not your code that's the problem -- I tried a different algorithm
    >> and got the same result.
    >> 
    >> Any clue what is going on here?  Or am I just hallucinating?
    >> 

Look at the 2 epsilon values:

(write (integer-decode-float 1.1107651257113995d-16) :base 16)
10020000000001
4505798650626049
* (write (integer-decode-float double-float-epsilon) :base 16)
10000000000001
4503599627370497
* 

So the two values basically differ by a single bit set near the MSB of
the number.  I think that if take the Q765 value and add it (by hand,
not machine) to 1 it will properly round to a number that is different
from 1.  There will be a bit at bit 53 and one at bit 63 so that when
that 80-bit number is rounded, it rounds up because of the bit at bit
63.  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.

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

Ray