From: Duncan Murdoch <murdoch_at_stats.uwo.ca>

Date: Mon 11 Sep 2006 - 12:46:25 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.
*

*>> > > f(1e300)
*

*>> > [1] 7.911257e+283
*

*>> >
*

*>> > 1e300*.Machine$double.eps
*

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

*>>
*

>

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 23:02:19 2006

Date: Mon 11 Sep 2006 - 12:46:25 GMT

On 9/11/2006 8:20 AM, Prof Brian Ripley wrote:

> On Mon, 11 Sep 2006, Duncan Murdoch wrote: >

>> On 9/11/2006 6:15 AM, Robin Hankin wrote:

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

>> > Then

>> >

> > (That is inaccurate) >

>> > is different from

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

I'm not sure what is going wrong, but that is too small (on my machine, at least):

> f1 <- function(x) .Machine$double.neg.eps * 2^floor(log2(x))
> f1(1e300)

[1] 7.435085e+283

> 1e300 + f1(1e300) == 1e300

**[1] TRUE
**
Notice the difference in the 3rd decimal place from the empirical answer
from my bisection search below.

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

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

I don't think that's a proof, since R level code can call C functions, and there are an awful lot of callable functions in R, but I don't have a counter-example.

Duncan Murdoch

>

>> But in answer to your question: a bisection search is what I'd use:

>

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 23:02:19 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 - 19:30:04 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.
*