>>>>> "Camm" == Camm Maguire <camm@enhanced.com> writes:
Camm> "Stavros Macrakis" <stavros.macrakis@verizon.net> 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