Re: [R] exactly representable numbers

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Mon 11 Sep 2006 - 11:12:03 GMT

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@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@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. Received on Mon Sep 11 21:16:04 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 11 Sep 2006 - 16:30:31 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.