From: Duncan Murdoch <murdoch_at_stats.uwo.ca>

Date: Mon 11 Sep 2006 - 13:40:28 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:
*

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

*>>
*

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

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

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:47:43 2006

Date: Mon 11 Sep 2006 - 13:40:28 GMT

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

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

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

> > I wasn't going into that much detail: just what accuracy are we looking > for here? .Machine$double.neg.eps isn't that accurate (as its help page > says). >

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

> > Oh, it is proof: the storage is declared as C double, and extended > precision double only exists on the chip. Since R has to look up 'x' > whenever it sees the symbol, it looks at the stored value, not on a > register. A compiler could do better, but the current code cannot.

It's a proof of something, but "apparently equivalent formulations" is a vague requirement. I would guess that if I did come up with a counter-example I would have to push the bounds a bit, and it would be questionable whether the formulations were equivalent.

For example, it would clearly be outside the bounds for me to write a function which did the entire computation in C (which would likely give a different answer than R gives). But what if some package already contains that function, and I just call it from R?

Duncan Murdoch

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:47:43 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 - 14: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.
*