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

[Stavros Macrakis]

Raymond--

I understand your theory, but I don't think it's correct.

Here's my understanding of IEEE floating point.  Suppose we're working
with a 5-bit stored mantissa and an 8-bit register mantissa.  When we
add .xxxxx and .00010000, the register result is .xxxxx100, which should
be rounded according to the 'even' rule.  Now what if we add .xxxxx and
.00010001?  Well, if there were no sticky bit, we would still get
.xxxxx100 in the register (since the final 1 has been shifted off),
which again would round to even when stored.  But if we have a sticky
bit (which remembers if anything non-zero has been shifted off), then
the register sum is .xxxxx101.  When we store .xxxxx101, we need to
round to five bits, and since 101 is greater than 100, we round up.  So
we shouldn't be losing the round-up quality when storing.  After all, it
would be a pity if despite the elaborate IEEE machinery, you couldn't
get a simple sum of two floats to round correctly except in intermediate
register results....

Or have I gotten something wrong here?

Anyway, your theory does not explain the Lisp facts.  How is it that
Emacs Lisp and GCL manage to round their results correctly -- remember,
they're written in C, too -- but native C does not?

I looked at the source of Emacs Lisp, and it takes no special
precautions.  I did notice, though, that its additions happen using
accum += term.  So I tried replacing z=x+y with z=x; z+=y; .  No luck.
I also tried different optimization levels, and different settings
for -ffloat-store2S

Perhaps there is some special gcc flag or IEEE floating global mode
involved?

       -s