[R] floating point question

[ripley@stats.ox.ac.uk]

On Fri, 31 Jan 2003, Thomas Lumley wrote:

> 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

(Sparc FPUs have some guard digits, as far as I recall, and the fine 
details do vary by hardware, or at least they used to and we still have 
some of the machines that varied running.)

> 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.

I used to think that!  But when it became possible to build R for Windows
not chopping in-FPU calculations to 53 bits, the results were more
consistent, and we changed over.  At that point most of the
inconsistencies were in printing, hence my remark earlier.  The print
routines in Windows `libc' (really msvcrt.dll) does lose a digit or two
quite often, quite possibly to avoid giving meaningless ones (which glibc
is happy to if asked).


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595