[R] floating point question

[Thomas Lumley]

On 31 Jan 2003, Peter Dalgaard BSA wrote:

> Deepayan Sarkar <deepayan at stat.wisc.edu> writes:
>
> > Might have something to do with .Machine$double.eps on the respective
> > machines.
> >
> > From help(.Machine),
> >
> > double.eps: the smallest positive floating-point number `x' such that
> >           `1 + x != 1'.  It equals `base^ulp.digits' if either `base'
> >           is 2 or `rounding' is 0;  otherwise, it is `(base^ulp.digits)
> >           / 2'.
>
> No, it's trickier than that. I think it's due to the guard digits that
> intel FPUs use. Both intel and Sun have 53bit mantissas in double pr.
> but intel stores intermediate results to 64 bit precision, so you get
> some double rounding effects. Essentially, during alignment, before
> adding to 1, 2^-53(1+2^-52) is getting rounded up to binary
>

There's some discussion of double rounding in a Appendix to David
Goldberg's "What every computer scientist should know about floating point
arithmetic", widely available on the Web.

However, as Brian Ripley points out, we have no control and there isn't
any point worrying about it.

Incidentally, this may well be a reason for the irritating habit of DLLs
in setting the floating point precision to 53 bits (well, the irritating
bit is not changing it back). This makes the results of computations
more nearly independent of where the numbers end up being stored.

	-thomas