[R] exactly representable numbers

[Prof Brian Ripley]

On Mon, 11 Sep 2006, Duncan Murdoch wrote:

> On 9/11/2006 6:15 AM, Robin Hankin wrote:
> > Hi Sundar
> > 
> > 
> > thanks for this.  But I didn't make it clear that I'm interested in  
> > extreme numbers
> > such as 1e300 and 1e-300.

That's not relevant, unless you are interested in extremely small numbers.

> > Then
> > 
> >  > f(1e300)
> > [1] 7.911257e+283

(That is inaccurate)

> > is different from
> > 
> > 1e300*.Machine$double.eps

Yes, since 1e300 is not a power of two.  However, Sundar is right in the 
sense that this is an upper bound for normalized numbers.

f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2.
This gives you the answer:  .Machine$double.neg.eps * 2^floor(log2(x))

Similarly for going below (but carefully as you get an extra halving on 
the powers of two).

These results hold for all but denormalized numbers (those below 1e-308).


> > [I'm interested in the gap between successive different exactly  
> > representable
> > numbers right across the IEEE range]
> 
> I'm not sure the result you're looking for is well defined, because on 
> at least the Windows platform, R makes use of 80 bit temporaries as well 
> as 64 bit double precision reals.  I don't know any, but would guess 
> there exist examples of apparently equivalent formulations of your 
> question that give different answers because one uses the temporaries 
> and the other doesn't.

Not at R level.  For something to get stored in a real vector, it will be 
a standard 64-bit double.

> But in answer to your question:  a bisection search is what I'd use: 
> you start with x+delta > x, and you know x+0 == x, so use 0 and delta as 
> bracketing points.  You should be able to find the value in about 50-60 
> bisections if you start with delta == x, many fewer if you make use of 
> the double.eps value.  Here's my version:  not tested too much.
> 
> f <- function(x) {
>    u <- x
>    l <- 0
>    mid <- u/2
>    while (l < mid && mid < u) {
>      if (x < x + mid) u <- mid
>      else l <- mid
>      mid <- (l + u)/2
>    }
>    u
> }
> 
>  > f(1e300)
> [1] 7.438715e+283
>  > 1e300 + 7.438715e+283  > 1e300
> [1] TRUE
>  > 1e300 + 7.438714e+283  > 1e300
> [1] FALSE
> 
> 
> Duncan Murdoch
> 
> > 
> > 
> > 
> > rksh
> > 
> > 
> > 
> > 
> > 
> > On 11 Sep 2006, at 10:50, Sundar Dorai-Raj wrote:
> > 
> >>
> >> Robin Hankin said the following on 9/11/2006 3:52 AM:
> >>> Hi
> >>> Given a real number x,  I want to know how accurately R can  
> >>> represent  numbers near x.
> >>> In particular, I want to know the infimum of exactly representable
> >>> numbers greater than x, and the supremum of exactly representable   
> >>> numbers
> >>> less than x.  And then the interesting thing is the difference   
> >>> between these two.
> >>> I have a little function that does some of this:
> >>> f <- function(x,FAC=1.1){
> >>>    delta <- x
> >>> while(x+delta > x){
> >>>    delta <- delta/FAC
> >>> }
> >>> return(delta*FAC)
> >>> }
> >>> But this can't be optimal.
> >>> Is there a better way?
> >> I believe this is what .Machine$double.eps is. From ?.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'.
> >>
> >> See also .Machine$double.neg.eps. Is this what you need?
> >>
> >> --sundar
> > 
> > --
> > Robin Hankin
> > Uncertainty Analyst
> > National Oceanography Centre, Southampton
> > European Way, Southampton SO14 3ZH, UK
> >   tel  023-8059-7743
> > 
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

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