[R] exactly representable numbers

[Duncan Murdoch]

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.
> 
> Then
> 
>  > f(1e300)
> [1] 7.911257e+283
> 
> is different from
> 
> 1e300*.Machine$double.eps
> 
> 
> [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.

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.